© ESO, 2014

Over the past two decades, hundreds of exoplanets have been discovered (Mayor & Queloz 1995; Wright et al. 2011; Schneider et al. 2011)¹, yet only few have been detected as part of binary systems although a large fraction of systems in the Galaxy are composed of multiple stars (Duquennoy & Mayor 1991; Tokovinin et al. 2006; Raghavan et al. 2010). It would be interesting to know if planet formation – and the planets’ subsequent orbital evolution (e.g. Lin et al. 1996; Rasio & Ford 1996; Baruteau et al. 2014; Davies et al. 2014) – happens in a similar manner around single stars and multiple stars. Any difference might offer crucial clues about which formation and evolution scenarios are dominant.

In some cases, planets orbit just one component of the binary, called an s-type orbit, which was highlighted in the recent discovery of Alpha Centauri Bb (Dumusque et al. 2012). In this paper we only consider planets orbiting the centre of mass of the binary, also called a circumbinary or p-type orbit. Schneider (1994) proposed such planets could be discovered using the transit method and, in the course of its photometric survey, NASA’s Kepler telescope discovered a handful of circumbinary systems, all of which are around eclipsing binaries and are nearly coplanar.

Because most binaries do not eclipse, present searches for circumbinary planets have been restricted to a very small fraction of the binary population. In this paper we explore the topic of circumbinary planets transiting non-eclipsing binaries. Whilst the exact number is not presently known, the complete binary population within the Kepler field invariably numbers in the tens of thousands. This is a significant fraction of all observed systems within the Kepler catalog and a large amount of potential planet hosts.

Our study is produced in the context of the Kepler telescope. The high quality of its photometric data and the fact that it performed its observations of the same field for four years provide a controlled sample within which our conclusions can be quickly tested. Our work nonetheless remains easily applicable to past surveys, for example CoRoT, and applicable to future surveys, such as TESS, NGTS, PLATO, and Kepler’s extended mission K2.

This paper is organised in the following way: we first summarise the current circumbinary discoveries (Sect. 1). We then move onto the concept of misaligned circumbinary systems. This is first motivated in Sect. 2 by rapidly reviewing observational and theoretical evidence for mutually inclined planes. Then, in Sect. 3 we introduce the concept of planets transiting non-eclipsing binaries. In Sect. 4 we explore some observational consequences.

In Sect. 5 we detail our method for constructing a synthetic distribution of circumbinary systems, drawing from a realistic distribution of binary stars taken from the Kepler and radial-velocity surveys. We then draw circumbinary planets using a number of possible distributions in mutual inclination and in orbital separation. Based on these distributions, we make predictions on the number of exoplanets transiting both eclipsing and non-eclipsing binaries (Sect. 6), before comparing these results with those found already in the Kepler data (Sect. 7). Finally, in Sect. 8 we discuss our results and ways of finding these planets using current and future telescopes, before concluding.

1. Current circumbinary discoveries

The Kepler telescope has, so far, discovered seven circumbinary systems, one of which has multiple planets (Kepler-47; Orosz et al. 2012b). In addition to the Kepler detections, there have been circumbinary planets discovered by the radial-velocity method (Correia et al. 2005) and through eclipse timing variations (e.g. Beuermann et al. 2010), although there are some doubts about some systems (Mustill et al. 2013).

The Kepler circumbinary planets were discovered by seeking repeated, quasi-periodic photometric dimmings, caused by the planet transiting one of the stars. There is a departure from strictly periodic transits, an effect known as transit timing variations, because of the binary motion. The requirement of semi-regular transits biased these searches to coplanar systems: inclined planets will generally produce transits that are either singular or irregularly spaced, something that we demonstrate throughout the paper.

Whilst the number of detected planets is currently low, two trends are emerging in the orbital distribution of circumbinary planets:

1) All the innermost planets reside close to the stability limit, as defined by Holman & Wiegert (1999), with planets closer to the binary being dynamically unstable (see also Dvorak 1986; Dvorak et al. 1989). The location of that limit is a function of the binary mass ratio and of the planet and binary eccentricities. It can be approximated by P_p ~ 4.5 P_bin. There is a weak dependence on the mutual orbital inclination, with inclined systems generally remaining stable at shorter orbits (Pilat-Lohinger et al. 2003; Doolin & Blundell 2011; Morais & Giuppone 2012).

The origin of this accumulation remains an open question. Inward, disc-driven orbital migration is possible (e.g. Pierens & Nelson 2013; Kley & Haghighipour 2014) and could preferentially place planets near the stability limit. Alternatively, this can result from an observing bias since, for a given binary, we have the highest probability of detecting the shortest stable planetary period. A recent commentary on circumbinary detections by Welsh et al. (2014) suggested that both explanations were possible. We tested both hypotheses, with results in Sect. 7.1 and discussion in Sect. 8.1.

2) There are no claims of planets transiting very short-period binaries (<5 days), although more than half of the eclipsing systems in the Kepler catalog have periods in this range. The cause for this dearth of planets is presently unknown.

The main orbital parameters of Kepler’s circumbinary systems are summarised in Table 1, with P_bin standing for the binary period; and e_bin for its eccentricity; P_p and e_p are the planets’ period and eccentricity; P_crit is the shortest, stable planetary orbit using the criterion from Holman & Wiegert (1999).

2. Observational and theoretical evidence of mutually inclined planes

Observations of circumbinary debris discs show they are largely found to be coplanar (Watson et al. 2011; Kennedy et al. 2012b; Greaves et al. 2014), but some misaligned discs do exist (e.g. 99 Herculis; Kennedy et al. 2012a).

The detections and monitoring of winking binaries, (KH 15D, Winn et al. 2004; WL 4, Plavchan et al. 2008; YLW 16A, Plavchan et al. 2013) indicate that they are composed of a binary surrounded by an optically thick and inclined protoplanetary disc. The disc is warped and precesses, as a result of distance-dependent torques exerted by the binary (Chiang & Murray-Clay 2004). The precession modulates the amount of light obscured by the disc and causes the winking phenomenon (Winn et al. 2006). Some binary stars possess rotation spins inclined with respect to their orbital spin (e.g. Albrecht et al. 2009; Zhou & Huang 2013; Albrecht et al. 2014) and may reflect the same processes that lead to inclined discs.

Proto-binaries formed by fragmentation must have an initial separation larger than 10 AU, with typical separations larger than 100 AU (Bate 2012). Migration of both components is therefore necessary in order for close binaries to form; this can occur under certain disc conditions for an accreting binary (Artymowicz & Lubow 1996). Foucart & Lai (2013) found that during mass accretion, angular momentum gets redistributed leading the binary and disc planes to align. Their corollary is that relatively short-period binaries (sub-AU) will only harbour coplanar circumbinary planets, because they would likely form after the binary migration and thus, after the disc has realigned.

An alternative result was proposed by Lodato & Facchini (2013), who realised that the warped disc can reach a misaligned steady state and subsequently precess as a right body, as it has been assumed by other studies dealing with similar situations (e.g. Terquem et al. 1999). To maintain equilibrium with the binary torques, the farther away the disc is from the binary, the more inclined it needs to be. This has been observationally confirmed (Capelo et al. 2012).

Earlier considerations aside, it is likely that the ultimate orientation of the disc does not completely determine the inclinations of planetary orbits, as epitomised by the existence of inclined and even retrograde hot Jupiters (Hébrard et al. 2008; Winn et al. 2009; Triaud et al. 2010; Schlaufman 2010; Triaud 2011a; Brown et al. 2012; Albrecht et al. 2012). Even if circumbinary planets can ultimately only form in coplanar discs, concluding that planets will also be coplanar does not take into account post-formation dynamical processes, such as planet-planet scattering. These processes likely exist in circumbinary systems just as they do in single-star systems, where scattering is an important mechanism for explaining eccentric and misaligned planets (Rasio & Ford 1996; Adams & Laughlin 2003; Jurić & Tremaine 2008; Chatterjee et al. 2008; Nagasawa et al. 2008; Matsumura et al. 2010a,b). These processes may be enhanced in circumbinary planetary systems. An inner binary induces large eccentricity variations in the orbits of its planets (Leung & Hoi Lee 2013), potentially leading to significant planet-planet scattering. Finally, multi-body dynamics can lead to a wide variety of configurations and we cannot exclude that more exotic scenarios such as an exchange of planets between components of a multiple stellar system would not provide a population of non-coplanar orbits (Moeckel & Veras 2012; Kratter & Perets 2012).

3. Thinking out of the plane

The concept of a planet transiting a single star is well understood and many astronomers have now developed a sort of intuitive knowledge of this technique. Some of these concepts do not directly translate to the case of circumbinary planets, particularly when there is a large mutual inclination between the binary and planet orbits. The purpose of this section is to establish a similar intuition for misaligned circumbinary planets.

3.1. Numerical method

All of the numerical simulations in this paper were conducted using an N-body code that integrates Newton’s equations of motion with a fourth-order Runge-Kutta algorithm. This algorithm generates a slight energy loss over time that is dependent on the size of the fixed integration time-step. Despite this limitation, the fourth-order Runge-Kutta algorithm is sufficient for our purpose: we only required short integrations (up to 200 years) where energy loss is insignificant. As a matter of caution, the energy integral was computed inside each simulation. Its variation was monitored and we made sure it stayed within an acceptable level (generally less than ~10^-7% over a standard four-year integration). The effects of spin, tides, and relativity were assumed negligible over these timescales, and consequently were not included in the model.

3.2. Transitability

To a given observer, a planet orbiting a single star will transit if its orbital inclination is very close to perpendicular to the plane of the sky. The probability that a planet will transit is a simple function of only the scaled radius R_⋆/a, where R_⋆ is the stellar radius, a is the separation between the planet and the star, and we neglect the radius of the planet (Seager & Mallén-Ornelas 2003; Winn 2010). The scaled radius corresponds to a circular area projected on the sky. This probability is purely geometric and is static over time; planets either always transit or never transit². The inclination of the planet with respect to the stellar spin axis does not affect the transit probability. This geometry is depicted in Fig. 1a.

A single star presents an axisymmetric geometry while a binary does not. This implies that estimating the transit probability of circumbinary planets is more complex; one must consider the inclination of the planet with respect to both the plane of the sky (the sky inclination) and the plane of the binary (the mutual inclination). In Fig. 1 we show the four possible transit arrangements for single and binary stars.

In case (b) the planet will transit both stars every orbit. In cases (c) and (d) the planet and binary orbits overlap, but transits will not occur at every period owing to the relative motion of the bodies. Any given transit requires a chance alignment between the star and the planet, which leads us to define the following term:

Transitability is an orbital configuration where the planet and binary orbits intersect on the sky, such that transits are possible but not guaranteed at every planetary orbit.

The concept of transitability leads to a description of each system and suggests when to expect and search for transits. The term is used several times in the paper. Finding planets requires transits, and not just when two orbits are crossing. For instance, both stellar radii have to be taken into account (the planet radius is generally negligible). Finding the fraction of systems that experience a transit is what we will now describe.

3.3. Three classes of transiting circumbinary planets

Throughout the paper we categorise transiting circumbinary planets as follows:

• consecutive transits on eclipsing binaries (EBs consecutive);

• infrequent transits on eclipsing binaries (EBs sparse);

• any transits on non-eclipsing binaries (NEBs).

Planets that consecutively transit eclipsing binaries are defined as presenting at least three transits on the primary star, each separated by approximately one planetary period. This particular class encompasses all circumbinary planets currently found by Kepler, and has a strong bias towards coplanar planets.

The remaining planets transiting eclipsing binaries are labelled as “sparse”. These planets are either misaligned, such that the transit sequence is irregular, or they are coplanar, but with a sufficiently long period that they have not been caught transiting three times within a given observing timeframe.

Planets transiting non-eclipsing binaries are accounted for if a planet transits either star at least once. This is the most relaxed criterion possible. No distinction is made between planets exhibiting consecutive or infrequent transits on non-eclipsing binaries, although consecutive transits will be rare. Confirmation of a planet transiting a non-eclipsing binary, like for most other transiting systems, will require follow-up observations. Our focus will remain solely on gas giants whose transit depths are large enough to be noticed in one event.

3.4. The dynamic nature of circumbinary orbits

Circumbinary orbits are not static. Perturbations from the binary cause the orbital elements³ of the planet to change on quick timescales (typically years for the close binaries we considered), whilst maintaining a stable orbit. This has been calculated analytically (Farago & Laskar 2010; Leung & Hoi Lee 2013) and numerically (Doolin & Blundell 2011). It has also been demonstrated observationally in the Kepler circumbinary discoveries. The variation that most significantly affects the transit probabilities is a precession of the longitude of the ascending node, Ω_p, and the mutual inclination, I_p.

This effect is illustrated in Fig. 2 for an example circumbinary system with binary masses 0.67 and 0.33 M_⊙, period P_bin = 5 days, and eccentricity e_bin = 0.5. The planet was given a period P_p = 37.5 days and eccentricity e_p = 0. A range of starting values for I_p and Ω_p were used, each corresponding to a different line in Fig. 2. On the timescale of years the planet mass has a negligible effect on the dynamics, and consequently was inserted as a test particle, as was done in Doolin & Blundell (2011). The integration was conducted over a full precession period. Our results match Fig. 2 of Doolin & Blundell (2011).

In addition to varying the alignment of the planetary orbit with respect to the binary orbit, precession changes the inclination of the planet’s orbital plane on the sky. This can cause the orbit to come in and out of transitability over time. To emphasise the importance of this effect, we simulated a new circumbinary system composed of a planet (outer red orbit) on a 55 day orbit inclined by 30° with respect to a binary of two 1 M_⊙ stars (inner black orbit) with a period of 10 days (Fig. 3). This system is similar to the existing Kepler circumbinary systems, except with a large mutual inclination. Five snapshots of the osculating Keplerian period are shown in Figs. 3a to e. These are labelled in Fig. 3f, which shows the variation of the planet’s orbital inclination on the sky, over time. The orbit of the planet precesses completely in slightly more than eight years. It enters intervals of transitability when the sky inclination is between 94° and 86°, which corresponds to ~10% of the precession period.

The existing Kepler circumbinary systems have slight but non-negligible mutual inclinations. This means that precession causes them to exit transitability for extended periods of time. Kepler-16b, for example, will start missing the primary star in 2018, after first stopping to transit the secondary in 2014. Calculations show that primary and secondary transits are calculated to return approximately in 2042 and 2049, respectively. A quicker effect is observed in Kepler-413, whose 4.1° mutual inclination is the largest found so far (Kostov et al. 2014).

The systems alternate between being in and out of transitablity, which implies that there are several systems (whose parameters may resemble very closely the current Kepler detections) that remained out of transitability for the duration of the survey. Any attempt to estimate occurrence rates from transiting circumbinary planets must therefore account for the probability that a system is in transitability and for the time it spends in that configuration.

It is important to note that not every misaligned circumbinary planet will pass through a region of transitability. In Fig. 4a we present the same planet that we simulated in Fig. 3, but showing its trace over an entire precession period for a binary plane inclined on the sky by 80°. This graph also shows an inherent property of an inclined circumbinary planet: its probability of exhibiting transitability, in other words, the probability that its orbit crosses the binary’s, is linked to the scaled extent of the planet’s orbit on the sky. This is much larger than the scaled radius of the star in an equivalent circumstellar case. Figure 4b represents exactly the same system with the difference that it is tilted by 50° with respect to the plane of the sky. The planet no longer crosses any region of transitability, and will never transit.

One could wonder whether a planet can exhibit permanent transitability on a non-eclipsing binary. We find that this is possible, but only in an extreme scenario of perfectly perpendicular planet and binary orbits, seen by an observer where I_bin = 0° and I_p = 90°. One might also consider a resonant system where P_p = NP_bin for some integer N, and the orbits are perfectly aligned such that the planet transits on every passing of the binary. Owing to precession and other short-term variations in the orbital elements (e.g. Leung & Hoi Lee 2013), such a configuration is not sustainable.

4. Observational consequences

It is important to understand how an inclined planet population would affect the rate of detectable transiting planets. In this section, hypothetical misaligned circumbinary systems are used for illustration. With Newtonian physics being scalable, the results were easily generalised in Sect. 5 and then applied to Kepler to allow their empirical verification.

4.1. Transit timing

A single planet transiting a single star does so with strict periodicity. In the presence of multiple planets, each contributes what is known as a wobble to the star, which consequently causes each planet to transit with slight aperiodicity (Miralda-Escudé 2002; Agol et al. 2005; Holman & Murray 2005). An additional effect comes from planet-planet interactions that add an extra perturbation on their orbits. The resulting transit timing variations (TTVs) are generally of the order of seconds or minutes, and have become detectable with the advent of the Kepler telescope (e.g. Holman et al. 2010). In the context of circumbinary planets, the main contribution to the TTVs is not from dynamical interactions but from simple geometry.

In the case of a coplanar planet transiting an eclipsing binary, both stars move significantly about their centre of mass. The position of a star at the time of transit varies of the order of twice the binary semimajor axis. This induces a geometric TTV effect , which is generally on the order of days (Armstrong et al. 2013). Perturbations on a planet’s orbit due to planet-binary interactions or with another planet still take place, but are of smaller amplitude.

In the case of a planet inclined with respect to the binary plane, its precessing orbit (Fig. 3) adds an extra complication to the estimation of transit timings. One must consider three important periods:

• the planet period, in which transits can only occur at approximately integer multiples of the planet period;

• the binary period, which affects how much the transit times can deviate from integer multiples;

• the precession period, which determines the frequency and length of transitability (Fig. 3f).

In Fig. 5 we illustrate the transit timing of a circumbinary planet transiting a non-eclipsing binary, extracted from one of our synthetic populations described in Sect. 5. The system is composed of a 35.1 day planet orbiting a 7.1 day binary with masses 1.06 and 0.44 M_⊙. The orbit of the planet, with respect to the binary, is nearly polar: I_p = 102°. Transits, represented by horizontal black lines, occur within regions of transitability (grey area) on the primary (A) and secondary (B) stars. The time spent in transitability for the secondary star is longer because, being less massive, it has a wider orbit.

There is only a 1° difference in the mutual inclination between the central system and the systems on either side. Despite such a small difference, the transit sequences vary significantly. This is due to the relative motion between one of the stars and the planet, which move in different directions on the sky. Any slight changes to the orbital parameters can cause the planet to be early or late enough such that it misses a transit, for instance the transit on the secondary star near 2.5 years; it is present in the system on the left-hand side, but missing in the middle and on the right.

At other times, the same change does not cause particular transits to be missed: there is a pair of transits on the primary, near two years, that is present in all three variations of the system. Whilst the transits still occur on these particular orbits, the exact mid-transit times and durations vary by hours, although this is not discernible on the plot.

By being so sensitive to the orbital parameters, one can obtain tight constraints on a system’s parameters, provided that there are enough transits. Inversely, it implies that a high precision is required to predict future transits accurately.

4.2. Impact of the observing timespan

Unlike the transit probabilities around single stars, which are static, the probability of observing a transit of a misaligned circumbinary system increases with time. There are two principal contributions:

• For circumbinary systems exhibiting transitability, any given period has a transit probability between 0 and 1. Longer observations mean that the planet covers more orbits, and hence has more opportunities to transit.

• The state of transitability is time-dependent because of precession (Fig. 3). Longer observations yield more time spent in transitability, and consequently a greater chance of catching a transit.

To demonstrate how transit probabilities evolve with time, and explore the link with mutual inclination, we took the Kepler-16 system as a template. We reproduced its binary and planet masses, radii, periods, and eccentricities. The angles within the orbital planes, ω_bin, M_bin, ω_p, and M_p, where M denotes the initial mean anomaly, were randomised between 0° and 360°. For the inclinations we drew I_p from three uniform distributions: 0° and 15°, 0° and 30°, and 0° and 45°. The nodal angle Ω_p was randomly drawn between 0° and 360°, since observations could coincide with any point on the precession period. The 3D orientation of the circumbinary system on the sky, with respect to a given observer, was randomised by applying a uniform 3D rotation algorithm by Arvo (1992).

We simulated 10 000 randomly drawn systems over four years, and recorded stellar eclipses and planet transits. Circumbinary systems were then discard as non-transiting or classified as transiting in one of the three categories described in Sect. 3.3. In Fig. 6 we plot the fraction of systems in each of the three classes as a function of time.

The fraction in the category of “EBs consecutive” remains at zero until the planet has completed at least two orbits (1.25 years), because of the requirement of three consecutive primary transits. The small jump corresponds to the small fraction of systems transiting eclipsing binaries. Once three orbits (1.87 years) have been completed, there are no additional detections.

The fraction for the other two categories can increase from zero immediately since we set a requirement of a single transit. Three trends are immediately recognisable:

• There is a higher fraction of systems in the “EBs sparse” class, than in the EBs consecutive class, because a planet inclined with respect to the binary plane precesses, and more systems enter into transitability.

• The fraction of systems whose planet transits an “NEB” is larger than in the EBs sparse category, because of a property that is seen throughout the entire paper: whilst the probability of detecting a transit in any given non-eclipsing system is low, non-eclipsing systems make up the large majority of all binaries (98% in this particular example).

• The total fraction of systems whose planet transits an NEB increases with the spread in orbital inclination. This indicates that binary systems at higher orbital inclination on the sky are reached, increasing the number of considered systems (see Sect. 4.3.1 for details).

We were interested to see what would happen if we extended the simulations of Fig. 6 to a hypothetical 200-year continuous observing run. The results, provided in Fig. 7, indicate that the fraction of transiting systems asymptotes. After such a long observing time, all systems exhibiting transitability have had at least one transit occur: it becomes unlikely that a planet would pass the binary orbit so many times and avoid transiting. The fraction of planets transiting eclipsing binaries is barely discernible from zero on this scale.

The asymptote is the total fraction of misaligned circumbinary systems that could ever be observed to transit, at some remote point in time. It is a significant departure from traditional considerations about transit probabilities around single stars, particularly for planets with periods greater than 200 days. This confirms the interest of defining the concept of transitability: the probability that a given system will go through transitability is static. It is a property linked to a system’s orbital and physical parameters whose only remaining unknown is the inclination of the system on the sky (see Fig. 4). The probability of transit in the case of a planet revolving around a single star is also static.

4.3. The effects of orbital parameters

The current section will explore how different parameters change the dynamical and geometric properties, and ultimately affect the likelihood of transit. To demonstrate this we took the same template as in Sect. 4.2, a Kepler-16 like system. Individual parameters were then altered in specific tests.

4.3.1. Mutual inclination and orbital periods

The parameters with the highest impact on transit probabilities are the mutual inclination between the planetary and binary planes, and their respective orbital periods. In this section we quantify their effects on the expected frequency of transiting systems.

First, for a Kepler-16 analog, we simulated a suite of uniformly drawn mutual inclination distributions whose upper bounds ranged from 0° (strictly coplanar) to 45°, with 5° increments. The results, depicted in Fig. 8a, show the fraction of each of the three types of transiting systems. The simulations were then re-run with the binary and planet periods reduced by a factor of eight (Fig. 8b). Despite some kinks in the plots that reflect statistical noise, our simulations are sufficiently sampled to illustrate some important trends.

Figure. 8a reveals that the consecutive-transit population is maximised in the coplanar case. Its percentage is similar to its equivalent single-star system: 0.43%. For the other two categories there are, as expected, no transits. Moving away from the strictly coplanar case, the frequency of planets transiting consecutively eclipsing binaries decreases and reaches a low plateau. If a planet has a significant mutual inclination it can still transit eclipsing binaries but these systems will more likely be in the sparse category. Except in the coplanar case, most transiting systems are found around non-eclipsing binaries. As in Sect. 4.2, this is due to the large majority of binaries being non-eclipsing. Aside from some statistical noise between 30° and 40°, the number of transiting systems increases with mutual inclination spread, which is consistent with Sect. 4.2.

As one would expect, shortening the orbital periods produces an increase in the fraction of transiting systems (Fig. 8b, we note the higher vertical scale). An additional trend has also appeared. In Fig. 8a consecutive transits of eclipsing binaries are maximised for completely coplanar systems. For shorter period systems in Fig. 8b. the peak instead corresponds to a spread in mutual inclinations. To determine the cause of this, in Fig. 9 we have plotted the cumulative distribution of the sky inclination of the binaries being transited, taken from the simulation in Fig. 8b (P_bin = 5.1 days). For comparison, we have also included the cumulative distribution of all eclipsing binaries in our simulation, irrespective of transits. The binaries eclipse if their sky inclination is between 85° to 95°, and the inclination distribution is approximately uniform. When planets are all coplanar, however, only eclipsing binaries with inclinations between 89° and 91° can be transited consecutively. This is only approximately 20% of all eclipsing systems. When the mutual inclination spread is increased to 5°, the planet is able consecutively transit across eclipsing binaries of all inclinations on the sky. Increasing the mutual inclination spread to an excessive level causes many planets to be too misaligned to transit consecutively.

A spread in mutual inclination can produce transits on the full spectrum of possible eclipsing binaries. A similar effect happens on the probability of a second planet transiting in a system where one is already known to be (Gillon et al. 2011; Figueira et al. 2012).

4.3.2. Binary mass ratio

To test the effect of the binary mass ratio we constructed four test systems, again based on Kepler-16, with primary masses of 1, 1.25, 1.5, and 1.75 M_⊙. The secondary masses were assigned such that the total mass was kept constant at 2 M_⊙. A radius was calculated for each star using the mass-radius relation (Kippenhahn & Weigert 1994), (1)where ξ = 0.57 for M ≥ 1 M_⊙ and ξ = 0.8 for M< 1 M_⊙. The mutual inclination was drawn from a uniform distribution between 0° and 30°. For each test system, 10 000 simulations were run for four years. Transits on both stars in non-eclipsing binaries were counted. The results are shown in Fig. 10a. This process was repeated with binary and planet periods reduced by a factor of 8 (Fig. 10b).

For an equal mass binary the chance of transiting either star is the same, assuming that the radii are the same (Eq. (1)). By varying the mass ratio away from unity, two effects are observed. Transits have an increased occurrence on the primary compared to the secondary simply because the primary star has a larger radius. In Fig. 10b a sharp decrease in transits occurs between mass ratios of 0.33 and 0.14. The precession period increases with decreasing mass ratio (Doolin & Blundell 2011) such that the precession period of the lower mass ratio systems becomes longer than the four years of observations that set the boundary of our simulation. Some systems do not have enough time to move into a window of transitability and consequently, the fraction of transiting systems drops.

Knowing the rate at which planets precess in and out of view is vital for calculating their abundance. The binary mass ratio is a crucial component to estimating the occurrence rates of transiting circumbinary planets for surveys limited in time.

4.3.3. Other trends

There are too many degrees of freedom in a circumbinary system to analyse all of them in depth. The most obvious that we have not investigated in detail is the impact of having eccentric binaries and planets. The eccentricity affects orbital precession (Farago & Laskar 2010; Doolin & Blundell 2011) and stability (Dvorak 1986; Dvorak et al. 1989; Holman & Wiegert 1999). Its effects are, however, included by construction in the population synthesis that will be described in Sect. 5.

4.4. Varied transit shapes

The transit chords of a coplanar planet across an eclipsing binary are constant, and so are its two transit depths. The transit width can vary significantly between consecutive transits, since it is a function of the planet-star relative velocity. A small mutual inclination will make the transit chord vary from transit to transit, because of precession. An example is shown in Fig. 11a, where colour corresponds to the transit time over a four-year simulation (starting black and finishing red).

If there is a significant mutual inclination, the planet is likely to fall into the sparse transiter category. In addition to many transits being missed, the transits that do occur can have significantly different chords across the star, as exhibited in Fig. 11b. Consequently, there can be significant variations in the transit depth, width, and asymmetry.

Moving to the case of non-eclipsing binaries, the only possible transiting planets are misaligned ones. Highly variable transit shapes are expected. An example is given in Fig. 11c for a highly inclined system. Slight curvature is visible in some of the transit chords. The complexity in the transit shapes could make detecting these systems more difficult.

Chords spanning an entire surface, at different angles and for several phases of the stellar rotation would allow us to construct tomographic maps and collect a detailed knowledge of the shape of the stars. This would be particularly interesting for studying tidal deformation in close binaries.

5. Creating synthetic distributions

In the previous section we demonstrated how different orbital parameters affect the fraction and type of transiting systems. We now extend our simulations from a small set of illustrative examples to a complete synthetic population of binary stars (Sect. 5.1) and orbiting planets (Sect. 5.2). This population forms the basis for Sect. 6, where we predict the number of transiting circumbinary planets in the Kepler data, around both eclipsing and non-eclipsing binaries. The methods of forming the population are general enough such that they could be applied to any telescope.

5.1. The binary distribution

Halbwachs et al. (2003) conducted a volume-limited radial-velocity survey of binary stars in the solar neighbourhood, building upon the work of Duquennoy & Mayor (1991). Halbwachs et al. estimated that % of main sequence stars exist in multiple systems, with orbital periods ranging from 1 day to 10 years. We have used their observed distribution to draw our binary systems. Our reasons for choosing a volume-limited radial-velocity survey are that we perceived this sample is well controlled and defined, especially at periods longer than typically found in Kepler.

The distinction between the Halbwachs et al. (2003) and Kepler binaries is that Kepler’s only contains eclipsing binaries, whilst Halbwachs et al. (2003) observed binaries of any orientation. Consequently, Kepler contains a strong geometric observing bias towards short periods, whilst Halbwachs et al. (2003) does not. We artificially applied this observing bias to Halbwachs et al. (2003), to make the two distributions comparable. This is a way to validate our approach, but also to verify for the first time whether the binary star population is similar throughout the Galaxy, by comparing the solar neighbourhood to the Kepler field.

First, we needed to estimate the size of the Kepler sample. We searched the Kepler target catalog and selected all objects having a Data Availability Flag = 2. This is a flag marking all stars that have been observed and logged. As of May 10, 2013, there were over 200 000 entries in the Kepler target catalog, but only 189 440 had measured effective temperatures (T_eff), surface gravities (log g), and metallicities ([Fe/H]). We used the Torres relation (Torres et al. 2010) to estimate masses and radii for each of these entries. We only kept objects whose mass was estimated within: 0.6 <M_⋆< 1.3 M_⊙. This restriction is necessary to coincide with the sample that has been collected by Halbwachs et al. In addition, later, we will need radii estimates and those are better defined using the Torres relation over the range that we just chose. This leaves us with a list of 149 526 systems. Therefore, using Halbwachs et al., we expect 20 186 binaries in Kepler field, with periods between 1 day and 10 years.

We then assigned a mass and a radius to the primary stars by picking 20 186 entries randomly within our selection of 149 526. The secondaries were created using the mass ratio distribution from Halbwachs et al. Their radii were computed using the same mass-radius relation shown in Eq. (1). The need to use different methods to assign the primary and secondary radii led to occasional inconsistencies, where the lower-mass secondary star was assigned a slightly larger radius. In these cases we scaled the primary star radius up to that of the secondary. The binary periods and eccentricities were drawn from distributions presented by Halbwachs et al. We repeated the procedure for the upper and lower 1σ and 2σ errors of the Halbwachs et al. binary abundance.

The 3D orientation was randomised using the same N-body code as earlier. We counted all eclipsing systems, including those that showed only the primary or secondary eclipse. The simulated eclipsing systems were then compared with those in Villanova’s Kepler eclipsing binary catalog⁴, which we downloaded on May 10, 2013. This catalog is an up-to-date version of those presented in Prša et al. (2011) and Slawson et al. (2011). It contained 2400 eclipsing binaries. By cross-referencing with the Kepler target catalog, 2191 of the eclipsing binaries have measurements for T_eff, log g and [Fe/H]. We used the Torres relation to calculate the primary star mass and found that 1684 binaries have a primary mass falling between 0.6 M_⊙ and 1.3 M_⊙.

The simulations cannot be meaningfully compared with the catalog for all period ranges. Because of the nature of Kepler’s observations and its very steep sensitivity with orbital distance, the eclipsing binary catalog contains a large number of very close and contact binaries, many of which have periods less than a day, whereas Halbwachs et al. (2003) contains no binaries closer than 1 day. At longer periods beyond 200 days, the Kepler catalog is likely incomplete since the data analysis and binary identification is still ongoing. We think the most accurate comparison is for orbital periods in the range 5−200 days.

Table 2 contains numbers comparing eclipsing binaries in the Kepler catalog to our simulated systems. Results are presented for several period ranges. Between 5 and 200 days the Kepler catalog almost matches the upper 95% confidence interval predicted by the Halbwachs et al. distribution, but the means differ by about 25%. There are several possibilities for the discrepancy:

• The catalog could contain some false positives due to blended background stars or mis-identification.

• Both methods have different biases; if close binaries are part of triples (as observed by Tokovinin et al. 2006) then the photometry would detect the close pair, while the radial- velocity will be more sensitive to the long period (especially for low mass, synchronously rotating binaries in a triple system).

• Halbwachs et al. (2003) may have surveyed a slightly binary-deficient region of the Galaxy (or Kepler a binary-rich region).

• Our method for assigning and calculating masses and radii might be inaccurate, which would affect the eclipse probabilities (e.g. Dressing & Charbonneau 2013; Batalha et al. 2013).

• The mass range we chose does not exactly correspond the spectral type range defined in Halbwachs et al. (2003).

• The period range used for comparison may not be the best.

A combination of the above is conceivable, and the last three options would be the easiest to manipulate in order to match the eclipse numbers. While attempting to make the numbers compatible would prove to be an interesting study, it is beyond what is necessary for our current exercise.

Further comparisons were made between the two distributions. In Fig. 12 we compare the eclipsing binary data from the simulated and catalog distributions by plotting the cumulative probabilities for the binary period (a); primary mass (b); and primary radius (c). The period and mass distribution match very well over the range 5 to 200 days. This suggests that whilst the raw numbers in Table 2 may contain some discrepancies, the underlying distributions are the same. Our approach complements the upcoming work of Orosz et al. (in prep), who numerically invert the Kepler eclipsing binary catalog to determine the underlying distribution.

According to the plot for the primary radius (Fig. 12c), the Kepler catalog contains a slightly lower abundance of stars above 1.5 R_⊙ than we simulate. The affected population is of the order of 10% and does not have a large impact on our simulations (Figs. 12a and 12b).

As a result of these tests, the binary distribution was considered to be sufficiently realistic to be used in the wider simulations. The lower eclipsing binary counts in the simulated distribution can only act to slightly underestimate predicted planet numbers, something that can be corrected by scaling them up.

5.2. The planet distribution

The binary distribution (Sect. 5.1) could be created with reasonable confidence since binary stars have well studied properties. Circumbinary planets, however, have only been conclusively detected within the past few years. Because of the limited knowledge of how circumbinary planets are distributed, we tested a wide range of possible orbital configurations, constructed from observational and theoretical results.

We dedicate sections to the two parameters that have the largest impact on transit detections: a planet’s inclination with respect to the binary plane, Sect. 5.2.1, and its orbital period, Sects. 5.2.2 to 5.2.4.

We restricted our simulations to gas giants, i.e. planets with masses larger than 50 M_⊕ (~0.15 M_Jup). This choice is not arbitrary, but is based on the fact that the orbital distribution of planets with smaller mass at orbital periods larger than 100 days is largely uncertain (Mayor et al. 2011; Howard et al. 2010, 2012) and would be hard to simulate without making many ad hoc assumptions. Furthermore, large planets cause deep transits that are more easily identified from few events, as expected in the inclined case.

The abundance of planets above 0.15 M_Jup orbiting single, solar-like stars, with periods less than 10 years is 13.9 ± 1.5% (Mayor et al. 2011). We used this abundance in the construction of our distributions.

5.2.1. Mutual inclination distributions

The distribution in mutual inclination has a large effect on the amount of transits observed (Sect. 4.3.1) which is currently an unknown quantity. We built a number of distributions that are described by five mathematical functions: uniform, Gaussian, Rayleigh, isotropic, and coplanar. For the uniform, Gaussian, and Rayleigh distributions, respectively, the upper boundary, standard deviation, and Rayleigh parameter received the values 1, 5, and 20°.

Some of these descriptions are supported by observational and theoretical studies. Foucart & Lai (2013) theorise that circumbinary planets should be very close to coplanar. Planet-planet scattering in single-star systems was tested in a number of numerical experiments by Chatterjee et al. (2008); the resulting inclination distribution for the inner planet (the one we could detect most easily) closely matches a 20° Gaussian distribution centred on zero. Hot Jupiters orbiting stars with mass >1.2 M_⊙ and younger than 2.5 Gyr have spin–orbit angles compatible with an isotropic distribution (Triaud 2011b). Finally, objects of the solar system follow a 1° Rayleigh distribution relative to the invariant plane (e.g. Lissauer et al. 2011; Clemence & Brouwer 1955).

We also tested a hybrid distribution that was constructed by combining two descriptions to take into account that two migration mechanisms can be responsible for shaping the orbital parameters of gas giants (Dawson & Murray-Clay 2013). Disc-driven migration and planet-planet scattering can both send planets within the snow line, but both together may be required to understand the existence of hot Jupiters, notably, their relation to metallicity (Dawson & Murray-Clay 2013) and the spread of their orbital obliquities (Triaud et al. 2010; Winn et al. 2010; Triaud 2011a; Albrecht et al. 2012). We therefore created a new distribution from the combination of a 1° Rayleigh with a 20° Gaussian, represented in equal proportions, that we called hybrid.

5.2.2. A featureless distribution of planets

We started with the simplest possible scenario: a uniform log ₁₀ planet distribution. We assigned planets on 13.9% of the synthetic binaries created in Sect. 5.1. From the flat log ₁₀ distribution we drew the planet period between 2.1 days and 8.6 years, where these bounds came from radial-velocity surveys⁵. The planet eccentricities, were also taken from the same surveys.

We removed any planets that were placed on dynamically unstable orbits, after computing the stability criterion defined in Holman & Wiegert (1999). This significantly reduced the simulated occurrence rate from 13.9% to 3.1%

For single stars, transit probability is proportional to , leading to an easily correctable observing bias towards short-period planets. The circumbinary geometry is more complex because there are two periods to consider: P_bin and P_p. Thanks to our choice of a featureless period distribution, we can investigate the geometric observing biases of circumbinary transiting planets.

5.2.3. Using planets from radial-velocity surveys

Our second set of simulations relies on the period distribution of radial-velocity-detected planets. While we may or may not expect that the circumbinary planet population should occupy a different orbital distribution from the circumstellar distribution, the current distribution is the only one that we know for sure has some basis in reality. Testing it can reveal similarities and distinctions between the populations orbiting single and binary stars.

Radial-velocity surveys have yielded planets with an orbital separation distribution that contains a number of structures, of which the most spectacular is an important and sudden rise in the number of objects beyond ~1 AU separation, dubbed the 1 AU bump. As it turns out, the current Kepler circumbinary detections straddle this bump.

We used exoplanets.org (Wright et al. 2011) to compile a list of planets, detected by the radial-velocity technique, with masses in the range M_p > 50 M_⊕. We only kept planets orbiting stars with masses between 0.6 < M_⋆ < 1.3 M_⊙ in order to be consistent with the synthetic binary distribution. We also removed periods longer than 10 years. We then probabilistically assigned planets to the binaries according to the 13.9% abundance and drew the period and eccentricity from the observed radial-velocity distributions.

After removing planets according to the Holman & Wiegert (1999) criterion, the occurrence rate is effectively reduced from 13.9% to 5.5%. The occurrence rate of gas giants around single stars reflects the efficiency of planet formation. Assuming this efficiency is the same in circumbinary proto-systems nevertheless leads to anticipating 60% fewer planets, most of them being lost within the instability region.

5.2.4. Accumulating planets near the stability limit

Our third and final period distribution was crafted to match the currently observed location of circumbinary planets. We use here the hypothesis that the accumulation of Kepler-discovered planets near the stability limit has a physical origin⁶. Some scenarios, including a halted disc-driven migration at the disc’s inner cavity, are discussed in Sect. 8.1.

In the previously simulated set-up of Sect. 5.2.3 we removed any planet that was placed within the instability region. We created a new set-up, similar to Sect. 5.2.3, with the exception that we rescue a fraction of otherwise lost systems. To do so, we reassigned 50% of the planets placed within instability, to a new, stable, near-circular orbit with a period randomly selected such that P_p/P_crit is between 1.1 and 1.43. These bounds correspond to the RMS of P_p/P_crit for the observed innermost circumbinary planets (Table 1). Rescuing these planets simulates a sort of stopping mechanism that would prevent some but not all planet loss.

In this pile-up scenario, the circumbinary planet occurrence rate fell from 13.9% to 9.5%.

5.3. Simulation set-up

Two binary populations were created that simulate the contents of the Kepler field, composed of 21 553 and 18 275 systems. For each binary distribution ten sets of planets were created. This was reproduced for all three planet period distributions. This means our results come from the combination of 20 different sets of circumbinary systems for each planet period distribution (i.e. 60 sets in total). For each of these 60 sets, the mutual inclination between the binary and planetary orbits was simulated using the 12 misalignment distributions described in Sect. 5.2.1 and presented in tables in the Appendix.

Finally, the 3D orientation of the circumbinary system in space was randomised. Using the N-body code, these systems were then simulated over a four-year timespan. All stellar eclipses and planetary transits were recorded.

6. Results

6.1. Simulated planet abundances

In Sects. 5.2.2 to 5.2.4 the effective planet abundances before the start of our simulations were 3.4%, 5.5%, and 9.7% for the log ₁₀, radial-velocity, and pile-up periods distributions. Our final numbers are affected by two additional factors. First, after running the N-body simulations a small number of unstable systems were discovered and discarded from our statistics. We missed them earlier because the stability criterion of Holman & Wiegert (1999) does not include the planet’s eccentricity. Second, we removed all binaries with periods below 5 days. Our justification for this is that no planets were found transiting Kepler’s very short-period binaries, where the a priori probability that they would transit is largest. We expand on this in Sect. 6.4. The final fraction of systems with gas giants are

• uniform log ₁₀ distribution: 2.5%;

• radial-velocity distribution: 4.5%;

• physical stability limit pile-up distribution: 9.0%.

6.2. Predicted transiting planets

In Tables 3 to 5 we show the predicted transiting planet numbers for each of our three period distributions. We give the mean number that we find across the 20 runs (including two binary populations) that compose each table entry. We also provide the minimum and maximum numbers of transiting systems across any of those runs. Results are separated in the three subclasses that we defined in Sect. 3.3. We further combine the EBs consecutive and EBs sparse systems in a broader bin enclosing all planets transiting eclipsing binaries⁷.

The most interesting result is that in most scenarios, and regardless of the period distribution we used, planetary transits are expected in non-eclipsing binary systems. Furthermore, it is striking to see that they may be many more of them than in the eclipsing case. This bodes well for the search for inclined circumbinary planets.

The number of planets expected to transit non-eclipsing binaries increases with the extent in the mutual inclination distribution. The number of transits in the EBs consecutive category generally decreases with the inclination spread, but the maximum does not necessarily correspond to the coplanar simulation. The number of transits in the EBs sparse category does not have a monotonic relationship with the inclination spread. All of these behaviours were announced in Sect. 4.3.1.

Current examination of the Kepler data focuses solely on the EBs consecutive category. However, we note that a planet occurrence rate cannot be securely estimated from this category alone: there is a degeneracy between the fraction of systems with planets and the underlying distribution in mutual inclinations.

The widest distribution in planetary orbital inclination will produce some systems that are near coplanarity and that have a significant chance to transit and produce consecutive events. Their presence does not say anything about all the other systems that could not be detected in transit. This means that a given number of detections is equally consistent between 1) a rare population of mostly coplanar circumbinary gas giants and 2) a population that is widely distributed in inclination and that has a generally higher planet frequency.

Only a minimum occurrence rate can be inferred from just considering planets transiting eclipsing binaries in a semi-regular fashion. Understanding and searching for planets transiting non-eclipsing binaries is required to help break this degeneracy.

6.3. Multi-transiting planets

Knowing that a planet transited once and that the system is a non-eclipsing binary (from spectroscopic follow-up for example) can be sufficient to know that the planet is on an inclined orbit. However, we would be left without much additional information about the orbit occupied by the planet. Because of the binary motion, the transit width, which can be used to obtained a vague idea of the orbital parameters on single stars (e.g. Dawson & Johnson 2012), becomes a degenerate parameter for binaries. Having two transits would help in many respects, from the identification of the system to its characterisation. Exploring our synthetic population, we estimated the fraction of systems that experience two transits compared to the total number of circumbinary transiting planets. Results are given in Table 6. Approximately 50% of systems that have transited once will transit at least a second time⁸.

6.4. Planets around very short-period binaries

The current lack of planets around short-period binaries (P_bin< 5 days) is considered surprising for two reasons. First, the median binary period in the Kepler eclipsing binary catalog is 2.9 days (as of May 10, 2013). Second, for a given ratio of planet and binary periods, the transit probability is maximised for a very close binary (Sect. 4.3.1). The inclusion of all binaries with periods under 5 days causes a significant increase in the amount of expected transiting planets. The number of transiters on eclipsing binaries more than doubled the results described in Sect. 6.2, whilst there was a ~50% increase in planets transiting non-eclipsing binaries. This result reaffirms the common thought that these planets should be numerously detected, if they exist at all.

7. Comparisons with Kepler

The Kepler circumbinary gas giants all fall in the EBs consecutive category. We can use our predictions and see if they reproduce existing discoveries. This comparison is able to distinguish between the three orbital distribution we investigated, to find that one is more likely to be true. This allowed us to reach a minimum occurrence rate compatible with current detection levels.

7.1. On the planet distribution

We chose the 1° Rayleigh mutual inclination distribution to compare our three proposed period distributions to Kepler’s discoveries. Our choice is based on the premises that 1) strict coplanarity is ruled out by the Kepler systems and that 2) the 1° Rayleigh distribution produces a relatively high number of transiting planets in the EBs consecutive category, which is the only type of system that can be compared with Kepler.

The empirical data and our simulations are compared to each others using the ratio P_p/P_crit. To improve our statistics, we increased the number of simulations from 20 (used for Sect. 6) to 1000. As in Sect. 6, only binaries with periods longer than 5 days were included. Their respective cumulative probability distributions are displayed in Fig. 13. The empirical cumulative distribution of P_p/P_crit, only includes the four Kepler discoveries that most closely match our tested mass range: Kepler-16b, -34b, -35b and -64b. A Kolmogorov-Smirnov (K-S) test was used to compare the simulated and observed systems; we tested the relative shapes of the distributions and not the planet’s numbers. Results of the K-S test are shown in Table 7.

7.1.1. The flat log ₁₀ period distribution

A uniform distribution in period is a poor match to the Kepler observations (Fig. 13a). This implies that the observed distribution of circumbinary gas giants is not featureless, but has some structures. Notably, observing biases appear to be insufficient to explain the stability limit pile-up in the Kepler results. This points to a physical origin for this over-density. Some possible explanations are discussed in Sect. 8.1.

7.1.2. The radial-velocity period distribution

The radial-velocity distribution of planets has an even poorer match to the Kepler data (Fig. 13b). We can draw from this that the circumbinary gas giant planet period distribution is inconsistent with the distribution of gas giants orbiting single stars. This is not surprising, since planet formation around single and binary stars is likely different enough to result in distinct distributions of planets. However, this had not been empirically demonstrated before.

7.1.3. The pile-up period distribution

The planet period distribution with a pre-imposed stability limit pile-up is by far the closest match to the Kepler data (Fig. 13c). Indeed, it is 120 times more likely to be true than the uniform, log ₁₀ distribution, and is more than 400 times better compared to the distribution of planets based on radial-velocity results. This comparison to empirical data leads us to conclude that the observed accumulation of gas giants near the stability limit is more likely to have a physical origin than to be an observing bias.

Despite its higher K-S test score, this distribution remains only marginally compatible with the known systems: our simulation predicts that we should find a higher fraction of transiting planets far from the stability limit than has been observed. This suggests that our mechanism, rescuing 50% of planets to populate the pile-up may not be sufficient, that the accumulation is likely to be more pronounced that what we set it to be.

Figure 14 illustrates our results. It is a 2D histogram of the planet and binary periods for systems in the EBs consecutive category only. Darker tones represent a higher probability. At the top and right are the marginalised probabilities. Like previous results shown, results for P_bin< 5 days are excluded. We delineated the approximate stability limit boundary at P_p/P_bin = 4.5. The innermost planets in each of the Kepler systems are overplotted, and approximately follow the relation P_p/P_bin ≈ 6.

Transits are more numerous when the planet and binary periods are short, as expected, but this is for two unrelated reasons:

• shorter period binaries are more likely to eclipse;

• shorter period planets are more likely to transit.

The circumbinary geometry does not make us preferentially observe planets at any special ratio of binary and planet periods.

7.2. Planet occurrence

Strictly within the bounds of our synthetic population (P_bin> 5 d; 0.6 <M_⋆< 1.3; M_p> 0.15 M_Jup)⁹, Kepler has found two systems: Kepler-16 and -34.

In Sect. 7.1, we demonstrated that the planet distribution with a pre-imposed stability limit pile-up is by far the closest match to the Kepler observations. From now on, we will only consider this distribution. The 1° Rayleigh inclination distribution predicts the most likely number of discoveries in the EBs consecutive category at 2.05 systems (Table 5). This is evidently compatible with Kepler’s two discoveries.

From these results, we can affirm that Kepler’s observations are compatible with a 9% abundance of circumbinary gas giants, whose distribution in period has a physical pile-up.

It is likely that this abundance is an underestimate, for the following reasons:

• Our simulated eclipsing binary numbers are on the low side.

• The Kepler discoveries are likely not exhaustive (e.g. Kostov et al. 2014).

• Two more systems (Kepler-35 and -64) lie just outside the range of our synthetic population, which if included would inflate our number.

• There is a degeneracy between abundance and mutual inclination spread, described earlier in Sect. 6.2, that can only lead to an increase in the occurrence rate.

Based on the physical parameters of three discovered systems, Welsh et al. (2012) had inferred early estimates of the occurrence rate between 1% and 10%. Recently, Armstrong et al. (2014), who used an independent method of analysing the Kepler light-curves to constrain the abundance, found a compatible estimate. It will also be interesting to see how our value compares with future studies.

8. Discussion

We conceptually and numerically explored what an inclined circumbinary planet implies for its probability to transit. Using this knowledge we constructed a series of simulations that built a synthetic population of circumbinary systems. We tested three planetary period distributions, on 12 different mutual inclination distributions for two binary populations. The reason we computed so many different alternatives was to tabulate enough results that can then be compared in the possible detection of planets transiting non-eclipsing binaries. The period distributions we tested were either mathematical or derived from observational evidence. The purpose of these choices was to make the work as model-independent as possible. It would, however, be interesting to repeat this work based on a theoretical population synthesis similar to those created for single stars (Benz et al. 2014, and references there-in).

The comparison of our distributions to empirical data collected by the Kepler spacecraft indicates that there is a physical pile-up of planets close to the stability limit as defined by Holman & Wiegert (1999). The number of detections implies that at least 9% of binaries host a circumbinary gas giant with an orbital period shorter than 10 years.

Once it is estimated that the Kepler data has been exhaustively analysed, around both eclipsing and non-eclipsing binaries, ratios between inclined and coplanar systems can be computed and compared to our results. Whether inclined planets are found or not would be interesting regardless, as it would constrain some formation and evolution models and place a limit on the maximum occurrence rate.

8.1. Implications for planet formation and evolution in circumbinary systems

Stellar formation produces binaries with separations usually larger than 10 AU whose separation can then decrease (e.g. Bate 2012). This orbital rearrangement could have been detrimental to the production of planets, yet as we have shown, reaching Kepler’s number of detections requires effective occurrence rates close to those found for single stars. From this, we can conclude that gas giant formation is probably as likely to take place in circumbinary discs as it does in circumstellar environments in opposition to some theoretical expectations (e.g. Meschiari 2012b). This may, however, exclude systems involving binaries on periods shorter than five days whose likely, late dynamical origin (Mazeh & Shaham 1979; Tokovinin et al. 2006; Fabrycky & Tremaine 2007) would probably lead to the disruption of any planetary system that would have formed earlier.

If gas giant formation rates are similar in binaries and single stars, the location of formation, or the planets’ subsequent orbital evolution appears different. Despite low numbers of detections, we are able to demonstrate that the observed over-density of objects found near the stability limit is difficult to explain by an observing bias selecting only the shortest possible planet period for a given binary period. We conclude that there is more likely a physical pile-up, which is reminiscent of the hot-Jupiter pile-up whose origin remains hotly debated (Kuchner & Lecar 2002; Eisner et al. 2005; Ford & Rasio 2006; Matsumura et al. 2010a; Guillochon et al. 2011; Dawson & Murray-Clay 2013; Plavchan & Bilinski 2013). There is, however, one important distinction: if the over-density of hot Jupiters can be explained by a special orbit towards which highly eccentric planets circularise (Ford & Rasio 2006; Matsumura et al. 2010a; Guillochon et al. 2011), the necessary tidal friction is absent in the circumbinary case.

Meschiari (2012a) studied planetesimal accretion in a circumbinary disc. His results imply that planet formation is likely inhibited between 1.75 AU and 4 AU around Kepler-16 meaning that its planet either has migrated (e.g. Goldreich & Tremaine 1979, 1980; Lin et al. 1996; Ward 1997; Rasio & Ford 1996; Baruteau et al. 2014; Davies et al. 2014) or has formed in situ. In situ formation would probably reproduce fairly well the pile-up that we observe however, gas giant formation inside the snow line requires specific conditions (Bodenheimer et al. 2000) and is not expected to take place in circumbinary discs (Paardekooper et al. 2012).

At this point it would seem more likely that the Kepler discoveries have formed beyond 4 AU and migrated in. Disc-driven migration was first studied for circumbinary planets by Nelson (2003), which, with more recent work (Pierens & Nelson 2007, 2008, 2013), shows that planets can migrate inwards before being halted near the inner cavity of the disc. The cavity’s location is close to the stability limit for planetary orbits. The wtork of these authors is consistent with the pile-up observed by Kepler, but is unable to accurately fit the circumbinary orbits of specific systems. The independent study of Kley & Haghighipour (2014), using models containing more sophisticated thermodynamics and boundary conditions, successfully reproduces the observed orbit of Kepler-38b, near the stability limit.

Gas giants could also migrate via an exchange of angular momentum with other planetary companions. From this postulate it is natural to expect that there should be inclined planets since planet-planet scattering has been shown to be an important mechanism for explaining the eccentricity and inclination of gas giants (Jurić & Tremaine 2008; Chatterjee et al. 2008; Nagasawa et al. 2008; Matsumura et al. 2010b). The inclination distribution may be enhanced by the fact that circumbinary discs are more likely to be inclined with respect to the inner binary at large separations than close-in (Lodato & Facchini 2013; Plavchan et al. 2013).

Instead of involving one orbital distribution, a proper expectation of circumbinary planets may involve two. For instance the presence of a pile-up at the stability limit could be linked to the coplanarity of the system. In this case there would be no reason to expect that inclined planets would also show a pile-up.

Undoubtedly more detections are needed to better understand planet formation and orbital evolution in the circumbinary context.

8.2. Towards the detection of planets transiting non-eclipsing binaries

We have shown that only a minimum occurrence rate can be extracted from only considering systems with semi-regular transits in a search of the Kepler data. To place a limit on the maximum occurrence rate, we must find if any planets transit non-eclipsing binaries, and if so, how many. Here, we explore potential search strategies to be applied in the search for these as yet non-considered planets.

One strategy is to first search for non-eclipsing binaries, which can be spotted from a variety of sources. Reconnaissance spectra of the Kepler field may already show that a number of targets are double line binaries. Non-eclipsing binaries can also be identified thanks to their photometric variability that combines ellipsoidal, reflection, and Doppler beaming effects (e.g. Faigler & Mazeh 2011; Shporer et al. 2011). This method is sensitive to Kepler binaries for periods as large as 20 days, and has recently been improved to include eccentric orbits (Mazeh, priv. com.). In the case of high eccentricities, binaries can exhibit a heartbeat-like photometric signature; over 100 of these have already been noted in the Kepler data. Once a binary is confirmed, a search of the light curve can be examined for transit like signatures. The quality of Kepler’s photometry implies that the transit of a Jovian planet (the objects we made predictions for) can be found visually without the need for a sophisticated algorithm.

Alternatively, one could search the Kepler data without any a-priori binary knowledge, looking for either single or infrequent transit-like signatures. The standard Kepler pipeline is not sensitive to these as it requires quasi-periodicity. This type of search may be appropriate for the Planet Hunters team¹⁰ (Fischer et al. 2012), which has already had success finding a circumbinary planet (Schwamb et al. 2013). Algorithms with the ability to identify single transit events would tremendously help the search. Although they would also likely pick up many false astrophysical positives (Brown 2003) and data artefacts. This risk can be reduced by insisting that at least two events must be present in the data; we have shown that about 50% of systems will do so. After finding some candidates in the photometry, a spectroscopic follow-up would easily reveal most single- or double-line binaries (e.g. Santerne et al. 2012; Díaz et al. 2013) fairly inexpensively. The width and depth of a circumbinary planet are often different depending on the orbital phase of the star that is being crossed. A transit whose shape implies a mean stellar density (Seager & Mallén-Ornelas 2003) compatible with either of the two components would become an interesting candidate for additional follow-up. Other usual techniques, such as searching for potential blends using adaptive optical systems, or using statistical algorithms like Blender (Fressin et al. 2013) can help remove many false positives.

If enough events are captured and follow-up constrains the physical parameters of the system sufficiently, it could become possible to predict when the next transit events would happen. A hit or a miss can refine the parameters to impressive precision as illustrated in Fig. 5. With predictive capacity, transits can be recorded, for example, to study the atmospheres of fairly cool gas giants (e.g. Seager & Deming 2010).

The above strategies can be immediately applied to the Kepler telescope since the data is already freely available. They can also be applied to future photometric surveys. Transits of non-eclipsing binaries may take on a larger importance in NASA’s upcoming TESS mission. The short observing timespan of TESS (one month for most of the sky, one year at the ecliptic poles, Ricker, priv. comm.) means that the detection strategy of circumbinary planets will have to evolve from that of Kepler’s, since getting multiple transits of circumbinary planets will be rare, even when coplanar. Sparse and single transits will need to be identified in the light curves. Confirmation would almost certainly require follow-up spectroscopy, which will be made easier by TESS’s focus on bright targets across both hemispheres. Identifying non-eclipsing binaries within the field will also become possible thanks to ESA’s Gaia satellite or even prior to launch from dedicated ground-based surveys. It is possible that inclined circumbinary planets will provide the bulk of transiting cold Jupiters in the TESS program. Additionally, the time-dependent transit probabilities outlined in Sect. 4.2 support an extended TESS mission, which is already under consideration. The proposed Kepler follow-up mission K2 is similar to TESS, but on a smaller scale.

In addition to the interest of detecting more circumbinary planets, there is the appeal of knowing whether misaligned systems can exist in orbiting binaries, like they do around single stars. Several techniques have been developed to measure the spin–orbit angle: spectroscopically, via the Rossiter-McLaughlin effect (Holt 1893; Rossiter 1924; McLaughlin 1924; Queloz et al. 2000; Gaudi & Winn 2007; Boué et al. 2013) and the Doppler shadow (Collier Cameron et al. 2010); statistically (Schlaufman 2010); photometrically thanks to gravity darkening (Szabó et al. 2011) or through spot-crossing events (Nutzman et al. 2011; Sanchis-Ojeda et al. 2012); astrometrically (Sahlmann et al. 2011); using asteroseismologic mode splitting (Chaplin et al. 2013); or interferometrically (Le Bouquin et al. 2009). Searching for transits on non-eclipsing binaries would become another way to detect inclined orbits by spatially resolving the misalignment.

9. Conclusions

We started this paper as an enquiry into whether any circumbinary planets could be caught transiting a non-eclipsing system. After studying and describing how the dynamics and how various physical and orbital parameters affect the probability of transit, we constructed synthetic populations of circumbinary systems. By simulating the orbits over 4 years and checking for transits, we reveal that Kepler could have caught such planets.

We urge a change in the protocols behind the search for circumbinary planets. Efforts should not just focus on eclipsing binaries since most binaries do not eclipse. We also find that applying a quasi-periodic transit signal criterion to find and confirm these planets severely restricts the number of objects that can be found, and biases them towards coplanarity.

Our simulations allow us to affirm that the binary mass and period distributions are similar in the solar neighbourhood and within in the Kepler field. This indicates the formation of binaries might be a universal process, probably linked to the initial mass function (e.g. Salpeter 1955; Kroupa 2001; Chabrier 2003; Bate 2012). We also find that it is likely that circumbinary planets are rare for binaries with periods <5 days lending further support for a dynamical origin for the closest of binaries (Mazeh & Shaham 1979; Fabrycky & Tremaine 2007; Tokovinin et al. 2006).

Comparing the existing detections with the output of our simulations, we remark that the occurrence rate of circumbinary gas giants may be close to that of single stars. Current models of gas giant formation (e.g. Pollack et al. 1996; Boss 2000; Helled et al. 2014) produce planets far from their star. To a first order it would seem that whether the central object is a singleton or a close pair does not have a large impact on planet formation efficiency. Inclined orbits are frequently found in single-star systems and have highlighted the important role dynamical interactions have in shaping exoplanetary systems. A similar formation rate would suggest similar mutual interactions happening in circumbinary systems. This further supports our claim that there should be a number of planets transiting non-eclipsing binaries waiting to be discovered within the Kepler timeseries.

The location of the current detections reveals that after formation, circumbinary gas giants have a migration history that is different from single-star systems, with planets likely piling-up at orbits near their binary’s stability limit. It remains to be seen whether disc-driven migration could place planets on such particular orbits without pushing them onto an unstable orbit, or whether a different mechanism, maybe specific to binaries, might be at work.

---

Acknowledgments

We would like particularly like to thank Dave Armstrong, Simchon Faigler, Veselin Kostov, Rosemary Mardling, Tsevi Mazeh and Stéphane Udry, for insightful and important discussions. We would also like to thank the many people who commented on the concept, or on specific parts of the manuscript. Those people include researchers working on exoplanets at MIT, notably Joshua Winn, as well as colleagues met during the recent second Kepler Science Conference at NASA’s Ames Research Center. We also acknowledge the work of two anonymous referees who helped improve the paper with their comments. D. V. Martin is funded by the Swiss National Science Foundation. A. H. M. J. Triaud is a Swiss National Science Foundation fellow under grant numbers PBGEP2-145594 and P300P2-147773. We acknowledge frequent visits of the exoplanet.eu (Schneider et al. 2011) and exoplanets.org (Wright et al. 2011) websites as well as an extensive use of the ADS and arXiv paper repositories. Finally, our work could not have happened without the important work and dedication of the Kepler team who compile updated and publicly available candidate lists and stellar catalogs on the MAST (https://archive.stsci.edu/) repository and at Villanova University.

References

All Tables

All Figures