Transit Probabilities in Secularly Evolving Planetary Systems

TL;DR

This paper develops a semi-analytical method to estimate transit probabilities in planetary systems influenced by inclined outer planets, helping to explain the Kepler Dichotomy and constrain outer planet inclinations.

Contribution

It introduces a new semi-analytical approach to assess how outer inclined planets affect inner planet transit probabilities and applies this to Kepler systems and the Kepler Dichotomy.

Findings

Outer inclined planets can significantly increase inner planet transit probabilities.

The Kepler Dichotomy can be explained by a combination of system-intrinsic and dynamical inclination effects.

Distribution of mutual inclinations constrains the fraction of systems with outer perturbers.

Abstract

This paper considers whether the population of known transiting exoplanets provides evidence for additional outer planets on inclined orbits, due to the perturbing effect of such planets on the orbits of inner planets. As such, we develop a semi-analytical method for calculating the probability that two mutually inclined planets are observed to transit. We subsequently derive a simplified analytical form to describe how the mutual inclination between two planets evolves due to secular interactions with a wide orbit inclined planet and use this to determine the mean probability that the two inner planets are observed to transit. From application to Kepler-48 and HD-106315 we constrain the inclinations of the outer planets in these systems (known from RV). We also apply this work to the so called Kepler Dichotomy, which describes the excess of single transiting systems observed by Kepler. We find 3 different ways of explaining this dichotomy: some systems could be inherently single, some multi-planet systems could have inherently large mutual inclinations, while some multi-planet systems could cyclically attain large mutual inclinations through interaction with an inclined outer planet. We show how the different mechanisms can be combined to fit the observed populations of Kepler systems with one and two transiting planets. We also show how the distribution of mutual inclinations of transiting two planet systems constrains the fraction of two planet systems that have perturbing outer planets, since such systems should be preferentially discovered by Kepler when the inner planets are coplanar due to an increased transit probability.

Equations

$$S={\int }_0^{{\theta }_0}{\int }_0^{{\varphi }_0}\sin {\theta }^{\prime }d{\theta }^{\prime }d{\varphi }^{\prime }={\left[1-\cos {\theta }^{\prime }\right]}_0^{{\theta }_0}{\left[{\varphi }^{\prime }\right]}_0^{{\varphi }_0},$$

$$-\sin \Delta i\sin \theta \sin \varphi +\cos \Delta i\cos \theta =0.$$

$$-\sin \Delta i\sin {\theta }_1\sin {\varphi }_1+\cos \Delta i\cos {\theta }_1=R_{\star }/a,$$

$$-\sin \Delta i\sin {\theta }_2\sin {\varphi }_2+\cos \Delta i\cos {\theta }_2=-R_{\star }/a.$$

$$\sin I_1=-{\kappa }_2(1-{\kappa }_1^2{)}^{1/2}+{\kappa }_1(1-{\kappa }_2^2{)}^{1/2},$$

$$\sin I_2={\kappa }_2(1-{\kappa }_1^2{)}^{1/2}+{\kappa }_1(1-{\kappa }_2^2{)}^{1/2}.$$

$$P=\frac{2R_{\star }^2}{\pi a_1{a}_2\sin \Delta i}.$$

$$\dot{\mathbf{y}}=i\mathbf{B}\mathbf{y},$$

$$\begin{array}{rl}& B_{jk}=\frac{1}{4}{n}_j\left(\frac{m_k}{M_{\star }+m_j}\right){\alpha }_{jk}{\tilde{\alpha }}_{jk}{b}_{3/2}^{(1)}({\alpha }_{jk})(j\ne k)\\ & B_{jj}=-\sum _{k=1,j\ne k}^N{B}_{jk},\end{array}$$

$$b_s^{(\nu )}(\alpha )=\frac{1}{\pi }{\int }_0^{2\pi }\frac{\cos (\nu x)dx}{(1-2\alpha \cos (x)+{\alpha }^2{)}^s}\alpha <1.$$

$$y_j(t)=\sum _{k=1}^N{I}_{jk}{e}^{i(f_{k}t+{\gamma }_k)},$$

$$i_3=\arctan \left[\frac{(L_1+L_2)\sin \Delta i}{L_3+(L_1+L_2)\cos \Delta i}\right]$$

$$\begin{array}{rl}& y_1=I_{11}{e}^{i(f_{1}t+{\gamma }_1)}+I_{12}{e}^{i(f_{2}t+{\gamma }_2)}\\ & y_2=I_{21}{e}^{i(f_{1}t+{\gamma }_1)}+I_{22}{e}^{i(f_{2}t+{\gamma }_2)},\end{array}$$

$$y_1-y_2=(I_{11}-I_{21})e^{i(f_{1}t+{\gamma }_1)}+(I_{12}-I_{22})e^{i(f_{2}t+{\gamma }_2)}.$$

$$y_1-y_2=(I_{12}-I_{22})\left(e^{i(f_{1}t+\pi )}+e^{i{f}_{2}t}\right).$$

$$y_1-y_2=\Delta iK\left[e^{i(f_{1}t+\pi )}+e^{i{f}_{2}t}\right],$$

$$y_1-y_2\approx \Delta i{K}_{\mathrm{simp}}\left[e^{i(f_{1}t+\pi )}+e^{i{f}_{2}t}\right],$$

$$K_{\mathrm{simp}}=\frac{3m_3{a}_2^{7/2}}{m_2{a}_1^{1/2}{a}_3^3}\frac{1}{b_{3/2}^1\left(\frac{a_1}{a_2}\right)\left(1+(L_1/L_2)\right)}.$$

$$\begin{array}{rl}& f_1\approx -\frac{\pi m_2{a}_1^{1/2}}{2M_{\star }^{1/2}{a}_2^2}{b}_{3/2}^1\left(\frac{a_1}{a_2}\right)\left(1+L_1/L_2\right),\\ & f_2\approx -\frac{3\pi m_3{a}_2^{3/2}}{2M_{\star }^{1/2}{a}_3^3}\frac{1}{1+L_1/L_2}.\end{array}$$

$$\begin{array}{rl}\max |\Delta i_{12}|& \approx 2\Delta i{K}_{\mathrm{simp}}\text{for }{K}_{\mathrm{simp}}<1,\\ & \approx 2\Delta i\text{otherwise. }\end{array}$$

$$\begin{array}{rl}⟨P⟩& \approx R_{\star }/a_2\text{for }\max |\Delta i_{12}|<R_{\star }/a_1-R_{\star }/a_2\\ & \approx \frac{2R_{\star }^2}{\pi a_1{a}_2\sin (\max |\Delta i_{12}|)}\text{otherwise, }\end{array}$$

$$\begin{array}{rl}\Delta i_{\mathrm{crit}}& \approx \frac{R_{\star }/a_1-R_{\star }/a_2}{2K_{\mathrm{simp}}}\text{for }{K}_{\mathrm{simp}}<1\\ & \approx \frac{R_{\star }/a_1-R_{\star }/a_2}{2}\text{otherwise, }\end{array}$$

$$\alpha =\left(\frac{m_{\mathrm{p}}}{M_{\star }}\right)\left(\frac{a_{\mathrm{p}}}{1\mathrm{au}}\right){\left(\frac{d}{1\mathrm{pc}}\right)}^{-1}\mathrm{arcsec},$$

$$i_{\mathrm{e}}=\arctan \left(\frac{\sin \Delta i\left(\sum _{n=1}^3{L}_n\right)}{L_{\mathrm{e}}+\cos \Delta i\left(\sum _{n=1}^3{L}_n\right)}\right),$$

$$K_{\mathrm{simp}}=\frac{3m_{\mathrm{e}}{a}_{\mathrm{d}}^{7/2}}{m_{\mathrm{d}}{a}_{\mathrm{bc}}^{1/2}{a}_{\mathrm{e}}^3}\frac{1}{b_{3/2}^1\left(\frac{a_{\mathrm{bc}}}{a_{\mathrm{d}}}\right)\left(1+(L_{\mathrm{bc}}/L_{\mathrm{d}})\right)},$$

$$\begin{array}{r}{P}_{\mathrm{sing},1}=\frac{R_{\star }}{a_1}-P\\ P_{\mathrm{sing},2}=\frac{R_{\star }}{a_2}-P,\end{array}$$

$$R_H={\left(\frac{m_2+m_3}{3M_{\star }}\right)}^{1/3}\left(\frac{a_2+a_3}{2}\right),$$

$$N_{\mathrm{sing},\mathrm{Kep}}=N_{\mathrm{sing},\mathrm{inh}}+N_{\mathrm{sing},\Delta i_{12}}+N_{\mathrm{sing},\mathrm{planet}}.$$

$$\begin{array}{r}{N}_{\mathrm{sing},\mathrm{Kep}}=N_{\mathrm{sing},\mathrm{inh}}+f(N_{\mathrm{sing},\mathrm{Ndoub}=341}{)}_{\Delta i_{12}}\\ +(1-f)(N_{\mathrm{sing},\mathrm{Ndoub}=341}{)}_{\mathrm{planet}},\end{array}$$

$${\left(\frac{N_{\mathrm{mod}}}{N_{\mathrm{Kep}}}\right)}_{\mathrm{tot}}={\left(\frac{N_{\mathrm{mod}}}{N_{\mathrm{Kep}}}\right)}_{\mathrm{inh}}+f{\left(\frac{N_{\mathrm{mod}}}{N_{\mathrm{Kep}}}\right)}_{\Delta i_{12}}+(1-f){\left(\frac{N_{\mathrm{mod}}}{N_{\mathrm{Kep}}}\right)}_{\mathrm{planet}}.$$

Full text

††pubyear: 2017††pagerange: Transit Probabilities in Secularly Evolving Planetary Systems–Transit Probabilities in Secularly Evolving Planetary Systems

1 Introduction

Over the past 20 years the number of exoplanet detections has soared most notably due to contributions from the Kepler space telescope (Kepler herein). As of November 2016 Kepler has detected 3414 confirmed planets, with 575 existing in multi-planet systems (exoplanet.eu; Schneider et al. (2011)). Planet multiplicity provides information on the underlying architecture of planetary systems, such as expected orbital spacing, mutual inclinations and size distributions. For the multi-planet systems observed by Kepler, super Earth/mini Neptune type objects on tightly packed orbits inside of $\sim$200 days are common (Lissauer et al. (2011); Lissauer et al. (2014); Morton et al. (2016)). Moreover such systems are observed to have small inclination dispersions of $\lesssim$5∘ (Lissauer et al. (2011); Fang & Margot (2012); Figueira et al. (2012); Tremaine & Dong (2012); Marmier et al. (2013); Fabrycky et al. (2014)).

How representative Kepler multi-planet systems are of a common underlying planetary architecture however is impeded by Kepler preferentially detecting objects which orbit closest to the host star. To generalise Kepler systems to an underlying population, it is therefore necessary to account for the inherent probability that transiting systems are observed. Taking into account such probabilities, there appears to be an over-abundance of planetary systems with a single transiting planet (Lissauer et al. (2011); Youdin (2011); Johansen et al. (2012); Ballard & Johnson (2016)). This is commonly referred to as the ’Kepler Dichotomy’.

It is currently not known what causes this excess. Statistical and Spitzer confirmation studies all suggest that the false positive rate for single transiting objects with Rp $<$4R⊕ is low at $\lesssim$15% (Morton & Johnson (2011); Fressin et al. (2013); Coughlin et al. (2014); Désert et al. (2015)). Perhaps then, there are populations of inherently single planet systems in addition to multi-planet systems which are closely packed and have small inclination dispersions. However there may also be a population of multi-planet systems where the mutual inclination dispersion is large, such that only a single planet is observed to transit.

The presence of an outer planetary companion may drive this potential large spread in mutual inclinations. Recent N-body simulations show that the presence of a wide orbit planet in multi-planet systems can decrease the number of inner planets that are observed to transit, either through dynamical instability or inclination excitation (Mustill et al. (2016); Hansen (2017)). Beyond a few au, planetary transit probabilities drop to negligible values. It is possible therefore that additional wide orbit planets could indeed exist in multi-planet systems observed by Kepler. Giant planets at a few au have been detected around stars in the general stellar population by a number of radial velocity (RV) surveys (Marmier et al. (2013); Rowan et al. (2016); Wittenmyer et al. (2016); Bryan et al. (2016)), with suggested occurrence rates ranging from $\sim 10-50\%$ (Cumming et al. (2008); Mayor et al. (2011); Bryan et al. (2016)). Moreover, indirect evidence of undetected giant planets has also been suggested through apsidal alignment of inner RV detected planets (Dawson & Chiang (2014)). As RV studies are largely insensitive to planetary inclinations, it is possible that such wide orbit planets could be on mutually inclined orbits, which may arise from a warp in the disc (Fragner & Nelson (2010)) or due to an excitation by a stellar flyby (Zakamska & Tremaine (2004); Malmberg et al. (2011)).

Calculating transit probabilities of multi-planet systems is complex, often requiring computationally exhaustive numerical methods such as Monte Carlo techniques (e.g. Lissauer et al. (2011); Johansen et al. (2012); Becker & Adams (2016); Mustill et al. (2016); Hansen (2017)). However analytical methods can offer a significantly more efficient route for this calculation and allows for coupling with other fundamental analytical theory, such as for the expected dynamical evolution of the system from inter-planet interactions. Despite this however, analytical investigations into the transit probabilities of multi-planet systems for this purpose are relatively sparse (e.g. Ragozzine & Holman (2010); Brakensiek & Ragozzine (2016)). Recently Brakensiek & Ragozzine (2016) showed how differential geometry techniques can be used to calculate multi-planet transit probabilities by mapping transits onto a celestial sphere. In this paper we perform a similar analysis, however we focus on regions where pairs of planets can be observed to transit. We also give an explicit analytical form using simple vector relations to describe the boundaries of such transit regions.

The multi-planet systems observed by Kepler appear to be mostly stable on long timescales (Lissauer et al. (2011); Pu & Wu (2015)). Dynamical interactions with a potential outer planet on an inclined orbit would therefore be expected to occur on secular timescales. Recent analytical work by Lai & Pu (2017) suggests that such interactions can lead to large mutual inclinations in an inner planetary system, assuming that the direction of the angular momentum vector of the outer planet is fixed. We build on this work by deriving analytical relations for the mutual inclination that can be induced in an inner planetary system by a general planetary companion. We then simplify this result specifically for when the companion is on a wide orbit. Combining this result with our robust analytical treatment of transit probabilities, we can then derive a simple relation describing how the presence of an outer planetary companion affects the transit probability of an inner system due to long term interactions.

We also complement recent N-body simulations of Kepler-like systems interacting with an inclined outer planetary companion shown in Mustill et al. (2016) and Hansen (2017) by using our robust treatment of transit probabilities to consider whether an outer planet with a range of masses, semi-major axes and inclinations can reduce an underlying population of Kepler double transiting systems enough to recover the observed number of single transiting systems through long term interactions only. We also investigate whether the presence of specific wide orbit planets in multi-planet systems preferentially predicts single transiting planets with a given distribution of radii and semi-major axes.

In §2 we overview our semi-analytical method for calculating the transit probability of two mutually inclined planets. In §3 we derive a simplified form to describe the evolution of the mutual inclination between two planets due to presence of an outer planetary companion. We show how this mutual inclination affects the transit probability of the two inner planets in §4. In §5 we apply this work to Kepler-56, Kepler-68, HD 106315 and Kepler-48 to place constraints on the inclination of the outer planets in these systems. In §6 we investigate whether a wide orbit planet in Kepler systems can decrease the number of observed two planet transiting systems enough to recover the observed abundances of single transiting systems. We finally discuss this work in §7 and conclude in §8.

2 Semi-analytical Transit Probability

A planet on a circular orbit with a semi-major axis $a$ and radius $R_{\mathrm{p}}$ subtends a band of shadow across the celestial sphere due its orbital motion. We refer to this band of shadow as the transit region (Ragozzine & Holman (2010); Brakensiek & Ragozzine (2016)). The probability that an observer will view an individual transit event of this planet, assuming that the system is viewed for long enough, is equal to the number of viewing vectors that intersect the transit region, divided by the total number of possible viewing vectors. Perhaps more intuitively, this is equivalent to the surface area of the transit region divided by the total surface area of the celestial sphere.

To calculate the area of a transit region on the celestial sphere first consider that the area of a given surface element ($S$) on a unit sphere is equal to

[TABLE]

where $\theta^{\prime}$ is the polar angle and $\phi^{\prime}$ is the azimuthal angle. A given area on the celestial sphere can therefore be represented on a 2d plane of $1-\cos\theta^{\prime}$ vs. $\phi^{\prime}$, from 0 $\rightarrow$ 2 and 0 $\rightarrow 2\pi$ respectively, such that the 2d plane has a total surface area of $4\pi$. Below we show how the boundaries of a given transit region traverses this 2d plane. This allows for the area contained within these boundaries and therefore the associated transit probability to be calculated.

2.1 Single Planet Case

Consider some fixed reference plane where [$\hat{\mathrm{\mathbf{X}}},\hat{\mathrm{\mathbf{Y}}}$] define a pair of orthogonal directions in this plane, and $\hat{\mathrm{\mathbf{Z}}}$ defines a direction orthogonal to this plane as shown in Figure 1. The fixed reference frame in Figure 1 is assumed to be centred on a host star with radius $R_{\star}$. The line of sight of an observer is considered to be randomly oriented over the surface of a celestial sphere with respect to this fixed reference plane. Now consider that the orbital plane of a planet with a semi-major axis $a$ and radius $R_{\mathrm{p}}$, is inclined to the fixed reference plane by $\Delta i$, with the intersection between the two planes occurring along the $\hat{\mathrm{\mathbf{X}}}$ direction. The direction of the normal of the orbital plane is given by $\hat{\mathrm{\mathbf{n}}}$. The position of a planet in the orbital plane is defined by the direction $\hat{\mathrm{\mathbf{r}}}$ which makes the angles $\theta$ and $\phi$ with the $\hat{\mathrm{\mathbf{Z}}}$ and $\hat{\mathrm{\mathbf{X}}}$ directions respectively. Hence $\hat{\mathrm{\mathbf{r}}}$ traces the centre of the transit region with respect to the fixed reference plane. As $\hat{\mathrm{\mathbf{n}}}\cdot\hat{\mathrm{\mathbf{r}}}=0$, where $\hat{\mathrm{\mathbf{n}}}=[0,\sin\Delta i,\cos\Delta i]$ and $\hat{\mathrm{\mathbf{r}}}=[\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta]$ it follows that

[TABLE]

Hence eq. (2) defines how the centre of a transit region inclined to a fixed reference plane by $\Delta i$ traverses a celestial sphere. This is shown by the dashed lines in Figure 2 for different values of $\Delta i$, where the surface area of the celestial sphere is shown on a 2d plane defined by eq. (1). We note that at the special case where $\Delta i$ = 90∘, $\phi$ can only take values of 0 or $\pi$.

Similarly the directions that define the boundaries of the transit region can be given by $\hat{\mathrm{\mathbf{r}}}_{1}$ and $\hat{\mathrm{\mathbf{r}}}_{2}$ which makes the angles $\theta_{1},\theta_{2}$ and $\phi_{1},\phi_{2}$ with the $\hat{\mathrm{\mathbf{Z}}}$ and $\hat{\mathrm{\mathbf{X}}}$ directions respectively, shown in Figure 1. The boundaries of the transit region also subtend an angle $\pm\theta_{\mathrm{sub}}$ from the orbital plane where $\sin\theta_{\mathrm{sub}}=R_{\star}/a$ assuming $R_{\star}\gg R_{\mathrm{p}}$ (Borucki & Summers (1984)). As $\hat{\mathrm{\mathbf{r}}}_{1}=[\sin\theta_{1}\cos\phi_{1},\sin\theta_{1}\sin\phi_{1},\cos\theta_{1}]$, $\hat{\mathrm{\mathbf{r}}}_{2}=[\sin\theta_{2}\cos\phi_{2},\sin\theta_{2}\sin\phi_{2},\cos\theta_{2}]$ and $\hat{\mathrm{\mathbf{n}}}\cdot\hat{\mathrm{\mathbf{r}}}_{1}=R_{\star}/a$ and $\hat{\mathrm{\mathbf{n}}}\cdot\hat{\mathrm{\mathbf{r}}}_{2}=-R_{\star}/a$, it follows that

[TABLE]

[TABLE]

Hence eq. (3) and eq. (4) describe how the lower and upper boundaries of the transit region for a planet inclined to a fixed reference plane by $\Delta i$ traverse a celestial sphere. The solid lines in Figure 2 show these boundaries for different values of $\Delta i$, where $R_{\star}/a$ = 0.25. This value of $R_{\star}/a$ might be considered to be unrealistically large and is used for demonstration purposes only. In Appendix A we further discuss how the values of ($\theta_{1},\phi_{1}$) and ($\theta_{2},\phi_{2}$) in eq. (3) and eq. (4) respectively would be expected to change as $\Delta i$ is increased from $\Delta i=0\rightarrow 90^{\circ}$.

An integration between the upper and lower boundaries of a transit region divided by the total surface area of the celestial sphere gives the associated single transit probability of the planet ($R_{\star}/a$, Borucki & Summers (1984)). All of the transit regions shown in Figure 2 for different $\Delta i$ therefore contain identical areas and hence have identical single transit probabilities equal to 0.25. We note that if the planet has a non-negligible radius then the single transit probability becomes $(R\pm R_{\mathrm{p}})/a$ for grazing and full transits respectively. Throughout this work however we assume that $R_{\mathrm{p}}\ll R_{\star}$.

2.2 Two Planet Case

Consider now a system containing two planets, both of which are on circular orbits with semi-major axes and radii of $a_{1}$, $a_{2}$ and $R_{\mathrm{{p_{1}}}},R_{\mathrm{{p_{2}}}}$ respectively, where $a_{1}<a_{2}$ and the orbital planes are mutually inclined by $\Delta i$ (we give an exact definition for mutual inclination in §3). The probability that a randomly oriented observer will view both planets to transit (assuming the system is observed for long enough) is equal to the overlap area between the transit regions of both planets, divided by the total area of the celestial sphere. We refer to this probability as the double transit probability.

Therefore, using eq. (3) and eq. (4) to find where the boundaries of the transit regions of each planet intersect, an outline of the overlap between the transit regions can be determined. The area of this overlap can subsequently be calculated by an appropriate integration, which when divided by 4$\pi$ gives the double transit probability. How the double transit probability changes as a function of $\Delta i$ is shown by the blue line in Figure 3, for when $R_{\star}/a_{1}=0.2$ and $R_{\star}/a_{2}=0.1$. We note that this result is consistent regardless of the choice of reference plane and the orientation of the orbital planes of both planets with respect to this reference plane (see Ragozzine & Holman (2010) for a further discussion). That is, the double transit probability depends on the mutual inclination between the two planets only (in addition to the physical size of the respective transit regions).

Depending on the value of $\Delta i$, the double transit probability ($P$ herein) can be split into three regimes (also discussed in Ragozzine & Holman (2010); Brakensiek & Ragozzine (2016)).

(1) For low values of $\Delta i$, the transit region of the outer planet is enclosed within that of the inner planet. The double transit probability is therefore equal to $R_{\star}/a_{2}$.

(2) $\Delta i$ is large enough that the transit region of one planet is no longer fully enclosed inside the other, however there is still partial overlap for all azimuthal angles on the celestial sphere. The transition to this regime occurs for a value of $\Delta i=I_{1}$, which causes $\theta_{1}$ in eq. (3) for both planets to be equal at $\phi_{1}=\pi/2$. Evaluating eq. (3) at this point gives

[TABLE]

where $\kappa_{1}=R_{\star}/a_{1}$ and $\kappa_{2}=R_{\star}/a_{2}$ for simplicity. We note that determining the overlap area of the two transit regions with an exact analytical expression in this regime is difficult and is commonly calculated by Monte Carlo techniques (e.g. Ragozzine & Holman (2010); Johansen et al. (2012); Becker & Adams (2016); Mustill et al. (2016); Hansen (2017)).

(3) For large $\Delta i$, the transit regions only overlap at the intersection of the two orbital planes. The transition to this regime occurs when $\Delta i=I_{2}$, where $\theta_{1}$ for the inner planet is equal to $\theta_{2}$ for the outer planet at $\phi_{1}=\phi_{2}=\pi/2$. Evaluating eq. (3) and (4) here gives

[TABLE]

The values of $I_{1}$ and $I_{2}$ are shown by the green and red lines respectively in Figure 3. If it is assumed that the transit region overlap in regime 3 can be represented as a 2d parallelogram, Ragozzine & Holman (2010) showed the double transit probability can be approximated by111For greater accuracy, we include a 2/$\pi$ factor here that is not included in Ragozzine & Holman (2010).

[TABLE]

For large $\Delta i$ therefore, the double transit probability predicted by eq. (7) tends to a value of 2$R_{\star}^{2}/\pi a_{1}a_{2}$. We show eq. (7) as the black dashed line in Figure 3. We note that in Ragozzine & Holman (2010) it was assumed that the double transit probability transitions straight from regime (1) to (3) at $\Delta i=\arcsin\left(\frac{2}{\pi}\cdot\mathrm{min}(R_{\star}/a_{1},R_{\star}/a_{2})\right)$1.

For $\Delta i>I_{2}$ our method predicts a double transit probability that agrees well with the analytical estimate from Ragozzine & Holman (2010). However there is a clear discrepancy for $I_{1}<\Delta i<I_{2}$, for when there is partial overlap between the transit regions at all azimuthal angles. This highlights the need for semi-analytical methods like the one suggested here over purely analytical relations, to robustly calculate double transit probabilities at all values of $\Delta i$. We note that our method also agrees well with the Monte Carlo treatment of double transit probabilities shown in Ragozzine & Holman (2010).

Calculating transit probabilities using the method outlined here is significantly more computationally efficient than equivalent Monte Carlo methods, as it is only necessary to solve combinations of eq. (3) and (4) for different planets to find where transit regions overlap. From integrating around this overlap, the associated double transit probability is also exact and not subject to Monte Carlo noise effects from under-sampling the total number of line of sight vectors.

3 Secular Interactions

3.1 N planet system

Consider a system of $N$ secularly interacting planets in which planet $j$ has a semi-major axis $a_{j}$ and mass $m_{j}$. The inclination and longitude of ascending node of planet $j$ are given by $I_{j}$ and $\Omega_{j}$ respectively, and can be combined into the associated complex inclination $y_{j}=I_{j}e^{i\Omega_{j}}$. Assuming that the vector involving all the planet’s orbital planes is given by $\bm{y}$ = [$y_{1},y_{2},..y_{N}$], the evolution of complex inclinations in the low inclination and eccentricity limit can be given by Laplace - Lagrange theory in the form

[TABLE]

(Murray & Dermott (1999)) where $\mathbf{B}$ is a matrix with elements given by

[TABLE]

$j$ and $k$ are integers associated with each planet, $M_{\star}$ and $m_{i}$ are the masses of the star and planet $i$, $n_{j}$ is the mean motion of planet $j$ where $n_{j}^{2}a_{j}^{3}=G(M_{\star}+m_{j})$, $\alpha_{jk}=\tilde{\alpha}_{jk}$ = $a_{j}/a_{i}$ for $a_{j}<a_{k}$ and $\alpha_{jk}=a_{k}/a_{j}$ and $\tilde{\alpha}_{jk}=1$ otherwise, and $b^{(1)}_{3/2}(\alpha_{jk})$ corresponds to a Laplace coefficient given by

[TABLE]

Eq. (8) can be solved to show that the evolution of $\bm{y}$ is given by a superposition of eigenmodes associated with each eigenfrequency $f_{i}$ of the matrix $\mathbf{B}$

[TABLE]

where Ijk are the eigenvectors of $\mathbf{B}$ scaled to initial boundary conditions and $\gamma_{k}$ is an initial phase term. If it is assumed that all objects are spherically symmetric, additional terms in the diagonal elements of $\mathbf{B}$ in eq. (9) (e.g. stellar oblateness) need not be included. A choice of reference frame for the inclination also becomes arbitrary, leading to one of the eigenfrequencies equalling zero (c.f. Murray & Dermott (1999)). It is only meaningful therefore to describe a mutual inclination between pairs of planets, with the invariable plane commonly being chosen as a reference plane. The invariable plane is defined as being perpendicular to the total angular momentum vector of a system. The mutual inclination is then the angle between individual angular momentum vectors of a pair of planets. The inclination solution described by eq. (11) also becomes simplified when the invariable plane is taken as a reference plane as the eigenvector associated with the zero value eigenfrequency is also equal to zero.

3.2 Two planet system with an inclined companion

Consider the same general two planet system from §2.2. Assume that the two planets are initially coplanar. Consider now a third planet on an external circular orbit, with a mass and semi-major axis of $m_{3}$ and $a_{3}$ respectively such that $a_{3}>a_{2}$. The orbital plane of this external planet is initially mutually inclined to the inner planets by $\Delta i$. We assume that each of the planets interact through secular interactions only and that inclinations and eccentricities remain small, allowing for application of Laplace - Lagrange theory. Assuming that the invariable plane is taken as a fixed reference plane, the initial inclination of the third planet $i_{3}$ is given by

[TABLE]

where $L_{j}=m_{j}a_{j}^{1/2}$ and is proportional to the angular momentum in the low eccentricity limit. The initial inclination of the inner planets with respect to the invariable plane is therefore $i_{1}=\Delta i-i_{3}$.

From eq. (11) the complex inclination of each of the inner two planets with respect to the invariable plane evolves in the form of

[TABLE]

where $y_{1}$ and $y_{2}$ are the complex inclinations of the innermost and second innermost planet respectively. The evolution of the mutual inclination between the inner pair of planets is hence given by

[TABLE]

The $t=0$ boundary conditions give $\gamma_{1}=\pi$ and $\gamma_{2}=0$. Also as $y_{1}$($t$ = 0) = $y_{2}$($t$ = 0) = $i_{1}$, it follows from eq. (12) that $I_{11}-I_{21}$ = $I_{12}-I_{22}$. The evolution of the mutual inclination from eq. (13) is therefore is equivalent to

[TABLE]

Hence the evolution of the instantaneous mutual inclination between the inner pair of planets, $\Delta i_{12}=|y_{1}-y_{2}|$, can be calculated if the first and second elements of the eigenvector associated with the $f_{2}$ eigenfrequency are known. In Appendix B we fully solve eq. (8) to give $I_{12}$ and $I_{22}$ in terms of physical variables. Here we simply say that

[TABLE]

where $K$ is dependant on the masses and semi-major axes of the three planets, shown explicitly in Appendix B. We note that the maximum value of $K\approx 1$, implying that the maximum value of the mutual inclination between the inner pair of planets from eq. (14) is twice the initial mutual inclination with the external third planet i.e. max|$\Delta i_{12}$| = 2$\Delta i$. For given values of masses and semi-major axes of the inner pair of planets therefore, the evolution of the mutual inclination between them is dependant on three quantities, $a_{3}$, $m_{3}$ and $\Delta i$.

The left panels of Figure 4 show how max|$\Delta i_{12}$| changes as a function of different combinations of $a_{3}$, $m_{3}$ and $\Delta i$ in eq. (15) for an example system where $a_{1}$, $a_{2}$ = 0.2, 0.5au and $m_{1}$, $m_{2}$ = 10M⊕ respectively. We note that the assumptions of Laplace-Lagrange theory are expected to break down when $\Delta i\gg 20^{\circ}$. Larger inclinations are included for demonstration purposes only. It is evident that as the third planet tends to a limit where it is on a wide orbit, with a low mass and low initial mutual inclination, the maximum mutual inclination between the inner pair of planets becomes small as one might expect.

3.3 Companion wide orbit approximation

In §4 we look to investigate how the evolving mutual inclination between the inner pair of planets affects the associated double transit probability, for the specific case where the external third planet is assumed to be on a wide orbit. For $a_{3}\gg a_{1},a_{2}$, certain individual and combinations of $\mathbf{B}$ matrix elements from eq. (9) become small and we find that eq. (15) can be simplified to

[TABLE]

where

[TABLE]

Here it is assumed that as $a_{3}\gg a_{1},a_{2}$, certain Laplace coefficients from the $\mathbf{B}$ matrix elements can be simplified, specifically $b^{1}_{3/2}(\alpha)\approx 3(\alpha)$ (Murray & Dermott (1999)). Similar simplifications can be made to each of the eigenfrequencies, for which

[TABLE]

As eq. (15) shows that the maximum value of the mutual inclination between the inner pair of planets cannot be larger than twice the initial mutual inclination with the wide orbit planet (max$|\Delta i_{12}|\ngtr 2\Delta i$) we assume that the maximum value of the mutual inclination between the inner two planets predicted by eq. (16) is

[TABLE]

The right panels of Figure 4 show max|$\Delta i_{12}$| predicted by eq. (19) and eq. (17) using the same planet parameters as shown in the left panels. We find that when $a_{3}$ $\gtrsim$1.25au, the simplified form for max|$\Delta i_{12}$| from eq. (19) and eq. (17) agrees with the full Laplace - Lagrange solution to within $\sim 25\%$ for all values of $m_{3}$ and $\Delta i$. For $a_{3}\sim$1au, the simplified form of max|$\Delta i_{12}$| begins to break down and eq. (19) can underestimate max|$\Delta i_{12}$| from the full Laplace - Lagrange solution by up to a factor of 2.

This estimate is similar to the result derived by Lai & Pu (2017), who assumed that the angular momentum vector direction of the outer inclined planet is fixed in time. They find that the maximum mutual inclination that can be induced in an inner pair of planets depends on the strength of the coupling between them (parametrized by $\epsilon_{12}$ in their eq. 12). Assuming inclinations are small we find eq. (19) agrees with the equivalent prediction of max|$\Delta i_{12}$| from Lai & Pu (2017) if $K_{\mathrm{simp}}=\epsilon_{12}$. Indeed, $K_{\mathrm{simp}}$ and $\epsilon_{12}$ are almost identical despite the different derivation techniques (e.g. we derive the full Laplace- Lagrange solution and then simplified assuming $a_{3}\gg a_{1},a_{2}$), apart from $K_{\mathrm{simp}}$ contains an additional factor of $a_{1}^{1/2}a_{2}^{3/2}$ whereas $\epsilon_{12}$ contains a factor of ($a_{2}^{2}-a_{1}^{2}$). By considering different combinations of $a_{1}$ and $a_{2}$ and comparing to the value of max|$\Delta i_{12}$| given by the full solution in Appendix B, we find that neither eq. (19) and (17) or the equivalent equation from Lai & Pu (2017) is favoured as a more accurate approximation, since which is closer to the full solution depends on the exact parameters.

4 Combining Transit Probabilities with Secular Theory

Considering two inner, initially coplanar planets and an outer inclined planetary companion, we combine the analysis of transit probabilities from §2 with secular interactions from §3 in two main ways. First in §4.1, we assume that the outer planet is not necessarily on a wide orbit. The evolution of the mutual inclination between the inner planets is therefore assumed to be given by the full Laplace-Lagrange solution derived in eq. (15). The double transit probability of the inner two planets during this evolution is then calculated through the method outlined in §2. This provides the most accurate prediction for how the double transit probability of two inner planets evolves (in the low inclination limit) considering a given outer planetary companion. We make use of this method for a detailed discussion of how an outer planet affects an inner population of Kepler systems in §6.

Second, in §4.2 we assume that the outer planetary companion is on a significantly wide orbit. The evolution of the mutual inclination between the inner two planets is therefore given by eq. (16) and eq. (17). Here we look to give a simple analytical form to describe the double transit probability of two inner planets, due to secular interactions with a given outer planetary companion. We make use therefore of simple analytical relations such as eq. (7) to describe double transit probabilities. Comparing with the work in §4.1 allows for the accuracy of these approximations to be judged. We demonstrate in §5 how simple constraints can be placed on the inclination of an outer companion in specific systems using this method.

4.1 Two planet system with an inclined companion

From Figure 3 it is clear that if the amplitude of the mutual inclination between the inner two planets is large, then the associated double transit probability, $P$, will only be at a maximum value for a small proportion of the secular evolution. The presence of an outer inclined planet may therefore result in a significant reduction in the mean double transit probability $\langle P\rangle$ on long timescales. Figure 5 shows how both the mutual inclination and the double transit probability evolve with time for two inner planets from Figure 4, which are perturbed by an outer planetary companion with a semi-major axis, mass and inclination of $a_{3}$ = 2au and $m_{3}$ = 1MJ and $\Delta i$ = 5∘ respectively. Indeed, $P$ is only at a maximum value for a small proportion of the secular evolution leading to a significant reduction in $\langle P\rangle$ compared with if the outer planet were not present.

Furthermore, the left panels of Figure 6 show how $\langle P\rangle$ changes due to perturbations from an outer planet with the same range of parameters considered in Figure 4. As one may expect, through comparing the left panels of Figures 4 and 6, an outer planet which induces a large value of max|$\Delta i_{12}$| also causes a significant reduction in the mean double transit probability of the inner two planets and vice versa for small values of max|$\Delta i_{12}$|.

The left panels of Figure 6 also suggest a clear boundary of $a_{3}$, $m_{3}$ and $\Delta i$, above which the outer planet causes $\langle P\rangle$ to be significantly reduced and below which $\langle P\rangle$ is unchanged. From Figure 3, the double transit probability of the two inner planets can be considered to be significantly reduced when $\Delta i_{12}>I_{1}$, where $I_{1}$ is given by eq. (5). We assume therefore that the boundary where $\langle P\rangle$ is significantly reduced occurs when max|$\Delta i_{12}$| $\approx$ $I_{1}$. The values of $a_{3}$, $m_{3}$ and $\Delta i$ which give this boundary are shown by the black lines in the left panels of Figure 6.

4.2 Companion wide orbit approximation

Considering the simplified evolution of the mutual inclination from eq. (16) and (17) for when $a_{3}\gg a_{1},a_{2}$, here we estimate the value of the mean double transit probability itself. We assume that $\langle P\rangle$ is dominated by the maximum or minimum value of the double transit probability, $P_{\mathrm{max}}$ and $P_{\mathrm{min}}$ respectively, depending on whether max|$\Delta i_{12}$| is greater than $I_{1}$. We assume that $I_{1}\approx R_{\star}/a_{1}-R_{\star}/a_{2}$ from eq. (5) for $R_{\star}/a_{1}$, $R_{\star}/a_{2}$ $\ll$ 1. From Figure 3, the value of $P_{\mathrm{max}}$ = $R_{\star}/a_{2}$, however a value of $P_{\mathrm{min}}$ is more difficult as no specific analytical estimate exists. We therefore assume $P_{\mathrm{min}}$ can be given by the estimate from Ragozzine & Holman (2010) shown by eq. (7). We note that this approximation for $P_{\mathrm{min}}$ would be expected to break down if max|$\Delta i_{12}$| predicts partial overlap between the transit regions of the inner planets for all azimuthal angles (see Figure 3). Assuming that the masses and semi-major axes of all the planets are known, in addition to the inclination of the outer planet and that max|$\Delta i_{12}$| is given by the simplified Laplace - Lagrange solution from eq. (19), $\langle P\rangle$ can be estimated by

[TABLE]

The right panels of Figure 6 show the value of $\langle P\rangle$ predicted by eq. (20), using the same planet parameters as those in the left panel. The black lines are identical to those in the left panels of Figure 6 and are included to aid comparison between both sides of the Figure.

The above assumptions bias the double transit probability toward spending a greater proportion of the secular evolution at $P_{\rm{min}}$. As such, eq. (20) can under predict $\langle P\rangle$, by a factor of up to 4 when comparing the left and right panels of Figure 6. We suggest therefore that eq. (20) should be used as a first order approximation of $\langle P\rangle$ only.

5 Application to specific systems

Here we consider real systems observed to have both transiting planets and an additional outer, non-transiting planet. Due to the inherent faintness of Kepler stars, follow up observations to detect non-transiting planets, namely by RV studies, are challenging. Thus the number of systems observed with such architectures are relatively sparse. We consider three of these systems: Kepler-56, Kepler-68 and Kepler-48 in addition to HD 106315. As RV surveys are largely insensitive to planetary inclinations, we apply eq. (20) with eq. (17) to place constraints on the inclination of the non-transiting planets in these systems.

Assume that, as the transiting planets are indeed transiting, the mean double transit probability is at a maximum. Rearranging eq. (20) one finds

[TABLE]

where $\Delta i_{\mathrm{crit}}$ is the inclination of the non-transiting planet required to significantly reduce the mean probability that the inner planets are observed to transit due to secular interactions. We note that eq. (21) assumes that the transiting planets are initially coplanar. However if these planets were initially mutually inclined by a small amount, a smaller secular perturbation from the outer planet would be required to significantly reduce the mean probability that the inner planets are observed to transit. In this case, $i_{\mathrm{crit}}$ from eq. (21) would be reduced.

5.1 Kepler-56

Kepler-56 is a red giant star with a mass and radius of M⊙ = 1.32 $\pm$ 0.13 M⊙ and R⊙ = 4.23 $\pm$ 0.15R⊙ respectively (Huber et al. (2013)), which is observed to host three planets. Interestingly, Kepler-56 represents one of the few red giant stars observed to host a planetary system (Lillo-Box et al. (2014); Ciceri et al. (2015); Quinn et al. (2015); Pepper et al. (2016)). The two inner planets (b, c) are observed to transit with periods of 10.5 and 21.4 days respectively (Borucki et al. (2011); Steffen et al. (2013); Huber et al. (2013); Hadden & Lithwick (2014); Holczer et al. (2016); Morton et al. (2016)) and have masses of 22.1${}^{+3.9}_{-3.6}$M⊕ and 181${}^{+21}_{-19}$M⊕ respectively (Huber et al. (2013)). Keck/HIRES and HARPS-North observations have revealed a non-transiting giant planet (d) with a period of 1002$\pm$5 days and minimum mass of 5.62$\pm$0.38MJ (Huber et al. (2013); Otor et al. (2016)). An interesting quirk of this system is that the transiting planets, while being roughly coplanar, are misaligned to the stellar spin axis by $\sim$40∘ (Huber et al. (2013)). It is unclear if this large obliquity is caused by long term dynamical interactions with a highly inclined companion, such as Kepler-56d, or from the star being inherently tilted to the disk from which the planets formed (Li et al. (2014)).

Applying eq. (21), we find that $i_{\mathrm{crit}}=704^{\circ}$. This unphysically large value means that, regardless of how Kepler-56d is inclined in this system, the mean double transit probability of the inner two transiting planets cannot be significantly reduced. That is, we suggest that the transiting planets in Kepler-56 are not strongly affected by the secular perturbations of Kepler-56d, regardless of its mutual inclination. This is a similar result to that found in Lai & Pu (2017) who also find that the inner planets are strongly coupled against external secular interactions. We therefore cannot place any constraint on the inclination of Kepler-56d using this method. We note however that this does not preclude that the 40∘ misalignment from the stellar spin axis comes from an inclined outer companion, since both inner planets could be inclined together without significant mutual inclination.

5.2 Kepler-68

Kepler-68 is a roughly solar type star with a mass and radius of 1.08$\pm$0.05M⊙ and 1.24$\pm$0.02R⊙ respectively (Gilliland et al. (2013); Marcy et al. (2014)). It hosts two transiting planets (b, c) with periods of 5.4 and 9.6 days respectively (Gilliland et al. (2013); Marcy et al. (2014); Van Eylen & Albrecht (2015); Holczer et al. (2016); Morton et al. (2016)) and fitted masses of 5.97$\pm$1.70 and 2.18$\pm$3.5M⊕ respectively (Marcy et al. (2014)). Keck/HIRES RV follow up of this system detected a non-transiting planet (d) with a period of 625$\pm$16 days with a fitted mass of 267$\pm$16M⊕ (Marcy et al. (2014)).

Applying eq. (21) we find $i_{\mathrm{crit}}=244^{\circ}$. Similar to Kepler-56 therefore, regardless of the mutual inclination of Kepler-68d, the mean double transit probability of the inner two transiting planets cannot be significantly reduced by secular perturbations. We therefore cannot place a constraint on the inclination of Kepler-68d using this method. We note that Kepler-68d can indeed have a large inclination without affecting the overall stability of the system according to a suite of N-body simulations, which suggest that Kepler-68d is inclined by $\Delta i<85^{\circ}$ (Kane (2015)).

5.3 HD 106315

HD 106315 is a bright F dwarf star at a distance $d=107.3\pm 3.9$pc (Gaia Collaboration et al. (2016)) with mass and radius of 1.07$\pm$0.03M⊙ and 1.18$\pm$0.11R⊙ respectively (Morton (2012); Petigura (2015); Crossfield et al. (2017)). Recent K2 observations detect two transiting planets (b, c) with periods of 9.55 and 21.06 days respectively and radii of 2.23${}^{+0.30}_{-0.25}$ and 3.95${}^{+0.42}_{-0.39}R_{\oplus}$ respectively (Crossfield et al. (2017); Rodriguez et al. (2017)). Mass-radius relationships suggest these planets have masses of 8 and 20M⊕ respectively (Weiss et al. (2016); Wolfgang et al. (2016); Crossfield et al. (2017)). Further Keck/HIRES RV observations also indicate the presence of a third outer companion planet (d) with a period of $P_{\mathrm{d}}\gtrsim 80$ days, which has a mass of $m_{\mathrm{d}}\gtrsim$1MJ (Crossfield et al. (2017)). As the exact period of this outer planet is unknown we consider two possibilities where the outer planet has a period of $P_{\mathrm{d}}=80$ days and $P_{\mathrm{d}}=365$ days respectively. Assuming $P_{\mathrm{d}}=80$ days implies a mass of $m_{\mathrm{d}}=1$MJ (Winn et al. (2009); Crossfield et al. (2017)). Applying eq. (21) with this outer planet gives $i_{\mathrm{crit}}=1.1^{\circ}$. This suggests that if the outer planet had a period of $P_{\mathrm{d}}=80$ days, it must have an inclination of $\Delta i\lesssim 1.1^{\circ}$, otherwise the mean probability of observing the inner two planets to transit would be significantly reduced due to the secular interaction. Conversely, if the outer planet is assumed to be further out with $P_{\mathrm{d}}=365$ days, implying a mass of $\sim$7MJ, eq. (21) suggests that $i_{\mathrm{crit}}=2.4^{\circ}$. That is, if the outer planet has a period of $P_{\mathrm{d}}=365$ days, it must have an inclination of $\Delta i\lesssim 2.4^{\circ}$, otherwise the secular interaction would significantly reduce the mean probability that the inner planets are observed to transit.

The mutual inclination of the outer planet might also be constrained through astrometric observations of HD 106315 with ESA’s Gaia mission (Perryman et al. (2001); Casertano et al. (2008); Sozzetti et al. (2014); Perryman et al. (2014); Sahlmann et al. (2015)). The astrometric displacement of the host star due to the presence of a planet is defined by

[TABLE]

with the astrometric signal-to-noise equal to $S/N=\alpha\sqrt{N_{\mathrm{obs}}}/\sigma$, where $N_{\mathrm{obs}}$ is the scheduled number of astrometric measurements ($N_{\mathrm{obs}}$ = 36 for HD 106315222http://gaia.esac.esa.int/gost/) with typical uncertainties of $\sigma=40\mathrm{\mu as}$ (de Bruijne (2012)). If $S/N>20$, the orbital inclination can be constrained to a precision of $<10^{\circ}$ (Sahlmann et al. (2015)). We find that for the example periods and masses considered above for HD 106315d that $S/N<10$. We therefore expect that the inclination of the above examples of HD 106315d cannot be constrained using Gaia astrometry. However if HD 106315d is outside of $\sim$1.3au, (implying a mass of $\gtrsim 12\mathrm{M}_{\mathrm{J}}$) eq. (22) suggests that $S/N>20$ such that the inclination of HD 106315d should be constrained by Gaia astrometry. Further RV follow-up of this system will allow for greater constraints to be placed on the mass and the orbit of HD 106315d, which in turn allow for greater constraints to be placed on the inclination, either through potential astrometry measurements or through our model represented by eq. (21).

5.4 Systems with three transiting planets and a wide orbit companion

Here we generalise the affect a wide orbit planet has on the transit probabilities of three inner transiting planets. Consider Kepler-48 as an example of such a system. Kepler-48 has a mass and radius of M⋆ = 0.88$\pm$0.06M⊙ and R⋆ = 0.89$\pm$0.05R⊙ respectively. It hosts three transiting planets (b,c,d) with periods of 4.78, 9.67 and 42.9 days and fitted masses of 3.94$\pm$2.10, 14.61$\pm$2.30 and 7.93$\pm$4.6M⊕ respectively (Steffen et al. (2013); Marcy et al. (2014); Hadden & Lithwick (2014); Holczer et al. (2016); Morton et al. (2016)). Keck/HIRES RV analysis also detects a non-transiting planet (e) with a period and fitted mass of 982$\pm$8 days and 657$\pm$ 25M⊕ respectively (Marcy et al. (2014)).

Returning to the derivation of the secular interaction in §3, the initial inclination of the non-transiting planet, $i_{e}$, with respect to the invariable plane can be generalised to

[TABLE]

where $L_{\mathrm{e}}=m_{e}a_{e}^{1/2}$ and is proportional to the angular momentum of Kepler-48e in the low eccentricity limit and $L_{n}=m_{n}a_{n}^{1/2}$ for either Kepler-48b, c, or d. The initial inclination of the transiting planets is therefore equal to $\Delta i-i_{\mathrm{e}}$.

As the strength of the secular interaction between planets largely depends on their separation (e.g. eq. (19)) we assume that Kepler-48d will be affected most by perturbations from the non-transiting planet. We demonstrate this in Figure 7, which shows how the mutual inclination between each of the transiting planets evolves assuming Laplace - Lagrange theory (eq. (11)) and that Kepler-48e is initially mutually inclined by $\Delta i=10^{\circ}$. The red line shows the mutual inclination between Kepler-48b and c ($\Delta i_{\mathrm{bc}}$), the blue between b and d ($\Delta i_{\mathrm{bd}}$), and the green between c and d ($\Delta i_{\mathrm{cd}}$). The mutual inclination between Kepler-48b and c is largely unchanged and they remain roughly coplanar. Conversely the mutual inclination between b and d and c and d is significant and roughly equal throughout the secular evolution. It can be assumed for Kepler-48 therefore that the inner two transiting planets are largely unaffected by the secular perturbations of Kepler-48e, but both can become significantly mutually inclined to the outer transiting planet.

As such, we assume that Kepler-48b and c can be treated as a single body whose angular momentum is the sum of Kepler48-b and c, reducing the system to a total of three planets. With this approximation, the evolution of the mutual inclination between Kepler-48b and c with d ($\Delta i_{\mathrm{bc,d}}$) is shown by the dashed black line in Figure 7. It can be seen that this way of treating Kepler-48b and c as a single body gives a good approximation for the evolution of the mutual inclination between Kepler-48b, c with d.

The initial mutual inclination of Kepler-48e which causes a significant reduction in the mean probability of the inner planets transiting, $\Delta i_{\mathrm{crit}}$, can therefore be approximated by eq. (21), where the value of $K_{\mathrm{simp}}$ becomes

[TABLE]

with the subscripts referring to a respective planet and the subscript ’bc’ to the planet which has the same total angular momentum as Kepler-48b and c.

We find that $\Delta i_{\mathrm{crit}}$ = 3.7∘. This suggests therefore that the inclination of Kepler-48e $\Delta i\lesssim 3.7^{\circ}$, otherwise the secular interaction would cause a significant reduction in the mean probability that all three inner planets are observed to transit. Under the simpler assumption that max|$\Delta i_{\mathrm{bc,d}}$| $\lesssim$ $R_{\star}/a_{\mathrm{d}}$, Lai & Pu (2017) also find that the inclination of Kepler-48e, considering secular interactions only, must also be small with $\Delta i\lesssim 2.3^{\circ}$.

6 Application to the Kepler Dichotomy

As discussed in §1, Kepler has observed an excess of single transiting systems which cannot be explained by geometric effects alone, commonly referred to as the Kepler dichotomy (Lissauer et al. (2011); Youdin (2011); Johansen et al. (2012); Ballard & Johnson (2016)). This may suggest that there is a population of inherently single transiting systems in addition to a population of multi-planet systems with small inclination dispersions. However there may also be a population of multi-planet systems where the mutual inclination dispersion is large, increasing the probability that only a single planet is observed to transit. Here we investigate whether both these types of multi-planet systems can significantly contribute to the abundance of systems observed by Kepler to have one and two transiting planets respectively.

The Kepler systems we consider are discussed in §6.1. A method for debiasing Kepler systems to a general population of planetary systems is described in §6.2. We consider the scenario where planets share some inherently fixed mutual inclination in §6.3, before considering when this mutual inclination is evolving due to the presence of an outer inclined planetary companion in §6.4. We note from the outset that we do not consider Kepler systems observed to have more than two planets. Instead we look to explore what effects an outer planet might have on observables of a subset of Kepler like systems, rather than observables of the whole Kepler population. We discuss this assumption further in §7.6.

6.1 Kepler Candidate Sample

We select planet candidates from the cumulative Kepler objects of interest (KOI) table from the NASA exoplanet archive333exoplanetarchive.ipac.caltech.edu, accessed on 13/09/16. The vast majority of the KOIs ($\sim 97\%$) that survive our cuts detailed below, to make it into our final sample are listed as being taken from the most recent Q1-17 DR24 data release. This data release is of particular note as it incorporates an automated processing of all KOIs (Coughlin et al. (2016)).

Out of the initial 8826 KOIs we consider those which orbit solar type stars, with surface temperatures and surface gravities between 4200K $<T<$ 7000K and 4.0 $<$ log($g$) $<4.9$ respectively. This reduces the total number of KOIs to 7446. We also find the total number of unique Kepler stars within this range (discussed in §7) is 164966 from the ’Kepler Stellar data’ table. We next remove false positives, which refer to KOI light curves that are indicative of either an eclipsing binary, having significant contamination from a background eclipsing binary, showing significant stellar variability which mimics a planetary transit or where instrument artefacts have produced a transit like signal (see Coughlin et al. (2014); Rowe et al. (2014); Rowe & Thompson (2015); Seader et al. (2015); Coughlin et al. (2016)). This reduces our sample of KOIs (candidates herein) to 4072 objects. We subsequently remove non planetary-like objects with radii $>$22.4$R_{\oplus}$ (Borucki et al. (2011)), leaving 3757 objects, after which we remove candidates with a SNR $<10$ reducing the possibility that a transit signal is caused by systematic background noise (Morton et al. (2016)), leaving 3327 objects. Finally we remove candidates listed as not having a satisfactory fit to the transit signal (Rowe et al. (2014); Rowe & Thompson (2015)). This gives our final sample of 3255 objects. We note that our choice of cuts means that KOI systems can become reduced in multiplicity. We find that our final sample includes systems which contain 1-6 candidates with $N_{i}$ = (1, 2, 3, 4, 5, 6) = (1951, 341, 117, 43, 15, 4) e.g. 1951 systems with a single candidate, 341 systems with two candidates etc. Herein, we consider the 1951 systems observed by Kepler to have a single transiting planet and the 341 systems observed to have two.

The smoothed distribution of the semi-major axis and planetary radii for the single and double planet transiting systems are shown in Figure 8. Comparing the left and right panels of Figure 8, there are types of planets which are only present in single transiting systems. We briefly discuss these differences here for future reference. Large planets with short periods i.e. Hot Jupiters, are not present in Kepler systems with two transiting planets. Indeed, investigations into the formation processes of Hot Jupiters predict a lack of close companions (Wright et al. (2009); Steffen et al. (2012); Mustill et al. (2015); Huang et al. (2016), see WASP-47 for an exception, Becker et al. (2015); Almenara et al. (2016)). Long period planets are also more abundant in the population of single transiting systems. This may not necessarily indicate that long period planets inherently favour being in single transiting systems, but instead they might be the inner planet of a higher multiplicity system where the outer planets are on too long a period to produce a significant transit signal.

Finally there appears to be an over abundance in the population of single transiting systems for planets with $R_{\mathrm{p}}\lesssim 2R_{\oplus}$ at periods $P<10$ days ($\lesssim 0.03$au) (see Lissauer et al. (2011); Johansen et al. (2012); Steffen & Coughlin (2016); Lopez & Rice (2016)). The formation processes which lead to these types of planets are unclear. It is also unknown if these objects are inherently rocky planets, or are the cores of Neptune sized planets whose envelopes have been irradiated (Dressing & Charbonneau (2015); Rogers (2015); Lopez & Rice (2016)). If these outlying systems are largely ignored, the question remains of whether the remaining planets in single transiting systems are part of the same underlying distribution of higher order planetary systems; i.e. could these single transiting systems contain similar planets which are not observed to transit?

For our dynamical analysis it is not the radii of these planets which is of relevance, rather their masses. We estimate the masses of planets according to the following mass-radius relations. For radii less than 1.5$R_{\oplus}$ we use the rocky planet mass-radius relation from Weiss & Marcy (2014), where density ($\rho_{\mathrm{p}}$) is related to radii ($R_{\mathrm{p}}$) through $\rho_{\mathrm{p}}=2.43+3.39(R_{\mathrm{p}}/R_{\oplus})$gcm*-3*. For radii 1.5 $\leq$ $R_{\mathrm{p}}$ $\leq$ 4$R_{\oplus}$, we use the deterministic version of the probabilistic mass-radius relation for sub-Neptune objects from Wolfgang et al. (2016), where mass ($M_{\mathrm{p}}$) is given by $M_{\mathrm{p}}/M_{\oplus}=2.7(R_{\mathrm{p}}/R_{\oplus})^{1.3}$. Once radii become $R_{\mathrm{p}}\gtrsim$4$R_{\oplus}$ deterministic mass-radius relations become uncertain due to the onset of planetary contraction under self-gravity (see Chen & Kipping (2017)). From the mass-radius relations detailed in Chen & Kipping (2017), we find their ’Neptunian worlds’ deterministic relation of $M_{\mathrm{p}}/M_{\oplus}=(1.23R_{\mathrm{p}}/R_{\oplus})^{1.7}$ gives the most sensible masses for all planets with $R_{\mathrm{p}}>4R_{\oplus}$.

6.2 De-biasing the Kepler population

As previously alluded to, Kepler only observes planetary systems that have their orbital planes aligned with our line of sight. It is therefore sensible to suggest that there is a much larger, underlying population of planetary systems within which only some are observed to transit. We refer to this underlying population of planetary systems as the model population. Conversely, we refer to the population of planetary systems actually observed by Kepler as the Kepler population. We assume that Kepler systems are representative of planetary systems in the model population once geometrical biases have been taken into account.

To construct an underlying model population, our primary goal is for this to predict the correct number and planet parameter distribution seen in the Kepler population for systems with two transiting planets (Figure 8 right). To achieve this we first assume that all stars either have two or zero planets. Any system which hosts two planets is assumed to be identical to one of the 341 double transiting systems observed by Kepler. We assume the abundance of a specific Kepler-like system in the model population is equal to the inverse of the mean of the double transit probability calculated by the method outlined in §2. Systems with inherently low mean double transit probabilities, are therefore probabilistically assumed to be more numerous in the model population. By definition therefore, each unique system in the model population would be expected to be observed with both planets transiting exactly once and so the model population predicts the correct distribution shown in the right panel of Figure 8. We note that a model population generated in this way is similar to the method described in Johansen et al. (2012), albeit with their work predicting the correct number and planet parameter distribution seen in the Kepler population for systems with three transiting planets.

The sum of the inversed mean double transit probabilities of all the 341 double transiting systems gives the total number of planetary systems in the model population. If we assume that all of the two planet systems are coplanar, we find the model population includes 16517 systems (the remaining 148449 systems observed by Kepler are assumed to have no planets).

Each system in the model population can be observed to have a single transiting planet, depending on the viewing angle. The sum of the mean single transit probabilities for each of the 16517 systems in the coplanar model population gives the total number of single transiting planets, $N_{\mathrm{sing}}$, that would be expected to be observed. Here the mean single transit probability for a given system is equal to $R_{\star}/a_{1}-R_{\star}/a_{2}$, where $a_{1}$, $a_{2}$ are the semi-major axes of each planet when $a_{2}>a_{1}$ and $R_{\star}$ is the radius of the host star. We find $N_{\mathrm{sing}}$ = 589, which clearly underestimates the 1951 single transiting systems in the observed Kepler population, by a factor of $\sim 3$. This is the Kepler dichotomy discussed in §1. We show the smoothed distribution of the semi-major axes and planet radii for these 589 predicted single transiting planets in the top left panel of Figure 9, which when compared with the left panel of Figure 8 clearly shows an under-prediction of the single transiting planets observed by Kepler.

6.3 Inherently inclined multi-planet systems

From transit duration variation (TDV) studies, the mutual inclinations of planets in multi-transiting systems are small at $\lesssim 2-3^{\circ}$ (Fang & Margot (2012); Fabrycky et al. (2014)). We note that this mutual inclination also best fits the distribution of impact parameters in the Kepler population. Perhaps then, if two planets are assumed to be inherently mutually inclined by a small amount, this may account for the abundance of single transiting planets in the Kepler population. Consider a fixed mutual inclination $\Delta i_{12}$ between the two planets in each of the 341 double transiting systems. The mean single transit probability for each planet from a given system, $P_{\mathrm{sing,1}}$ and $P_{\mathrm{sing,2}}$ respectively where $P_{\mathrm{sing,1}}>P_{\mathrm{sing,2}}$, is now given by

[TABLE]

where $P$ is the mean double transit probability and $P_{\mathrm{sing,1}}+P_{\mathrm{sing,2}}$ is the total mean single transit probability for this system. As $\Delta i_{12}$ increases, the mean double transit probability decreases (Figure 3). Therefore for a fixed population of double transiting systems considered here, the expected abundance of single transiting systems increases. Figure 10 shows how $N_{\mathrm{sing}}$ increases with $\Delta i_{12}$ for when the number of double transiting systems is kept constant at 341 systems. If $\Delta i_{12}=4.4^{\circ}$, we find $N_{\mathrm{sing}}=1951$, i.e. the number of single transiting planets expected to be observed from the model population is equal to the number in the observed Kepler population. This suggests that mutual inclinations in Kepler systems observed with two planets must be less than 4.4∘, or the number of single planet systems observed by Kepler would be too large relative to the number of doubles.

We show the distribution of the semi-major axes and radii of the expected single transiting planets for when $\Delta i=4.4^{\circ}$ in the top right panel of Figure 9. Comparing with the left panel of Figure 8, there is an over abundance of predicted single transiting planets with radii of $\sim 2.5$R⊕ and semi-major axes of $\sim$0.15au. This is due to the model population compensating for not being able to reproduce all types of single transiting planets in the Kepler population (e.g. Hot Jupiters discussed in §6.1). Herein therefore when discussing how well a model population predicts the Kepler population of single transiting planets we refer to how well the types of planets from each population compare, rather than the total number. That is, we look to find which value of $\Delta i_{12}$ causes the associated version of the top right panel of Figure 9 to be most like the left panel of Figure 8.

We judge the success of this comparison using a modified $\chi^{2}$ minimisation test, in which we simply sum the square of the difference between the number of singles with a given radius and semi-major axis expected from the model population, with that of the observed Kepler population. Varying $\Delta i_{12}$ we therefore look to identify a minimum in the modified $\chi^{2}$ space without caring for the modified $\chi^{2}$ value itself. We show this in Figure 10, with the modified $\chi^{2}$ minimum occurring for $\Delta i_{12}=3.6^{\circ}$. The distribution of the single transiting planets expected from the model population for this mutual inclination is shown in the bottom left panel of Figure 9. Comparing with the left panel of Figure 8, these single transiting planets share a stronger agreement with those in the Kepler population, compared with when the outer planet predicted $N_{\mathrm{sing}}=1951$ (e.g. top right panel of Figure 9). We note that the total number of single transiting planets expected from the model population for $\Delta i_{12}=3.6^{\circ}$ is 1504. We assume therefore that the remaining 1951-1504 = 447 single transiting transiting planets in the Kepler population not fit by this model population are inherently single planet systems.

Despite the model population for $\Delta i_{12}=3.6^{\circ}$ giving the lowest modified $\chi^{2}$ value, this mutual inclination is perhaps larger than that suggested by TDV studies. We note however that mutual inclination estimates from TDV studies consider a range of planet multiplicities. For example Fang & Margot (2012) consider a model population of planetary systems with 1-7+ planets and predict that $\sim$50% of observed planetary systems should contain a single planet, with the remaining systems containing multiple planets with mutual inclinations of $\lesssim 3^{\circ}$. In order to properly predict the inherent mutual inclination in the multi-planet systems considered in this work therefore, it would be necessary to simultaneously model the TDV data directly. We consider such an analysis as part of future work. Instead in §6.4 we consider the possibility that Kepler planets form coplanar, but end up mutually inclined due to perturbations from an outer planetary companion on an inclined orbit. This may provide another way to predict the correct abundance of single transiting systems observed by Kepler, and also result in a low mutual inclination for those systems with two transiting planets.

6.4 Including an inclined planetary companion

We now consider the effects of a hypothetical outer planet in each of the systems in the model population. We first amend the assumption from §6.2 and assume that all stars either host three or zero planets. Any system which hosts three planets is assumed to be identical to one of the 341 double transiting systems from the Kepler population plus an additional outer planet. The outer planet is assumed to have the same mass and semi-major axis in all systems and starts on an inclination to the inner planets when these are coplanar, causing the mutual inclination between the inner planets to evolve according to eq. (15). We assume that the outer planet satisfies the Hill stability criterion of $\Delta$ = 2$\sqrt{3}$ (Chambers (1999)) with the outer of the inner two planets for all 341 considered systems, where $\Delta=(a_{3}-a_{2})/R_{H}$ and

[TABLE]

where $M_{\star}$ is the stellar mass. If this criterion is not satisfied, we move the outer planet for this specific system until it is. For example, when the outer planet is assumed to have a semi-major axis and mass of 1au and 1$\mathrm{M}_{\oplus}$ respectively, we find 6 of the 341 systems do not satisfy this stability criterion and the outer planet needs to be moved to a mean semi-major axis of 1.2au. When the outer planet has a semi-major axis and mass of 1au and 10$\mathrm{M}_{\mathrm{J}}$ respectively we find 22 of the 341 systems do not satisfy the stability criterion and the outer planet needs to be moved to a mean semi-major axis of 1.4au.

Each one of the 341 systems is again replicated enough times in the model population to be expected to be observed exactly once. That is, the inverse of the mean double transit probability of the inner two planets, gives the abundance of each of the 341 systems in the model population. The associated mean single transit probabilities for each of the inner two planets is of the same form as eq. (25). The sum of the mean single transit probabilities for every system in the model population therefore again gives the abundance of a given single transiting planet that would be expected to be observed from the model population that also fits the number of double transiting systems.

Similarly to the modelling approach in §6.3, we look to identify which mass ($m_{3}$), semi-major axis ($a_{3}$) and initial inclination ($\Delta i$) of the outer planet causes the types of single transiting systems expected from the associated model population to be most like those in the observed Kepler population. For a given combination of $a_{3}$, $m_{3}$ and $\Delta i$ we therefore calculate a modified $\chi^{2}$ value described in §6.3. We show these modified $\chi^{2}$ values in Figure 11 for an outer planet with $\Delta i$ = 10∘ (top panel), $m_{3}$ = 1MJ (middle panel) and $a_{3}$ = 2au (bottom panel). Inclinations of $\Delta i\gg 20^{\circ}$ where eq. (15) is expected to break down are included for completeness.

From the top panel in Figure 11, it is clear that there is a ’valley’ of semi-major axes and masses of the outer planet which causes a significantly lower modified $\chi^{2}$ value. It can be assumed therefore that such an additional planet predicts single transiting systems whose radii and semi-major axes better fit those in the Kepler population. However there is also a distinct minimum in the modified $\chi^{2}$ space when the outer planet has a semi-major axis of $\sim$1au for a mass of $\sim$30M⊕. Similarly in the other panels of Figure 11 there appear to be distinct minima. For the middle panel this occurs for an outer planet (of $m_{3}$ = 1MJ) with a semi-major axis of 1.38au, initially inclined to the inner planets by $\Delta i=5.7^{\circ}$. Finally for the bottom panel, this minimum occurs for a mass of $\sim$6MJ and inclination of 6∘ (where $a_{3}=2$au). Generally, we find the distribution of single transiting planets expected from the model population is more representative of those in the Kepler population for $3\lesssim\Delta i\lesssim 10^{\circ}$.

The bottom right panel of Figure 9 gives the distribution of single transiting planets expected from the model population when the outer planet exists in a minimum of the modified $\chi^{2}$ space with $a_{3}$ = 1.07au, $m_{3}$ = 24M⊕ and $\Delta i$=10∘ (white circle in the top panel of Figure 11). We note that the total number of single transiting planets expected from this model population is $1564$. The outer planet parameters which predict $N_{\mathrm{sing}}=1564$ are shown by the white lines in Figure 11. This line highlights that while many outer planet parameters can predict $N_{\mathrm{sing}}=1564$, some predict single transiting planets which are more representative of those in the Kepler population. We note that $N_{\mathrm{sing}}$ predicted by the same range of outer planet parameters from Figure 11 is shown in Appendix C.

7 Discussion

7.1 Combining inherently mutually inclined and outer planet populations

In reality it is likely that the total number of single planet transiting systems observed by Kepler ($N_{\mathrm{sing,Kep}}=1951$) is contributed to by different populations of planetary systems. These may include a number of inherently single planet systems ($N_{\mathrm{sing,inh}}$) in addition to a number of single transiting planets observed from a population of two planet systems which have a fixed mutual inclination of $\Delta i_{12}$ ($N_{\mathrm{sing,\Delta i_{12}}}$). They may also include a number of single transiting planets which are observed from a population of initially coplanar two planet systems interacting with an inclined planetary companion ($N_{\mathrm{sing,planet}}$). Hence in general, it can be considered that

[TABLE]

Here we make the assumption that the total number of double transiting systems observed by Kepler ($N_{\mathrm{doub,Kep}}=341$) is made up of a fraction $f$ that are two planet systems with an inherent mutual inclination and a fraction (1-$f$) that are two planet systems with an inclined outer companion. We can thus rewrite eq. (26) as

[TABLE]

where ($N_{\mathrm{sing,Ndoub=341}}$)${}_{\Delta i_{12}}$ is the number of singles that would have been produced from the population of two planet systems with a fixed mutual inclination of $\Delta i_{12}$, had it been numerous enough to reproduce the 341 double transiting Kepler systems (which is shown in Figure 10 as a function of $\Delta i_{12}$). Conversely ($N_{\mathrm{sing,Ndoub=341}}$)planet is the number of singles that would have been produced from the population of two planet systems which are perturbed by an outer companion, had it been numerous enough to reproduce the 341 double transiting systems. We estimate the number of inherently single planet systems to be $N_{\mathrm{sing,inh}}$ = 447 from §6.3. We note that $N_{\mathrm{sing,inh}}$ will change for different values of $\Delta i_{12}$, however for simplicity we keep it constant at 447.

For the assumed $N_{\mathrm{sing,inh}}$ and an assumed fixed mutual inclination for the fraction of the double transiting systems that are inherently inclined ($f$), eq. (27) means that the number of single transiting systems observed by Kepler can be reproduced by specific combination with the fraction of double transiting systems that have an outer planet ($1-f$) and the properties of these planetary systems which determine the ratio of single to double transiting systems from this population (i.e. ($N_{\mathrm{sing,Ndoub=341}}$)planet). This combination is plotted in Figure 12, which can be read alongside Figure 14 to determine the outer planet parameters required to reproduce the required ($N_{\mathrm{sing,Ndoub=341}}$)planet. For example, for $f=0.2$ and $\Delta i_{12}=2^{\circ}$, ($N_{\mathrm{sing,Ndoub=341}}$)planet = 1676 from Figure 12, which from Figure 14 would be reproduced by an outer planet with $a_{3}=2$au, $m_{3}=132$M⊕ and $\Delta i$ = 10∘. For $f=0.5$, ($N_{\mathrm{sing,Ndoub=341}}$)planet is increased to 2192 requiring the mass of this outer planet to be increased to $m_{3}=955$M⊕ (for $a_{3}=2$au and $\Delta i$ = 10∘). The outer planet parameters required to produce ($N_{\mathrm{sing,Ndoub=341}}$)planet are therefore extremely sensitive to the value of $f$. However, increasing the value of $\Delta i_{12}$ for a given value of $f$ increases the value of ($N_{\mathrm{sing,Ndoub=341}}$)${}_{\Delta i_{12}}$ and hence decreases ($N_{\mathrm{sing,Ndoub=341}}$)planet as seen in Figure 12, requiring an outer planet which is a weaker perturber of the inner planets.

It should be noted that $f$ and $1-f$ are not equivalent to the underlying fraction of stars that host a two planet system with a fixed mutual inclination, or a two planet system with an outer companion respectively. However if $f$ is known, such fractions for the underlying population of stars can be estimated through occurrence rate calculations. We discuss such calculations of occurrence rates in §7.3, however it is first necessary to estimate a value for $f$, which we discuss below.

7.2 Comparing inherently mutually inclined and outer planet populations

From §6.3 a sole population of two planet systems which are inherently mutually inclined by $\Delta i=3.6^{\circ}$ (i.e. when $f=1$) can reproduce a population of single and double transiting systems representative of those observed by Kepler (Figure 9). However from §6.4 a sole population of two planet systems with an outer planet (i.e. $f=0$) can also reproduce a population of single and double transiting systems representative of those observed by Kepler (Figure 11). Here we look to differentiate between these two models by considering the predicted distribution of mutual inclinations that would be observed in the two planet populations for each model. We note that combining these two models in a way described in §7.1 (i.e. when $0<f<1$) would then give some intermediate distribution of mutual inclinations between the overall two planet population.

For the model in which the two planets have an inherent mutual inclination of $\Delta i=3.6^{\circ}$, that distribution is narrowly peaked at 3.6∘ (see Figure 13). In contrast, for the model in which two planets are perturbed by an inclined outer planet, the distribution of mutual inclinations is biased toward coplanar systems. This is because, while the outer planet induces a significant mutual inclination between the inner planets, as required to reproduce the correct ratio of single to double transiting systems, the inclination is not always large (see Figure 5) and the probability of witnessing a double transit system is much higher when their mutual inclination is low. Consider an outer companion with $m_{3}=24$M⊕, $a_{3}=1.07$au and $\Delta i=10^{\circ}$, which was in a minimum of the modified $\chi^{2}$ space (white circle, Figure 11 top). Weighting the secularly evolving mutual inclinations between the inner two planets in the 341 considered systems by the associated double transit probability gives the predicted distribution of mutual inclinations which are most likely to be observed. This distribution is shown by the black line in Figure 13. It is clear that the most likely observed mutual inclination is when the inner two planets are coplanar. Moreover the number of systems expected to be observed with mutual inclinations beyond 0.5∘ drops to negligible values.

From transit duration variation studies, the distribution of mutual inclinations between planets in multi-planet Kepler systems is peaked at $\sim 2^{\circ}$ (Fang & Margot (2012); Fabrycky et al. (2014)), noting however that these works consider different planet populations to those considered here as discussed in §6.3. Combining the two above models to produce a similar distribution in mutual inclinations may therefore allow for $f$ to be determined. We look to combine the two models in such a way to predict a value of $f$, as well as modelling the TDVs of the planetary systems considered in this work directly to predict the distribution of inherent mutual inclinations, as part of future work. For example if a fraction of two planet systems observed by Kepler are considered to have a fixed mutual inclination of $\Delta i_{12}=4^{\circ}$, then in order to reproduce a distribution of mutual inclinations that peaks at $\sim$2∘ from modelling of TDVs, it might be expected that $f\sim 0.5$.

An additional method to estimate $f$ might be to consider whether hypothetical outer planets considered in this work would have been detectable by other means. It is expected that RV studies would be most sensitive to such outer planetary companions. On Figure 11 we plot a rough constraint from RV studies, shown by the red dashed lines, assuming a detection threshold of $\sim$2m/s. Outside of 5au we assume RV studies are not sensitive to planets due to long periods. Planets above or to the left of these lines would therefore be detectable with this level of RV precision. This detection threshold suggests that a wide orbit planet located in the minima of the modified $\chi^{2}$ values in Figure 11 (white circle) should be just detectable by RV studies. This would assume however that all Kepler systems with two planets host this outer companion i.e. $f=0$. From Figure 12 and highlighted in §7.1, if $f>0$ a planet with a larger mass, shorter period or larger inclination is required to reproduce the total number of single transiting systems observed by Kepler. Such outer planets should be readily detectable by RV surveys. For example, for the values of $f=0.2$ and $f=0.5$ for $\Delta i_{12}=2^{\circ}$ considered in §7.1, both of the outer planets in these cases would be expected to be detectable by RV surveys. Due to the inherent faintness of Kepler stars, few have been extensively studied for wide orbit planets. We suggest therefore that detailed follow-up RV studies of Kepler systems would allow for $f$ to be constrained. Generally for example, a low yield of outer planets in RV studies would suggest that $f$ is low and vice versa.

7.3 Occurrence Rates

Similar to that discussed specifically for Kepler systems in §7.1, consider that the underlying population of planetary systems contains three possible types of planetary systems. These include inherently single planet systems, two planet systems which have a fixed mutual inclination of $\Delta i_{12}$ and two planet systems which are being perturbed by an inclined outer planet. In §7.1 it was shown that combining these systems with a free parameter $f$, which describes the fraction of the observed double transiting population that are two planet systems with a fixed mutual inclination, recovers the total number of single and double transiting systems observed by Kepler.

However this value of $f$ is not the same as the fraction of the underlying population of stars that have two planets that are inherently mutually inclined. Here we define the occurrence rate of a given population to be the fraction of stars which would be expected to host such systems. Occurrence rates in this work can be estimated by taking the ratio of the number of systems in a given model population ($N_{\mathrm{mod}}$) to the total number of stars observed by Kepler ($N_{\mathrm{Kep}}$). The individual occurrence rates for the inherently single planet systems is therefore given by ($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)inh, for the two planet systems with the fixed mutual inclination of $\Delta i_{12}$ by ($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)${}_{\Delta i_{12}}$ and for the two planet systems being perturbed by an inclined outer planet by ($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)planet. For example, for the population of two planet systems which were inherently mutually inclined by 3.6∘ (for when $f=1$), i.e. those which predicted a population of single transiting planets representative of those observed by Kepler (§6.3), the number of systems in the model population was equal to 43807. From §6.1 the total number of Kepler stars was 164966. Therefore the occurrence rate for this type of system, ($N_{\mathrm{mod}}/N{{}_{\mathrm{Kep}}}$)${}_{\Delta i_{12}}$ = 27%. Conversely, considering the population of two planet systems which were perturbed by an outer companion with $m_{3}=24$M⊕, $a_{3}=1.07$au and $\Delta i=10^{\circ}$ (white circle Figure 11 top) for when $f=0$, predicted 42733 systems in the associated model population. Therefore the associated occurrence rate of this type of system ($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)planet = 26%.

The calculation of the occurrence rate for the population of inherently single planet systems is slightly different to that described above. From §6.3, assume that there are 447 inherently single planet systems (noting that this is subject to the value of $\Delta i_{12}$). The distribution of the semi-major axes of these 447 planets is equal to the difference between the distributions of semi-major axes for the single transiting systems observed by Kepler and those predicted by the population of two planet systems with a fixed mutual inclination of $\Delta i_{12}=3.6^{\circ}$, i.e. the difference between the left panel of Figure 8 and the bottom left panel of Figure 9. The number of inherently single planet systems in a model population is then the sum of the inverse of the single transit probabilities ($R_{\star}/a$) of all these 447 planets. We find this model population contains 15852 systems, predicting an occurrence rate of inherently single planet systems of 9.6%. This is large compared with the occurrence rate of Hot Jupiters ($\sim 1-2\%$ e.g. Marcy et al. (2005); Cumming et al. (2008); Mayor et al. (2011); Wright et al. (2012); Santerne et al. (2016)). We therefore expect that our population of inherently single planet systems is dominated by a different population, such as those described in §6.1 which are poorly constrained.

In a similar way to that described for eq. (27), the total occurrence rate of assumed planetary systems in the underlying population of planetary systems can be estimated to be

[TABLE]

Consider the example combination of systems from §7.2 for when $f=0.2$, $\Delta i_{12}=2^{\circ}$ and the outer planet parameters are $a_{3}=2$au, $m_{3}=132$M⊕ and $\Delta i=10^{\circ}$. Here $f$($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)${}_{\Delta i_{12}}$ $\sim$ 3% and ($1-f$)($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)${}_{\mathrm{planet}}\sim$21%. We note that $f$($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)${}_{\Delta i_{12}}$/($1-f$)($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)planet = 3/21 = 14%. This highlights that the occurrence rate of stars which have two planet systems with an inherent mutual inclination is similar to, but not the same as the parameter $f$.

Combining with the occurrence rate of inherently single planet systems estimated above, the total occurrence rate of planetary systems becomes 34%. This is similar to occurrence rates of $\sim 25\%-30\%$ for Kepler like planets derived from injection and recovery analysis of planet candidates from the Kepler pipeline (Petigura et al. (2013); Christiansen et al. (2015)).

Estimates of occurrence rates for planets similar to the outer planets considered in this work exist from RV studies. Cumming et al. (2008) suggest an occurrence rate of 7.0 $\pm$ 1.4% for planets with masses and semi-major axes of $m_{\mathrm{p}}$ = 1-10MJ and $\sim$1-5au respectively. Extrapolating this occurrence rate also predicts that 17-20% of stars have gas giants within 20au. Similarly Mayor et al. (2011) suggest an occurrence rate of 13.9 $\pm$ 1.7% for planets with masses and periods of $m_{\mathrm{p}}>50$M⊕ and $P<10$yrs respectively. More recently Bryan et al. (2016) suggest that for systems with 1 or 2 RV planets, the occurrence rate of an additional companion with a mass and semi-major axis of 1-20MJ and 5-20au respectively is as high as 52 $\pm$ 5%. The above example occurrence rate for the systems with an outer planet, i.e. ($1-f$)($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)${}_{\mathrm{planet}}\sim$21%, is then therefore not contradicted by these studies. However, this example assumed an estimated value of $f$. In addition to the methods described in §7.2, observationally estimated occurrence rates for outer planets may also be able to constrain the value of $f$. For example if it is assumed that the occurrence rate of the types of outer planets considered in this work is 13.9% (Mayor et al. (2011)), then it can be estimated that ($1-f$)($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)${}_{\mathrm{planet}}\sim 13.9\%$. As ($N_{\mathrm{mod}}/N_{\mathrm{Kep}}$)${}_{\mathrm{planet}}\ngtr$ 1 (i.e. it is unphysical that there are more stars in the model population than the number actually observed by Kepler), this results in an upper limit of $f\leq 0.86$. We suggest therefore that combining this method of placing constraints on $f$ with those described in §7.2 might provide a strong constraint on the percentage of planetary systems which may share a fixed mutual inclination compared with systems that may host an outer inclined planet.

7.4 Comparing with similar works

Whether an outer planet can reduce the multiplicity of expected transiting planets in an inner planetary system in the context of N-body simulations has recently been investigated by Hansen (2017). A notable example they include is the effect of a companion with a mass of 1MJ at 1au, which is inclined to an inner population of planetary systems with a variety of multiplicities by 10∘. They find the ratio of the total number of double to single transiting systems that Kepler would be expected to observe is 0.184 (i.e. $\sim$5 times more expected single than double transiting systems). We find an identical outer planetary companion in our work gives this ratio to be 0.14. We suggest this difference is caused by the population of inner planetary systems used. Hansen (2017) incorporate 50 model inner planetary systems with a range of multiplicities (the vast majority contained 3-6 planets at the end of their simulations), rather than the two planet Kepler systems considered in this work. Higher multiplicities increases the number of competing secular modes in the system, which can stabilise inner planets against the secular perturbations of an outer companion (e.g. Read & Wyatt (2016)). Such an example was shown in this work in §5 for application to Kepler-48. Perhaps then, mutual inclinations are more easily induced between inner planets in this work, increasing the predicted number of single transiting planets that Kepler would be expected to observe, relative to a fixed population of planetary systems.

Moreover compared with N-body simulations, our work does not allow for dynamical instability. If inclinations are large then they couple with eccentricity (Murray & Dermott (1999)), potentially causing orbital crossings between neighbouring planets leading to dynamical instabilities on short, non-secular timescales. Indeed Hansen (2017) find for the above mentioned outer planetary companion that roughly half of the 50 systems they consider lose at least one planet. Moreover Pu & Wu (2015) suggest that the abundance of single and double transiting systems might be the remains of higher order planetary systems that were once tightly packed and have since undergone dynamical instability. A detailed discussion on how dynamical stability would be expected to affect our results is difficult. Our choice that all planets must be initially Hill stable is by no means a robust constraint on the long term stability of all the planetary systems we consider during the secular interaction.

The effects of dynamical instability in tightly packed planet systems interacting with a wide orbit companion planet was also shown by Mustill et al. (2015). They find that an outer giant planet undergoing Kozai-Lidov interactions with a stellar binary (Kozai (1962); Lidov (1962)) can have an eccentricity which takes its orbit within the inner planets, leading to a significant reduction in planet multiplicity. Moreover more recent work in Mustill et al. (2016) suggests that these same interactions can cause $\sim$50% of Kepler like systems to lose a planet, either through collisions or ejections. If inclination is not completely decoupled with eccentricity then, these works suggest that dynamical instability plays a significant role in sculpting an inner planetary system.

7.5 Metallicity Distribution

The fraction of stars with gas giants increases with higher metal content (e.g. Gonzalez (1996); Thorngren et al. (2016)). However it is unclear if this relation extends to smaller planets with $R_{p}\lesssim 4R_{\oplus}$ (Mayor et al. (2011); Zhu et al. (2016)). If single transiting planets are in systems which contain an outer giant companion similar to that considered in this work then the transiting planet should follow a similar metallicity relation as the giant planet. If there is an inherent population of single planet systems with $R_{p}\lesssim 4R_{\oplus}$, in addition to a population of inherently mutually inclined double transiting systems, then these systems will follow a different metallicity relation. Therefore the population of single and double transiting systems observed by Kepler may contain a mixture of metallicity relations. If a distinction can be made between these different relations then this may place constraints on the presence of additional planets in Kepler systems with a single transiting planet.

7.6 Assumptions of this work

Throughout this work we have considered mutual inclinations evolve between two planets due to secular interactions with an outer planet. As stated above, increasing the multiplicity of planetary systems complicates the evolution of mutual inclinations. For application to the Kepler dichotomy, including higher multiplicity systems may cause proportionally fewer to be observed as single transiting systems. We look to investigate this as part of future work. Moreover higher multiplicity systems also allow for investigation into whether the presence of an outer planetary companion can explain the number of higher order systems observed by Kepler. This is of particular interest as Johansen et al. (2012) find that generating a model population which predicts the number of systems observed by Kepler with three transiting planets (with small inherent mutual inclinations and no outer companion) cannot simultaneously predict the number of systems with a single and two transiting planets observed by Kepler.

We have also assumed that the inner transiting planets interacting with an outer companion were initially coplanar. However these transiting planets would most likely also have a small inherent mutual inclination (e.g. Fang & Margot (2012); Fabrycky et al. (2014)) which in turn may affect the mean double transit probability.

8 Summary and Conclusions

In summary, during the first part of this work we developed a semi-analytical method for the calculation of transit probabilities by considering the area a transiting planet subtends on a celestial sphere (§2). Applying this method to a general two planet system, we showed how the probability that both planets are observed to transit changes as they become mutually inclined.

In §3 we discussed how the mutual inclination between two initially coplanar planets evolves due to secular interactions with an external mutually inclined planetary companion. We derived the full solution describing this evolution assuming that the mutual inclination remains small, before simplifying it under the assumption that the external planet was on a wide orbit. We found that the maximum mutual inclination between the inner two planets is approximately equal to twice the initial mutual inclination with the external planet. Below this the maximum mutual inclination between the inner two planets scales according to the mass, semi-major axis and inclination of the external planet by $\propto\Delta im_{3}/a_{3}^{3}$.

How the secular interaction causes the double transit probability of the inner two planets to evolve was shown in §4. Assuming that this double transit probability is significantly reduced when the maximum mutual inclination exceeds $\approx(R_{\star}/a_{1})+(R_{\star}/a_{2})$ we derived an expression for the mean of the double transit probability considering a given external planetary companion. This expression was applied to Kepler-56, Kepler-68, and Kepler-48 to place constraints on the inclination of the outer RV detected planets in these systems in §5. We found that the inner two transiting planets in Kepler-56 and Kepler-68 are not significantly secularly perturbed by the outer planets, regardless of their inclination. For HD 106315 we find that an outer planet inferred from recent RV analysis can cause a significant perturbation to the mutual inclination of two internal transiting planets. Moreover we find that if the outer planet is present within $\sim$1au, its inclination must be no more than $2.4^{\circ}$, otherwise the probability of observing both the inner planets to transit is significantly reduced. We also found that the RV detected planet in Kepler-48 needs to be inclined with respect to the inner planets by $\lesssim 3.7^{\circ}$, otherwise the probability that all the inner planets are observed to transit is significantly reduced. We conclude therefore that using the expression for the mean transit probability between inner planets from eq. (20) and (17) can be used to place significant constraints on the inclinations of RV detected planets, whose host systems also contain transiting planets.

We further applied our method of calculating transit probabilities to the Kepler population in §6. We found that relative to a fixed population of transiting systems with two planets on initially coplanar orbits, the expected number of single transiting systems can be significantly increased both by inherently inclining the two planets and by introducing an outer planetary companion. We found that an inherent mutual inclination of $\Delta i_{12}=3.6^{\circ}$ predicts a population of single transiting planets most representative of those in the Kepler population. Moreover, we found that outer planets initially inclined by $\sim 3-10^{\circ}$ to the inner planets also predict a representative population of single transiting systems. These outer planets should be detectable by RV studies.

However it is likely that planetary systems observed by Kepler may include a combination of systems which include inherently single planet systems, two planet systems which have some fixed mutual inclination and two planet systems interacting with an inclined outer planet. For two planet systems which are perturbed by an outer planet, the distribution of the mutual inclinations between the inner planets of such systems is biased toward coplanar systems. This is due to an increased probability of observing inner planets when coplanar compared with when mutual inclinations are larger. We suggest that combining populations of inherently mutually inclined two planet systems with two planet systems which are interacting with an outer planet may be able to reproduce the observed distribution of mutual inclinations between Kepler planets. In doing so, this may provide constraints on the presence of outer planets in the Kepler population. We suggest also that detailed follow-up of RV studies in Kepler systems will provide a more direct constraints on the presence of outer planets. There should also be a dichotomy in the number of transiting systems observed by the upcoming TESS mission (Ricker et al. (2014)), however for these systems astrometry and RV techniques will be able to be used to verify the presence and influence of outer planets.

Acknowledgements

We thank Simon Gibbons and Grant Kennedy for useful conversations regarding this work. MJR acknowledges support of an STFC studentship and MCW acknowledges the support from the European Union through grant number 279973. We also thank the reviewer for comments which were a great help in improving this paper. This research has also made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Appendix A Further discussion of Transit equations

A.1 Central Transit Line

The centre of the transit region is defined by eq. (2)

[TABLE]

Assuming that $\phi=0\rightarrow 2\pi$ and that a corresponding value of $\theta$ for each value of $\phi$ can be in the range of $0<\theta<\pi$, eq. (2) can be rearranged to give

[TABLE]

A.2 Upper Transit Boundary

The upper boundary of a transit region is given by eq. (4)

[TABLE]

A value of $\theta_{2}$ for a given $\phi_{2}$ can be calculated through solving a quadratic of the form

[TABLE]

where

[TABLE]

Depending on the value of $\Delta i$, the calculation of $\theta_{2}$ for a given value of $\phi_{2}$ can be grouped into three different regimes, (1) $\Delta i$ is small enough that the upper boundary of the transit region never crosses the fixed reference plane. (2) $\Delta i$ is large enough that the upper boundary of the transit region does cross the fixed reference plane. (3) For high values of $\Delta i$ the upper transit boundary only has values of $\theta_{2}$ for $0<\phi_{2}<\pi$. This can be thought of as the transit region going over the pole of the celestial sphere.

For regime (1), the value of $\theta_{2}$ for a given $\phi_{2}$ is equivalent to that obtained from the positive root of eq. (30), mirrored about $\pi/2$. The transition to regime (2) occurs for when the upper transit boundary first crosses the fixed reference plane. Here $\Delta i=\arcsin(\chi)$. As $\Delta i$ is increased beyond this value the intersection between the upper transit boundary and the fixed reference plane occurs at $\phi_{2}=\phi_{0}$ and $\phi_{2}=\pi-\phi_{0}$, for which $\theta_{2}=\pi/2$. From eq. (4) $\phi_{0}$ is given by $\phi_{0}=\arcsin(\chi/\sin\Delta i)$. Therefore $\theta_{2}<\pi/2$ for $\phi_{0}<\phi_{2}<\pi-\phi_{0}$ and $\theta_{2}>\pi/2$ otherwise. When $\phi_{0}<\phi_{2}<\pi-\phi_{0}$, $\theta_{2}$ is hence obtained from the positive solution of eq. (30) and by the positive solution mirrored about $\pi/2$ otherwise.

Finally the transition to regime (3) occurs when $\Delta i=\arccos(\chi)$. Similarly to regime (2) as $\Delta i$ is increased beyond this value, the upper transit boundary crosses the fixed reference plane at $\phi_{2}=\phi_{0}$ and $\phi_{2}=\pi-\phi_{0}$ and hence $\theta_{2}$ is only defined for when $\phi_{0}<\phi_{2}<\pi-\phi_{0}$. The solution from eq. (30) which gives the smaller value $\theta_{2}$ corresponds to $\theta_{2}>\pi/2$ values and needs to be mirrored about $\pi/2$, with the solution giving the larger value of $\theta_{2}$ corresponding to $\theta_{2}<\pi/2$ values.

To summarize consider that for a given value of $\phi_{2}$, eq. (30) gives two solutions for $\theta_{2}$, denoted as $\theta_{2}^{*1}$ and $\theta_{2}^{*2}$ respectively. For $\Delta i<\arcsin(\chi)$,

[TABLE]

where $\theta_{2}^{*1}>0$ and $\theta_{2}^{*2}<0$.

For $\arcsin(\chi)<\Delta i<\arccos(\chi)$,

[TABLE]

where $\theta_{2}^{*1}>0$, $\theta_{2}^{*2}<0$ and $\phi_{0}=\arcsin(\chi/\sin\Delta i)$.

For $\Delta i>\arccos(\chi)$,

[TABLE]

where $\theta_{2}^{*1}>0$, $\theta_{2}^{*2}>0$.

A.3 Lower Transit Boundary

The lower boundary of the transit region is given by eq. (3)

[TABLE]

Depending on the value of $\Delta i$, the calculation $\theta_{1}$ for a given $\phi_{1}$ can be grouped into the same regimes as described for the upper transit boundary. However now in regime (1), $\theta_{1}<\pi/2$ for $0<\phi_{1}<2\pi$, in regime (2) the lower transit boundary crosses the fixed reference plane at $\phi_{1}=\pi+\phi_{0}$ and $\phi_{1}=2\pi-\phi_{0}$ and in regime (3) $\theta_{1}$ is only defined for $\pi+\phi_{0}<\phi_{1}<2\pi-\phi_{0}$. Assuming that $\theta_{1}^{*1}$ and $\theta_{1}^{*2}$ are the solutions for $\theta_{1}$ for a given $\phi_{1}$ in the modified form of eq. (30), then following the same discussion as for the upper transit boundary it can be shown that for $\Delta i<\arcsin(\chi)$,

[TABLE]

where $\theta_{1}^{*1}>0$ and $\theta_{1}^{*2}<0$.

For $\arcsin(\chi)<\Delta i<\arccos(\chi)$,

[TABLE]

where $\theta_{1}^{*1}>0$, $\theta_{1}^{*2}<0$ and $\phi_{0}=\arcsin(\chi/\sin\Delta i)$.

For $\Delta i>\arccos(\chi)$,

[TABLE]

where $\theta_{1}^{*1}>0$, $\theta_{1}^{*2}>0$.

Appendix B Secular Solution for Mutual Inclination Evolution

From eq. (11) the evolution of complex inclinations according to Laplace-Lagrange theory is given by

[TABLE]

where Ijk are the eigenvectors of the matrix B from eq. (9) scaled to initial boundary conditions, $f_{i}$ are the eigenfrequencies of B and $\gamma_{k}$ are initial phase terms. If it is assumed that all the planets and the star are point masses and that the invariable plane is taken as a reference plane, it follows that $f_{3}=0$ and I${}_{j3}=0$. From the initial conditions |$y_{1}(0)$| = |$y_{2}(0)$| = $i_{1}$. Hence the complex inclinations of the inner two planets respectively are given by

[TABLE]

[TABLE]

Also from the initial conditions $-I_{11}+I_{12}=i_{1}$ and $-I_{21}+I_{22}=i_{1}$. The complex mutual inclination between the inner two planets is equivalent to

[TABLE]

Solving eq. (11), we propose a set of variables to represent the full solution of $I_{12}$ and $I_{22}$,

[TABLE]

where $m$ = 1, 2,

[TABLE]

[TABLE]

Hence the components of the eigenvector associated with the $f_{2}$ eigenfrequency are given by

[TABLE]

The non zero $f_{1}$ and $f_{2}$ eigenfrequencies of the matrix B from eq. (9) can be obtained by solving a quadratic of the form

[TABLE]

We note that the solution given by eq. (44) recovers exactly what is predicted when solving eq. (11) by numerical methods. The full solution which describes how the mutual inclination between the inner two planets according to Laplace - Lagrange theory is therefore given by

[TABLE]

with the variable $K$ used in §3.2 being equivalent to $(1-R_{12})(1-R_{11})/\epsilon$.

Appendix C Reproducing the total number of single transiting planets observed by Kepler

In §6 we considered Kepler systems with two transiting planets which are secularly interacting with an outer planet on an inclined orbit. We found that the number of single transiting systems Kepler would be expected to observe can be dramatically increased as a result of this interaction. Figure 14 shows the total number of single transiting objects Kepler would be expected to observe from the method outlined in §6.4, for when the outer planet has the same parameters as the respective panels of Figure 11. Again for $\Delta i\gg 20^{\circ}$, Laplace-Lagrange theory is expected to break down and is included for completeness. The white line gives where the total number of single transiting planets Kepler would be expected to observe from the model population is equal to the number in the Kepler population i.e. 1951. The red dashed lines give an estimate for an RV detection threshold.