The Number of Transits Per Epoch
for Transiting Misaligned Circumbinary Planets

As discussed in Section 1, this paper primarily concerns - how many transits per epoch does one expect for transiting misaligned circumbinary planets? The number of transits per epoch is perhaps the simplest observable one might reasonably define for such systems, yet we lack a good understanding of its behaviour. Figure 1 graphically illustrates how a circumbinary planet can cause zero, one or two transits. Indeed, more than two transits is also possible if the binary motion is sufficiently fast to “catch up” with the planet again during a given epoch. This case is challenging to graphically illustrate, and thus is not shown explicitly in Figure 1, but it is fully considered in what follows.

In light of this, a series of transit epochs will manifest a potentially variable number of transits each time. For example, the first epoch might yield two transits, the next epoch yields zero, the third epoch just a single dip, and so on. Such a series could be compactly represented using the digits “201…”. From the observer’s perspective, the true nature of the system is of course not immediately obvious; they are simply faced with a particular sequence and need to try to interpret the underlying cause.

The “Random Transiter” (Rappaport et al., 2019) again provides an analogous example of this situation. Rather than seeing a sequence of clustered events of varying counts (as is the case for circumbinary planets), Rappaport et al. (2019) instead report non-clustered transits scattered across the time series. Despite this difference, the problem is ostensibly similar to ours in that the lack of any clear pattern in the transits was confounding and unusual. A basic conclusion of their work is that the sequence was consistent with being produced from a random number generator. If indeed the transit times are truly stochastic, then it would be impossible to ever predict when the next transit would occur based on the previous data. This determination of randomness, if verified with future data, is important because one might hypothesise that complex, chaotic systems (e.g. stellar activity) to feasibly produce (outwardly) stochastic events like this. Thus, the determination of whether a given series of astrophysical observations are random or not is a crucial lens through which we can make progress in interpreting such data.

To this end, we here begin by asking if for circumbinary planets the transit counts per epoch, represented as a sequence of integers, yields a sequence that would (at least outwardly) appear random or not? To address this, we develop and test a memory slot prediction algorithm (see Section 2.4) and show that it outperforms random predictions for non-isotropic circumbinary planet populations. This establishes that non-isotropic circumbinary planet populations are non-random.

In order to investigate the randomness (or lack thereof) for circumbinary transit sequences, we first need a large, correctly-labelled population of such systems to study. Due to the limited number of detected circumbinary planets, and the need for accurate labels, we generated a population of artificial circumbinary planet systems by specifying parameters of the binary and the planetary orbit. In what follows, we describe how we generated each system.

We first randomly pick an observed Kepler target from the Data Release 25 (DR25) Kepler Input Catalog queried using the Mikulski Archive for Space Telescopes (MAST). The queried stellar mass is adopted as the true mass of the primary star in our circumbinary system. In this way, the mass of our primary stars are representative of the mass distribution of the Kepler sample.

Next, the mass ratio, $q$, of the binary stars is assumed to follow a $q^{0.5}$ distribution following (Parker & Reggiani, 2013). After assigning $q$, we are then able to trivially determine the mass of the second star. The stellar radii (which have little effect on our results) can then be approximated using the empirical mass-radius relation of main sequence stars given by Demircan & Kahraman (1991):

In what follows, it assumed that the gravitational influence of the planet on the binary is negligible and thus the planet is treated as a massless particle of one Earth radius. Further, it is assumed that the eccentricity of both the binary and the planetary orbit is zero.

The binary separation is randomly drawn from a log-uniform distribution between a predefined minimum and maximum value. As a starting point, the minimum binary separation among all the observed circumbinary systems is 0.08 AU and we adopt a minimum of 0.05 AU in what follows. For the maximum, we are only interested in stable planets around binary stars. Thus, the maximum binary separation is defined to be the separation such that a planet just outside the instability region defined by Holman & Wiegert (1998) has a period of one Kepler window (4.35 years). With a larger binary separation, even the innermost stable circumbinary planet cannot be detected by Kepler because its period is too long.

For the semi-major axis of the planet, we again set a minimum and maximum possible value and sample in between following a log-uniform distribution. The minimum semi-major axis of a stable planetary orbit is given by Holman & Wiegert (1998) as:

where $a_{\mathrm{bin}}$ is the semi-major axis of the binary orbit, $e_{\mathrm{bin}}$ is the eccentricity of the binary orbit (in our case fixed to zero), and $\mu=\frac{m_{B}}{m_{A}+m_{B}}$, where $m_{A}$ is the mass of the more massive star, and $m_{B}$ is the mass of the less massive star. On the other hand, the maximum semi-major axis is defined such that the corresponding planetary period is equal to one Kepler window. Beyond this maximum semi-major axis, the planet would be unobservable.

The binary star’s orbit should be distributed isotropically (i.e. no preferred orientation) because there is nothing special about the observer’s line of sight. To generate such an isotropically distributed population of binary star orbits, we take caution in specifying the binary orbit’s longitude of ascending node, $\Omega_{\mathrm{bin}}$, and inclination angle, $i_{\mathrm{bin}}$ (defined with respect to the observer’s line of sight).

First, consider the inclination angle. An isotropic population might naively appear to follow a uniform distribution, but this would in fact produce vectors non-uniformly distributed across the celestial sphere (due to the curved geometry). Instead, an isotropic distribution is truly described by a uniform distribution in $\cos i_{\mathrm{bin}}$. On the other hand, the longitude of ascending node of the binary orbit, $\Omega_{\mathrm{bin}}$, is sampled uniformly to yield an isotropic distribution. To ensure that we cover all possible orbital configurations, we uniformly sample $\cos i_{\mathrm{bin}}$ between -1 and 1 and uniformly sample $\Omega_{\mathrm{bin}}$ between $-\frac{\pi}{2}$ and $\frac{\pi}{2}$.

With the binary orbit’s longitude of ascending node, $\Omega_{\mathrm{bin}}$, and inclination angle, $i_{\mathrm{bin}}$, sampled as prescribed above, we now considered these two parameters for the planetary orbit - $i_{P}$ and $\Omega_{P}$. First consider $i_{P}$. The inclination of the planetary orbit with respect to the binary orbit (rather than the observer’s line of sight) characterise the coplanarity of the system, which is intimately related to the expected transit sequences. For example, if a population of circumbinary planets were perfectly coplanar with their host binary’s orbit, one would expect to get transits over every epoch if the system is favorably aligned in an "edge-on" orbit with respect to the observer. If, on the other hand, a population of circumbinary planets were perfectly isotropic (i.e. no preferred orientation) with respect to their host binary’s orbit, then one would expect transits to be less frequent and more erratic in nature. Next, consider $\Omega_{P}$. Whether or not $\Omega_{P}$ is related to coplanarity as well depends on how it’s defined. If we define $\Omega_{P}$ with respect to the binary orbit, which means we always rotate the planetary orbit around the axis perpendicular to the binary orbital plane by $\Omega_{P}$, then $\Omega_{P}$ is unrelated to coplanarity. That is, the orbits are coplanar for any value of $\Omega_{P}$ as long as $i_{P}$ and $i_{\mathrm{bin}}$ are identical. However, If we define $\Omega_{P}$ with respect to the observer’s line of sight, then the difference between $\Omega_{P}$ and $\Omega_{\mathrm{bin}}$ is closely related to coplanarity, where the case with $\Omega_{P}=\Omega_{\mathrm{bin}}$ corresponds to a coplanar orbit. To avoid confusion and for simpler implementation, we choose the second definition, which is to define $\Omega_{P}$ with respect to the observer’s line of sight. Finally, we note that the "coplanarity" being to referred to here is a characteristic of the population of circumbinary planets, not of an individual planet. To generate such a population, we require a flexible, analytic probability distribution describing $\Omega_{P}$ and $i_{P}$, from which one can sample.

As mentioned in the earlier discussion on $i_{\mathrm{bin}}$, an isotropic distribution is truly described by a uniform distribution in the cosine of the inclination angle. The same holds true if we want to produce an isotropic population of planets. However, as discussed before, it is the relative inclination angle that characterises coplanarity. For simplicity of notation, we denote the relative inclination angle as $\Delta i$, where $\Delta i=i_{P}-i_{\mathrm{bin}}$. A coplanar population is described by a highly compressed (approaching a Dirac Delta function) distribution centred on $\mathrm{E}[\Delta i]=0$, or equivalently, $\mathrm{E}[\frac{\pi}{2}-\Delta i]=\frac{\pi}{2}$ and $\mathrm{E}[\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}]=0$. On the other hand, we take a uniform distribution in the cosine of $\frac{\pi}{2}-\Delta i$ for an isotropic distribution, where we use the angle $\frac{\pi}{2}-\Delta i$ for consistency. Since a truncated normal has the flexibility to produce highly compressed, diffuse and even quasi-uniform distributions, we adopted this form in what follows. Accordingly, one may write that $\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}\sim\bar{\mathcal{N}}[0,\sigma_{\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}}]$, where $\bar{\mathcal{N}}$ denotes a truncated normal distribution (with support from $-1$ to $+1$ and centered at $0$) and $\sigma_{\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}}$ is the standard deviation, such that $\sigma_{\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}}\to 0$ is consistent with a coplanar population and $\sigma_{\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}}\to\infty$ is consistent with an isotropic population. Finally, to convert between $\Delta i$ and $i_{P}$ (inclination of the planetary orbit measure with respect to the observer’s line of sight), we simply use the relationship $i_{P}=i_{\mathrm{bin}}+\Delta i$, where $i_{\mathrm{bin}}$ is obtained from an isotropic distribution as discussed above.

For $\Omega_{P}$, it is again crucial to realize that it is the relative longitude of ascending node that characterises coplanarity. We denote the relative longitude of ascending node as $\Delta\Omega$ (where $\Delta\Omega=\Omega_{P}-\Omega_{\mathrm{bin}}$) and elect to also distribute the longitude of the ascending node of the planet following another truncated Gaussian. We limit ourselves to consider only prograde planets (planets orbiting in the same direction as the binary stars), which means we want to sample $\Delta\Omega$ between $-\frac{\pi}{2}$ and $\frac{\pi}{2}$ instead of the complete $-\pi$ to $\pi$ range. To do so, we sample from a truncated Gaussian such that $\Delta\Omega/(\frac{\pi}{2})\sim\bar{\mathcal{N}}[0,\sigma_{\Delta\Omega}]$ (with support from $-1$ to $+1$ and centered at $0$), where now by virtue of normalizing by $\frac{\pi}{2}$ we maintain the same truncation limits as with $\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}$. Finally, we can convert between $\Delta\Omega$ and $\Omega_{P}$ using the relationship $\Omega_{P}=\Omega_{\mathrm{bin}}+\Delta\Omega$, where $\Omega_{\mathrm{bin}}$ is obtained from an isotropic distribution as discussed above.

These two parameters, $\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}$ and $\Delta\Omega$, both describe the degree of non-coplanarity, but along different dimensions. For example, a coplanar population would correspond to $\sigma_{\\ cos\big{(}\frac{\pi}{2}-\Delta i\big{)}}\to 0$ and $\sigma_{\Delta\Omega}\to 0$. In contrast, an isotropic population corresponds to $\sigma_{\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}}\to\infty$ and $\sigma_{\Delta\Omega}\to\infty$. In this way, and because of the carefully chosen parameterisation that maintains the same truncation limits, the degree of coplanarity can be encapsulated with just a single term, $\sigma\equiv\sigma_{\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}}=\sigma_{\Delta\Omega}$. To provide some context, the $\sigma$ of all observed transiting circumbinary planets is $\sim 5\times 10^{-4}$ (noting that $\sigma$ is in fact a dimensionless quantity).

In Figure 2, we consider three representative values of $\sigma$: $\sigma=10$ (corresponding to an isotropic population), $\sigma=0.3$ (corresponding to an intermediate population between coplanar and isotropic), and $\sigma=10^{-2}$ (corresponding to a coplanar population). For each value of $\sigma$, we show the corresponding probability distributions for $\cos\big{(}\frac{\pi}{2}-\Delta i\big{)}$ and $\Delta\Omega$. One can see that these, as well as the distribution of normal vectors on the globes, are consistent with the expected behaviours in limiting cases.

Now that all the relevant orbital parameters are established, we make the additional note that these parameters only set the initial configuration of a circumbinary planet system. They are by no means constant and can evolve over time according to the dynamics of the system. This time-dependence of the orbital parameters is fully captured by our simulations.

We are now ready to generate a random circumbinary system. Because we are investigating the effect of changing $\sigma$, these systems are of course not representative of the underlying distribution, whose $\sigma$ is not well-characterised. But the stellar masses and binary separations are designed to be quasi-representative (especially once we filter out unphysical systems in the next subsection).

We use the rebound package (Rein & Liu, 2012) to evolve the circumbinary system over a time span much longer than one Kepler observing window and record the number of transits over each epoch. Note that we do not distinguish between the two stars when counting transits. In practice, it is often possible to distinguish the identity of the star being transited - especially for $q\ll 1$ since the depths become increasingly distinct (e.g. Doyle et al. 2011). However, the misaligned nature of our systems means that the impact parameters are also constantly changing, leading to V-shaped transits and variable depths even for U-shaped events due to limb darkening (Mandel & Agol, 2002). This complicates the process of flagging identities to each event and introduces the possibility of mis-labelled events. In contrast, we argue that the simple act of counting how many dips are present in the data is far less ambiguous and thus presents a more robust labelling system, at the expense of some compromise in information content. Of course, even here there is the potential for some degree of confusion due to overlapping transits, but broadly one should expect the number of dips to be interpretable in such cases (e.g. see Gillon et al. 2017).

After generating and evolving each of the circumbinary systems, we are left with a sequence of integers that denotes the number of transits per epoch for each planet. Next, we impose several filters to ensure that we are left with systems that are both physically sensible and observable.

Among all the observed circumbinary planets, the shortest period is ${\simeq}50$ days (Kepler-47b; Orosz et al. 2012). Accordingly, we impose this as a constraint in our artificial circumbinary population and filter out all planets with periods shorter than 50 days.

We can only extract information from circumbinary planets that are observable to us. In what follows, an “observable” planet is defined to be one that undergoes at least two epochs with at least one transit in one Kepler observing window. Thus, we need to filter out all unobservable planets.

To determine if a planet is “observable” or not, we first calculate the number of epochs it undertakes during one Kepler window and denote it as $n$. Then, we randomly extract a sub-sequence of length $n$ from the full transits per epoch sequence we obtained from our simulation. If this sub-sequence does not contain at least two non-zero entries, them we discard it because it would outwardly lack the two necessary to cross our threshold of being “observable”. We repeat this procedure for all possible sub-sequences with length $n$ and filter out all the unobservable sub-sequences.

For an isotropic circumbinary planet population, $\sim 8$% of the systems passes through the filters mentioned above and contain at least one length $n$ sub-sequence that is “observable”. On the other hand, only $\sim 0.4$% of the systems in a coplanar planet population contain observable sub-sequences. The fact that coplanar systems are less likely to be observable seems counter-intuitive, but we note that for a coplanar system to be observable, the inclination of its binary and planet orbit must be precisely aligned to the observer’s line of sight. On the other hand, since the inclination of the binary and planet orbit is not correlated in the isotropic case, the requirement for observable systems is less stringent.

Recall that the underlying question tackled in this section is - do the transits per epoch sequences of transiting misaligned circumbinary planets follow a pattern which is discernible from being purely random? The hypothesis of a stochastic sequence could be undermined if were possible to make predictions for the next number in a sequence. Even if the accuracy of said predictions are not perfect, the demonstrable ability to predict better than random would still establish that the sequence was not truly stochastic in nature. To this end, we introduce a simple algorithm for predicting the number of transits a circumbinary planet undergoes in a given epoch using information from the previous epochs (hence the name "memory slot“), and show that the predictive power of this algorithm is consistently better than random guesses, thus proving that the transit occurrences of circumbinary planets are non-random and tractable.

With the predictive model defined, we next turn to a point of comparison. We consider two possible "benchmark" models to this end, which we introduce in what follows. Further, we introduce the important concept of "broken sequences".

It is straightforward to apply the memory slot prediction algorithm described in Section 2.4 to the validation set and measure the corresponding prediction accuracy. However, this prediction accuracy in isolation means very little since it remains unclear whether our model establishes non-randomness to the sequences. To better understand and contextualise this prediction accuracy, it is useful to compare it to a benchmark accuracy, derived using a "naive" algorithm (following Sandford et al. 2021).

We considered two possible benchmarks. The first is obtained using a simplistic algorithm: predict the mode of the dataset every time. In this way, the model simply predicts the same guess every time once “trained”. As a second model, we try taking the observed sequence and simply drawing a random choice from it as the “prediction” for the next epoch. We dub these benchmark models as the “modal” and “random” models. Our prediction algorithm can claim to have predictive power only if it can out-perform these benchmarks.

In considering the types of sequences that might be observed, one concern is that many sequences might be discarded by automated surveys as seemingly impossible. This is particularly poignant for misaligned circumbinary systems, since the binary star may not be eclipsing and thus its binary nature may not even be known beforehand. A sequence such as “1111” would appear clearly reasonable, even with large timing variations, as a likely single planet transiting a single star with some additional hidden perturber (e.g. see Nesvorný et al. 2012). A sequence such as “1121” would also be quite reasonable, with the “2” event possibly caused by a second longer period transiting planet or even an exomoon (e.g. Teachey & Kipping 2018) - which we note need not transit each epoch (Martin et al., 2019). However, a sequence such as “1101” is far more challenging to understand and apparently without precedent (although the lack of precedence may simply be a product of the fact that such sequences are unexpected and thus possibly rejected). With the “1121” case, we can cluster the transits within each epoch together and still draw a quasi linear ephemeris through the times of the clustered events, consistent with the exomoon hypothesis for example. But a “1101” sequence produces a sparse chain that cannot be represented by even a quasi linear ephemeris. We consider here that such sequences, which we dub “broken sequences”, may trip over many automated searches and could plausibly be absent in future catalogs.

We define “broken sequences” as those for which the time gap between observable epochs (i.e. non-zero numbers of transits) is not constant (and thus an “unbroken sequence” is the opposite case).

For example, a transits per epoch sequence of “1011” or “100101” corresponds to a broken system, and a sequence of “1211” or “10201” corresponds to an unbroken system. The fraction of simulations that generate broken sequences is thus an important parameter when considering the occurrence rate and completeness estimates, in cases where this feature would indeed negatively impact the detectability of such sequences in surveys.

We find that the ratio of broken to unbroken sequences is strongly affected by the coplanarity of the circumbinary population generated. Varying the governing $\sigma$ parameter in Figure 4, we find a log-logistic type behaviour emerges. This reveals that for quasi-isotropic populations, some 30% of sequences may present broken sequences in Kepler-like observing windows. Accordingly, surveys designed to measure the occurrence rate and properties of transiting circumbinary planets need to either use algorithms capable of identifying broken sequences or correct their rates for this effect appropriately.

As noted, the fraction of broken systems versus $\log(\sigma)$ curve in Figure 4 resembles a logistic function. This can be understood by the fact that at extreme values of $\sigma$, the distribution saturates to either perfectly isotropic or perfectly coplanar.

Broken systems could potentially cause a survey to reject a system as unphysical, although this of course depends upon the search algorithms being utilised. A sequence such as “1101” could be plausibly reconstructed though and the correct period inferred. However, the broken sequence “1010001” would be challenging to assign the true period for without detailed modelling or additional information; in particular we only observe the 1’s and so thus the simplest assumption is that the underlying sequence is “1101”. In such a case, we would have assigned an orbital period which is twice as long as the true one - an aliased period. Of course, this aliasing issue can also arise on unbroken sequences, such as the sequence “10101”. Thus, when interpreting such systems, besides from dividing the sequences up between broken and unbroken, we can also delineate them between “aliased” and “unaliased”.

This aliasing is important to our predictive model because it effectively transforms sequences like “1010001” $\to$ “1101”, which will clearly influence the training and validation of our proposed algorithm. However, if the true period were known, perhaps through detailed modeling or additional observations such as radial velocity measurements, then this transformation would not occur and we could use our algorithm without modification.

As with broken sequences, the ratio of aliased and unaliased sequences exhibits a log-logistic type behaviour as we vary the governing $\sigma$ parameter in Figure 4. For quasi-isotropic populations, more than a quarter of the sequences are aliased. This not only affects the training and validation of our proposed algorithm, but also implies that extra caution needs to be applied when assigning a period to observed transiting circumbinary planets.

To summarise, we have proposed a simple predictive model for the number of transits per epoch based on the previous two epochs. In addition, we have proposed two benchmark predictive models that simply exploit the ensemble statistics rather than the last two epochs, one that always returns the mode and one that chooses a random sample from the list. In the results that follow, we always compare the performance of our full predictive model to these simpler two as a benchmark.

Further, we note that observed sequences can be delineated in two different ways. First, a sequence is “broken” if a quasi linear ephemeris cannot explained even the clustered sequence of transits (otherwise its unbroken). Such cases could potentially be missed by detection algorithms and thus deserve special consideration. Second, a sequence could lead one to assign an aliased period for the planet rather than the true one. In this way systems can be broken or unbroken, and also aliased or unaliased.

In evaluating the performance of our predictive model against the benchmark models, there are thus a total of 4 different scenarios that could be considered

• A]

  Broken+unbroken sequences, where the true period is correctly assigned.

• B]

  Broken+unbroken sequences, where the apparent period is assigned (which may be aliased).

• C]

  Unbroken sequences only, where the true period is correctly assigned.

• D]

  Unbroken sequences only, where the apparent period is assigned (which may be aliased).

For each of these four scenarios, we start from the transits per epoch sequences obtained from generating and evolving populations of circumbinary planets as prescribed in Section 2.2. Then we filter out the broken sequences and transform the aliased sequences to reflect the apparent period when necessary. The code we used to treat each of the four scenarios is made publicly available at this URL.

For each of the four scenarios, we compare the performance of our predictive model with the two benchmark models described in Section 2.5. Our results are summarised in Tables 1,  2,  3, and  4.

The results from this exercise reveal several important insights. First, the “random” version of the benchmark algorithm performs worse than the “modal” version for isotropic distributions, but the accuracies of the two benchmarks converge as we approach the coplanar distributions. For isotropic cases, the large of number of zeros clearly favours the modal method but as the sequences become more varied, the two methods converge in accuracy.

However, as a second key observation, we note that our two memory slot predictions almost always out-perform either benchmark. The only exception is in scenario D for the highly isotropic cases, where the modal case wins out. We note that the number of one’s is exceptionally high in this regime, reaching 94%, and thus the modal algorithm achieves excellent accuracy purely as a result of the monotonous nature of the sequences. For near coplanar cases, the two memory slot method reaches an accuracy of ${\sim}92$% where as the benchmark methods cannot do better than ${\sim}50$%.

As a result of this second observation, we can arrive at a conclusion concerning the original motivation behind this section. The number of transits per epoch for misaligned (or aligned) circumbinary transiting planets does not follow a purely random sequence and can be effectively predicted. As before, we emphasise that our prediction algorithm is not optimal and encourage work to develop improved iterations.

As before, we generate and evolve populations of circumbinary planets as described in Section 2.2 and apply the filters described in Section 2.3. Furthermore, we do not assume prior knowledge of the planetary period, and thus sequences filled with many zeros (that can be explained by longer, aliased periods) are transformed with our aliasing transformation prescription (see Section 2.7).

The results are shown in Table 1 - 4 for each of the four scenarios described in Section 2.8. Scenarios B and D are particularly important, since they correspond to the case where one is simply tallying transits/epoch and is not equipped with special knowledge of the true period (which we consider to be the typical case in an automated survey). Scenarios B and D recall are distinct as to whether broken systems are included (B, yes; D, no) which is dependent upon the details of the search algorithm in question.

The results reveal how coplanar systems are indeed dominated by epochs with transits, as expected, but highly misaligned systems are far more sporadic. For isotropic distributions, we find that of the non-zero events, over $90$% produce single transits. With such a high fraction of single-transit epochs, there is clearly a danger in such systems being misclassified as planets orbiting single stars - especially for scenario D where broken systems are discarded in the first place. Timing variations would of course also be present here but those could be interpretable as planet-planet interactions. In addition, such systems will exhibit depth variations for $q<1$. This suggests that it may be useful to survey the existing catalog of planets ostensibly orbiting single stars to see if their depth and timing variations might also be consistent with a misaligned circumbinary scenario.

In practice, one might conceive that inference of $\sigma$ from a population of circumbinary planets would proceed through a full consideration of the transit shapes, durations, timing, etc to derive individual system posteriors, and then these would be fed into a population inference framework, such as hierarchical Bayesian modelling (e.g. see Hogg et al. 2010; Chen & Kipping 2017). Indeed, population level analyses of the alignment of circumbinary transits have already been published, such as that of Li et al. (2016) who favour a nearly coplanar population. However, given the poor sensitivity to misaligned systems with current methods (Martin, 2018), the influence of selection bias may be driving this.

In what follows, we seek a simplified description for inferring $\sigma$ that depends on the simple and well-defined metric of transit counts per epoch. Whilst surely more precise inferences of the inclination distribution could be inferred by leveraging the finer grained information encoded in the light curve, this metric is attractive for its simplicity, enabling a relatively straight-forward inference procedure.

To accomplish this, we first need an interpolating function of the grid of results presented earlier. In what follows, we use the results from Scenario D, as the likely typical scenario, but the below can be repeated for any of the other scenarios following the same prescription. Note that in this scenario, broken sequences are removed and the remaining unbroken sequences undergo an aliasing transform. This means that there are no zeros in the final sequences.

As was found in Figure 4, we find that logistic functions work well in this problem when the independent variable is $\log\sigma$. This is because at extreme values of $\sigma$, the distribution saturates to either perfectly isotropic or perfectly coplanar. Thus, we regressed logistic function fits to the 4 non-zero columns in the middle panel of Table 4. Additionally, we impose the constraint that probabilities must be greater than 0, less than 1, and sum to 1. To impose the “sum to unity” constraint, we choose to fit columns 3, 4, 5 of the middle panel first under the constraint that they sum to less than unity, then calculate the value of column 2 by hand²²2We also tried other combinations of these but found this to provide the best match to our simulation results..

The resulting fits are shown in Figure 5, where one can see the close agreement achieved between the logistic functions and the simulation results. Note that even for column 2 where we manually computed the values instead of directly optimising the fit, the errors are still quite small, which validates our fitting strategy.

The logistic function fits are as follows, where $f_{j}(\sigma)$ is the fraction of epochs exhibiting $j$ transits:

and by definition

Naturally, the above could be easily repeated for the other scenarios, where we note that one would need to include an $f_{0}$ function to account for zero-transit epochs.

Consider that we have observed $N$ transit epochs of $M$ distinct transiting circumbinary planets. The $M$ objects form a population, for example a group of known circumbinary systems with similar stellar parameters. Let us assume that each member of the population has a unique $\Omega_{\mathrm{bin}}$ and $i_{\mathrm{bin}}$ value, and that that the ensemble of these values can be approximately described by the truncated normal distributions posited in Section 2.2 and depicted in Figure 2. In this way, each object has a unique orbital alignment but the population is characterised by a single parameter - $\sigma$.

From the sample of $N$ epochs, we can count the number of cases where just a single transit is observed, $n_{1}$, or two transits, $n_{2}$, etc. Defining $n_{4}$ as the number of of 4+ transits, we expect $N=n_{0}+n_{1}+n_{2}+n_{3}+n_{4}$. Note that, in general, $n_{0}$ is not directly observed, but rather inferred from the lack of a transit and an inferred ephemeris. However, in our case of assuming scenario D, it is not possible to have a zero-transit epoch as discussed earlier so we can discard and define $N=\sum_{k=1}^{4}n_{k}$. Equipped with the data $\mathcal{D}=\{n_{1},n_{2},n_{3},n_{4}\}$, we might ask - what is the allowed range for $\sigma$ that can accommodate these observations?

To make this inference, it is necessary to first define a likelihood function, $\mathcal{L}$. The likelihood of obtaining the data $\mathcal{D}$, given a choice of $\sigma$, can be expressed as a multinomial probability mass function, since the problem is discretised:

Since the binomial term in the above is simply a constant for a given set of data, $\mathcal{D}$, the objective function for inference can be expressed as

The likelihood function is typically a monotonic function of $\sigma$ and thus a central value of $\sigma$ is not well-defined (see Figure 6). Instead, it is typically more useful to state that $\sigma$ tends toward the coplanar/isotropic asymptote, and then assign a statistical upper/lower limit on $\sigma$ parameter, respectively.

To show an example, we set $\sigma=0.01$ and sampled $N=100$ epochs to produce the fake data set $\mathcal{D}=\{n_{1}=52,n_{2}=39,n_{3}=5,n_{4}=4\}$. Once again, we emphasise that because we are working in scenario D, there are no zero transit epochs by construction. The likelihood function is shown in the top panel of Figure 6, plateaus at small $\sigma$ and suggests a nearly coplanar solution. Credible intervals for 68.3% and 95.5% are shaded in hatch and cross-hatch respectively, which are found to be consistent with the injected truth. The exercise was repeated for $\sigma=10$ and produced another self-consistent recovery (see lower panel of Figure 6), despite the fact the likelihood function flips over and favours an isotropic solution here.

The inference model here is simplistic and could surely be improved by leveraging other information within each epoch (e.g. transit durations), labelling the transits by star and considering alternative population models beyond a truncated normal. But our simple demonstration establishes that transit sequences already provide constraining power on the population’s inclination distribution and the problem is tractable.