Transit Probabilities in Evolving Secular Systems

Transit Probabilities in Secularly Evolving Planetary Systems

Matthew J. Read Mark C. Wyatt and Amaury H. M. J. Triaud
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA
E-mail: mjr201@ast.cam.ac.uk
Accepted XXX. Received YYY; in original form ZZZ
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.

keywords:
planets and satellites: dynamical evolution and stability
pubyear: 2017pagerange: Transit Probabilities in Secularly Evolving Planetary SystemsTransit 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 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 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 R 4R is low at 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 (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 and radius 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 () on a unit sphere is equal to

(1)

where is the polar angle and is the azimuthal angle. A given area on the celestial sphere can therefore be represented on a 2d plane of vs. , from 0 2 and 0 respectively, such that the 2d plane has a total surface area of . 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

Figure 1: The coordinate system used to show how a transit region traverses the surface of a celestial sphere. The dashed line represents an orbital plane inclined to a fixed reference plane by . The direction is normal to the orbital plane. The directions , and trace the central, lower and upper boundaries of a transit region respectively.
Figure 2: The surface of a celestial sphere represented on a 2d plane. The dotted lines represent the centre of a transit region for a planet inclined to a fixed reference plane by . The solid lines refer to the boundaries of such transit regions for when . The area within these transit regions are identical, giving an identical single transit probability equal to 0.25.

Consider some fixed reference plane where [] define a pair of orthogonal directions in this plane, and 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 . 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 and radius , is inclined to the fixed reference plane by , with the intersection between the two planes occurring along the direction. The direction of the normal of the orbital plane is given by . The position of a planet in the orbital plane is defined by the direction which makes the angles and with the and directions respectively. Hence traces the centre of the transit region with respect to the fixed reference plane. As , where and it follows that

(2)

Hence eq. (2) defines how the centre of a transit region inclined to a fixed reference plane by traverses a celestial sphere. This is shown by the dashed lines in Figure 2 for different values of , 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 = 90, can only take values of 0 or .

Similarly the directions that define the boundaries of the transit region can be given by and which makes the angles and with the and directions respectively, shown in Figure 1. The boundaries of the transit region also subtend an angle from the orbital plane where assuming (Borucki & Summers (1984)). As , and and , it follows that

(3)
(4)

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 traverse a celestial sphere. The solid lines in Figure 2 show these boundaries for different values of , where = 0.25. This value of might be considered to be unrealistically large and is used for demonstration purposes only. In Appendix A we further discuss how the values of () and () in eq. (3) and eq. (4) respectively would be expected to change as is increased from .

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 (, Borucki & Summers (1984)). All of the transit regions shown in Figure 2 for different 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 for grazing and full transits respectively. Throughout this work however we assume that .

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 , and respectively, where and the orbital planes are mutually inclined by (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 gives the double transit probability. How the double transit probability changes as a function of is shown by the blue line in Figure 3, for when and . 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 , the double transit probability ( herein) can be split into three regimes (also discussed in Ragozzine & Holman (2010); Brakensiek & Ragozzine (2016)).

(1) For low values of , the transit region of the outer planet is enclosed within that of the inner planet. The double transit probability is therefore equal to .

(2) 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 , which causes in eq. (3) for both planets to be equal at . Evaluating eq. (3) at this point gives

(5)

where and 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 , the transit regions only overlap at the intersection of the two orbital planes. The transition to this regime occurs when , where for the inner planet is equal to for the outer planet at . Evaluating eq. (3) and (4) here gives

(6)

The values of and 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/ factor here that is not included in Ragozzine & Holman (2010).

(7)

For large therefore, the double transit probability predicted by eq. (7) tends to a value of 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 1.

For 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 , 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 . 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.

Figure 3: The double transit probability as a function of mutual inclination between two planets from our method (blue line) for when and . The dashed black lines represent the associated analytical estimate given by eq. (7). The green and red lines represent which inclination cause the double transit probability to go from regime 1 to 2 and regime 2 to 3, with the regimes being defined in §2.2.

3 Secular Interactions

3.1 N planet system

Consider a system of secularly interacting planets in which planet has a semi-major axis and mass . The inclination and longitude of ascending node of planet are given by and respectively, and can be combined into the associated complex inclination . Assuming that the vector involving all the planet’s orbital planes is given by = [], the evolution of complex inclinations in the low inclination and eccentricity limit can be given by Laplace - Lagrange theory in the form

(8)

(Murray & Dermott (1999)) where is a matrix with elements given by

(9)

and are integers associated with each planet, and are the masses of the star and planet , is the mean motion of planet where , = for and and otherwise, and corresponds to a Laplace coefficient given by

(10)

Eq. (8) can be solved to show that the evolution of is given by a superposition of eigenmodes associated with each eigenfrequency of the matrix

(11)

where I are the eigenvectors of scaled to initial boundary conditions and is an initial phase term. If it is assumed that all objects are spherically symmetric, additional terms in the diagonal elements of 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

Figure 4: The maximum mutual inclination, , between two planets on circular, initially coplanar orbits with semi-major axis of 0.2, 0.5au and masses of 10M respectively, from the secular interaction with an outer third planet. The value of calculated by the full Laplace-Lagrange solution from eq. (15) is given by the colour scale on the left panels. The right panel colour scales give calculated by the simplified Laplace-Lagrange solution for when , given by eq. (16) and eq. (17). For the top panels , for the middle panels M and for the bottom panels au. It is important to note that the assumptions of Laplace-Lagrange theory break down when . Larger inclinations are only included in this Figure to aid comparison between predicted by the full and simplified Laplace-Lagrange theory solutions.

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 and respectively such that . The orbital plane of this external planet is initially mutually inclined to the inner planets by . 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 is given by

where 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 .

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

(12)

where and 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

(13)

The boundary conditions give and . Also as ( = 0) = ( = 0) = , it follows from eq. (12) that = . The evolution of the mutual inclination from eq. (13) is therefore is equivalent to

(14)

Hence the evolution of the instantaneous mutual inclination between the inner pair of planets, , can be calculated if the first and second elements of the eigenvector associated with the eigenfrequency are known. In Appendix B we fully solve eq. (8) to give and in terms of physical variables. Here we simply say that

(15)

where 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 , 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|| = 2. 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, , and .

The left panels of Figure 4 show how max|| changes as a function of different combinations of , and in eq. (15) for an example system where , = 0.2, 0.5au and , = 10M respectively. We note that the assumptions of Laplace-Lagrange theory are expected to break down when . 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 , certain individual and combinations of matrix elements from eq. (9) become small and we find that eq. (15) can be simplified to

(16)

where

(17)

Here it is assumed that as , certain Laplace coefficients from the matrix elements can be simplified, specifically (Murray & Dermott (1999)). Similar simplifications can be made to each of the eigenfrequencies, for which

(18)

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) we assume that the maximum value of the mutual inclination between the inner two planets predicted by eq. (16) is

(19)
Figure 5: (left): The evolution of the mutual inclination of the two inner planets considered in Figure 4 due to secular interactions with a third planet with = 2au, = 1M and = 5. (right): The associated evolution of the double transit probability.

The right panels of Figure 4 show max|| predicted by eq. (19) and eq. (17) using the same planet parameters as shown in the left panels. We find that when 1.25au, the simplified form for max|| from eq. (19) and eq. (17) agrees with the full Laplace - Lagrange solution to within for all values of and . For 1au, the simplified form of max|| begins to break down and eq. (19) can underestimate max|| 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 in their eq. 12). Assuming inclinations are small we find eq. (19) agrees with the equivalent prediction of max|| from Lai & Pu (2017) if . Indeed, and are almost identical despite the different derivation techniques (e.g. we derive the full Laplace- Lagrange solution and then simplified assuming ), apart from contains an additional factor of whereas contains a factor of (). By considering different combinations of and and comparing to the value of max|| 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, , 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 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 = 2au and = 1M and = 5 respectively. Indeed, is only at a maximum value for a small proportion of the secular evolution leading to a significant reduction in compared with if the outer planet were not present.

Figure 6: The mean double transit probability of two planets from Figure 4, which are being secularly perturbed by a third planet on a mutually inclined orbit according to the full Laplace - Lagrange solution (left panels) and the simplified Laplace-Lagrange solution for when the third planet is assumed to be on a wide orbit. The black lines show the boundary where the maximum mutual inclination between the inner planets exceeds from eq. (5) and is assumed to be significantly reduced. The black lines on the respective left and right panels are identical and included to aid comparison. As noted in Figure 4, Laplace - Lagrange theory is expected to break down for . Larger inclinations are only included here for demonstration purposes only.

Furthermore, the left panels of Figure 6 show how 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|| also causes a significant reduction in the mean double transit probability of the inner two planets and vice versa for small values of max||.

The left panels of Figure 6 also suggest a clear boundary of , and , above which the outer planet causes to be significantly reduced and below which is unchanged. From Figure 3, the double transit probability of the two inner planets can be considered to be significantly reduced when , where is given by eq. (5). We assume therefore that the boundary where is significantly reduced occurs when max|| . The values of , and 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 , here we estimate the value of the mean double transit probability itself. We assume that is dominated by the maximum or minimum value of the double transit probability, and respectively, depending on whether max|| is greater than . We assume that from eq. (5) for , 1. From Figure 3, the value of = , however a value of is more difficult as no specific analytical estimate exists. We therefore assume can be given by the estimate from Ragozzine & Holman (2010) shown by eq. (7). We note that this approximation for would be expected to break down if max|| 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|| is given by the simplified Laplace - Lagrange solution from eq. (19), can be estimated by

(20)

The right panels of Figure 6 show the value of 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 . As such, eq. (20) can under predict , 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 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

(21)

where 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, 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 0.13 M and R = 4.23 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.1M and 181M respectively (Huber et al. (2013)). Keck/HIRES and HARPS-North observations have revealed a non-transiting giant planet (d) with a period of 10025 days and minimum mass of 5.620.38M (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 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 . 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.080.05M and 1.240.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.971.70 and 2.183.5M respectively (Marcy et al. (2014)). Keck/HIRES RV follow up of this system detected a non-transiting planet (d) with a period of 62516 days with a fitted mass of 26716M (Marcy et al. (2014)).

Applying eq. (21) we find . 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 (Kane (2015)).

5.3 Hd 106315

HD 106315 is a bright F dwarf star at a distance pc (Gaia Collaboration et al. (2016)) with mass and radius of 1.070.03M and 1.180.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 and 3.95 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 days, which has a mass of 1M (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 days and days respectively. Assuming days implies a mass of M (Winn et al. (2009); Crossfield et al. (2017)). Applying eq. (21) with this outer planet gives . This suggests that if the outer planet had a period of days, it must have an inclination of , 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 days, implying a mass of 7M, eq. (21) suggests that . That is, if the outer planet has a period of days, it must have an inclination of , 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

(22)

with the astrometric signal-to-noise equal to , where is the scheduled number of astrometric measurements ( = 36 for HD 106315222http://gaia.esac.esa.int/gost/) with typical uncertainties of (de Bruijne (2012)). If , the orbital inclination can be constrained to a precision of (Sahlmann et al. (2015)). We find that for the example periods and masses considered above for HD 106315d that . 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 1.3au, (implying a mass of ) eq. (22) suggests that 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.880.06M and R = 0.890.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.942.10, 14.612.30 and 7.934.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 9828 days and 657 25M respectively (Marcy et al. (2014)).

Figure 7: The mutual inclination between the respective planets in Kepler-48, when the non-transiting planet, Kepler-48e is initially mutually inclined by . The black dashed line shows the evolution of the mutual inclination between the inner two transiting planets with the outer transiting planet, for when the inner two planets are treated as a single body with an equal orbital angular momentum.

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

(23)

where and is proportional to the angular momentum of Kepler-48e in the low eccentricity limit and for either Kepler-48b, c, or d. The initial inclination of the transiting planets is therefore equal to .

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 . The red line shows the mutual inclination between Kepler-48b and c (), the blue between b and d (), and the green between c and d (). 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 () 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, , can therefore be approximated by eq. (21), where the value of becomes

(24)

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 = 3.7. This suggests therefore that the inclination of Kepler-48e , 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|| , Lai & Pu (2017) also find that the inclination of Kepler-48e, considering secular interactions only, must also be small with .

Figure 8: The smoothed distribution of the radii and the semi-major axes of planets observed by Kepler to be in systems with a single transiting planet (left) and in systems with two transiting planets (right). Pixel sizes are log() = 0.15 by log() = 0.1.

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 () 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 7000K and 4.0 log() 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 (Borucki et al. (2011)), leaving 3757 objects, after which we remove candidates with a SNR 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 = (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.

Figure 9: The distribution of the radii and semi-major axis of single transiting planets observed from the model population with: (top) no third planet. (middle) A third planet with = 1M, = 1.9au and = 10. The total number of single transiting planets predicted by the model population is equal to that observed by Kepler. The colour scale for this panel is saturated for ease of comparison. (bottom) A third planet with = 24M, = 1.07au and = 10. We find the 1564 single transiting planets predicted here are a best fit to those observed by Kepler (left panel of Figure 8). The contours show the distribution of single transiting planets from the Kepler population. Pixel sizes are log() = 0.15 by log() = 0.1.

Finally there appears to be an over abundance in the population of single transiting systems for planets with at periods days (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 we use the rocky planet mass-radius relation from Weiss & Marcy (2014), where density () is related to radii () through gcm. For radii 1.5 4, we use the deterministic version of the probabilistic mass-radius relation for sub-Neptune objects from Wolfgang et al. (2016), where mass () is given by . Once radii become 4 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 gives the most sensible masses for all planets with .

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, , that would be expected to be observed. Here the mean single transit probability for a given system is equal to , where , are the semi-major axes of each planet when and is the radius of the host star. We find = 589, which clearly underestimates the 1951 single transiting systems in the observed Kepler population, by a factor of . 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.

Figure 10: (top) The expected number of single transiting planets observed from a model population generated from Kepler systems with two planets that are mutually inclined by . The number of double transiting systems predicted by the model population is constant with 341 systems. (bottom) The associated modified comparing types of single transiting planets predicted by the model population with the Kepler population. The minimum modified value corresponds to .

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 (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 between the two planets in each of the 341 double transiting systems. The mean single transit probability for each planet from a given system, and respectively where , is now given by

(25)

where is the mean double transit probability and is the total mean single transit probability for this system. As 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 increases with for when the number of double transiting systems is kept constant at 341 systems. If , we find , 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 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 R and semi-major axes of 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 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 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 we therefore look to identify a minimum in the modified space without caring for the modified value itself. We show this in Figure 10, with the modified minimum occurring for . 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 (e.g. top right panel of Figure 9). We note that the total number of single transiting planets expected from the model population for 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 giving the lowest modified 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 50% of observed planetary systems should contain a single planet, with the remaining systems containing multiple planets with mutual inclinations of . 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 = 2 (Chambers (1999)) with the outer of the inner two planets for all 341 considered systems, where and

where 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 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 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 (), semi-major axis () and initial inclination () 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 , and we therefore calculate a modified value described in §6.3. We show these modified values in Figure 11 for an outer planet with = 10 (top panel), = 1M (middle panel) and = 2au (bottom panel). Inclinations of 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 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 space when the outer planet has a semi-major axis of 1au for a mass of 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 = 1M) with a semi-major axis of 1.38au, initially inclined to the inner planets by . Finally for the bottom panel, this minimum occurs for a mass of 6M and inclination of 6 (where 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 .

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 space with = 1.07au, = 24M and =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 . The outer planet parameters which predict are shown by the white lines in Figure 11. This line highlights that while many outer planet parameters can predict , some predict single transiting planets which are more representative of those in the Kepler population. We note that predicted by the same range of outer planet parameters from Figure 11 is shown in Appendix C.

Figure 11: Modified value comparing types of single transiting planets predicted by the model with Kepler population. For the top panel = 10, for the middle panel = 1 and for the bottom panel = 2au. Laplace-Lagrange theory is expected to break down for . The red dashed line refers to a rough RV detection threshold. The white line shows where the model population predicts . The white triangle and circle gives the third planet parameters used to produce the middle and bottom panels of Figure 9 respectively.

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 () is contributed to by different populations of planetary systems. These may include a number of inherently single planet systems () in addition to a number of single transiting planets observed from a population of two planet systems which have a fixed mutual inclination of (). 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 (). Hence in general, it can be considered that

(26)

Here we make the assumption that the total number of double transiting systems observed by Kepler () is made up of a fraction that are two planet systems with an inherent mutual inclination and a fraction (1-) that are two planet systems with an inclined outer companion. We can thus rewrite eq. (26) as

(27)

where () is the number of singles that would have been produced from the population of two planet systems with a fixed mutual inclination of , had it been numerous enough to reproduce the 341 double transiting Kepler systems (which is shown in Figure 10 as a function of ). Conversely () 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 = 447 from §6.3. We note that will change for different values of , however for simplicity we keep it constant at 447.

For the assumed and an assumed fixed mutual inclination for the fraction of the double transiting systems that are inherently inclined (), 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 () and the properties of these planetary systems which determine the ratio of single to double transiting systems from this population (i.e. ()). 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 (). For example, for and , () = 1676 from Figure 12, which from Figure 14 would be reproduced by an outer planet with au, M and = 10. For , () is increased to 2192 requiring the mass of this outer planet to be increased to M (for au and = 10). The outer planet parameters required to produce () are therefore extremely sensitive to the value of . However, increasing the value of for a given value of increases the value of () and hence decreases () as seen in Figure 12, requiring an outer planet which is a weaker perturber of the inner planets.

Figure 12: The number of single transiting planets needed to be predicted by a population of two planet systems with an outer planetary companion, assuming that () of observed Kepler systems host such systems. The remaining fraction of observed Kepler systems are assumed to be two planet systems inherently mutually inclined by .

It should be noted that and 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 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 , 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 (i.e. when ) 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. ) 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 ) 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 , 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, au and , which was in a minimum of the modified 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.

Figure 13: Predicted distribution of mutual inclinations between the two planets in the observed Kepler double transit population for different model populations that both produce the correct number of double and single transiting systems. The grey line refers to the model where the two planet are inherently inclined by =4.4. The black line refers to the model where two planets are secularly perturbed by a outer companion with = 1M, and =1.9au.

From transit duration variation studies, the distribution of mutual inclinations between planets in multi-planet Kepler systems is peaked at (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 to be determined. We look to combine the two models in such a way to predict a value of , 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 , then in order to reproduce a distribution of mutual inclinations that peaks at 2 from modelling of TDVs, it might be expected that .

An additional method to estimate 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 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 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. . From Figure 12 and highlighted in §7.1, if 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 and for 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 to be constrained. Generally for example, a low yield of outer planets in RV studies would suggest that 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 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 , 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 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 () to the total number of stars observed by Kepler (). The individual occurrence rates for the inherently single planet systems is therefore given by (), for the two planet systems with the fixed mutual inclination of by () and for the two planet systems being perturbed by an inclined outer planet by (). For example, for the population of two planet systems which were inherently mutually inclined by 3.6 (for when ), 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, () = 27%. Conversely, considering the population of two planet systems which were perturbed by an outer companion with M, au and (white circle Figure 11 top) for when , predicted 42733 systems in the associated model population. Therefore the associated occurrence rate of this type of system (