format: text

text:
# Survival and Extinction in Cyclic and Neutral Three–Species Systems

## 1 Introduction

Cyclic phenomena are often ignored when studying epidemiological and evolutionary processes, but nevertheless they can have important consequences, e.g. we contract most infective diseases only once in our lifetime, because our immune system has “memory”. Vaccines are designed based on this knowledge, and they work quite efficiently. However, it is widely accepted that some viruses, such as the flu virus, mutate in order to make themselves unrecognizable by our immune system, and thus be able to reinfect us. Immunity to pertussis (whooping cough) is temporary, and decreases as the time after the most recent pertussis infection increases. The chicken pox (agent VZV—varicella zoster virus) also repeats, because immunity wanes with time. Our model is motivated by this type of situation.

In epidemiological models populations are often categorized in three states: susceptible $`(S)`$, infected $`(I)`$ and recovered $`(R)`$. There exists a vast literature about the so-called $`SIR`$ models, where the “loss of immunity” step is not considered (e.g. the classic texts by Bailey , Anderson and May , the review by Hethcote ). We consider the case when the “loss of immunity” step is present, otherwise referred to as $`SIRS`$ models .

Recently it has been shown that many species of bacteria are able to produce toxic substances that are effective against bacteria of the same species which do not produce a resistance factor against the toxin . These bacteria exist in colonies of three possible types: sensitive $`(S)`$, killer $`(K)`$ and resistant $`(R)`$. The killer type produces the toxin, and the resistance factor that protects it from its own toxin, at a metabolism cost. The resistant strain only produces the resistance factor (at a smaller metabolism cost). The sensitive strain produces neither. Clearly, the $`S`$ colony can be invaded by a $`K`$ colony, the $`K`$ can be invaded by $`R`$, and the $`R`$ can be invaded by $`S`$.

In our first model the population is categorized into three “species”: $`A`$, $`B`$, $`C`$, and the rules are such that when an $`A`$ meets a $`B`$, it becomes $`B`$, when $`B`$ meets $`C`$, it becomes $`C`$, when $`C`$ meets $`A`$, it becomes $`A`$. The total population size is conserved. Otherwise this could be seen as a three-party voter model, when the follower of a certain party “converts” a follower of another certain party when they meet.

The reaction is similar to the cyclical “Rock-paper-scissors game”, of which a biological example has been found recently . It is played by the males of a lizard species that exist in three versions: the blue throat male defends a territory that contains one female, the orange throat male defends a territory with many females, and the male with a throat with yellow stripes does not defend his own territory, but can sneak into the territory of the orange males and mate with their females. Hence, a “blue” population can be invaded by “orange” males, while an “orange” population is vulnerable to “yellow” males, who on their turn are at an disadvantage against the “blue” males, who defend their territory very well. Many spatial models that relate to such systems have been built .

The second model is a three–component version of the famous Kimura model of neutral genetic drift . In that case, when individual from two species meet, the offspring may be either of the first or the second species, with equal probability.

## 2 Description of the Model

Consider a system in which three species $`A`$, $`B`$, $`C`$ are competing in a way described by the reaction: $`A+B2B`$, $`B+C2C`$, $`C+A2A`$.

The rate equations for this system will be:

$`N{\displaystyle \frac{dA}{dt}}=ACAB`$
$`N{\displaystyle \frac{dB}{dt}}=BABC`$ (1)
$`N{\displaystyle \frac{dC}{dt}}=CBCA`$
with $`A+B+C=N=`$ const. (we assume the rates are the same, in which case a time–rescale will remove them from the equations.) These equations can be rewritten as:

$`N{\displaystyle \frac{d}{dt}}\mathrm{ln}A=CB`$
$`N{\displaystyle \frac{d}{dt}}\mathrm{ln}B=AC`$ (2)
$`N{\displaystyle \frac{d}{dt}}\mathrm{ln}C=BA`$
which leads to the second conservation rule: $`ABC=H`$ = const.

The above model contrasts with the neutral drift model ($`A+B2A`$ or $`2B`$, $`B+C2B`$ or $`2C`$, $`C+A2C`$ or $`2A`$) for which:

$$\frac{dA}{dt}=\frac{dB}{dt}=\frac{dC}{dt}=0$$
(3)
Assuming that the system is subject to stochastic noise due to Poisson birth and death processes (intrinsic noise) we get the master equation for the cyclic model:

$$\frac{P(A,B,C,t)}{t}=\frac{1}{N}[(A1)(C+1)P(A1,B,C+1,t)ACP(A,B,C,t)+$$

$$+(A+1)(B1)P(A+1,B1,C,t)ABP(A,B,C,t)+$$

$$+(B+1)(C1)P(A,B+1,C1,t)BCP(A,B,C,t)]$$
(4)
For the neutral drift case the master equation is:

$$\frac{P(A,B,C,t)}{t}=\frac{1}{N}[\frac{1}{2}(A1)(C+1)P(A1,B,C+1,t)+$$

$$+\frac{1}{2}(A+1)(C1)P(A+1,B,C1,t)ACP(A,B,C,t)+$$

$$+\frac{1}{2}(A+1)(B1)P(A+1,B1,C,t)+\frac{1}{2}(A1)(B+1)P(A1,B+1,C,t)$$

$$ABP(A,B,C,t)+\frac{1}{2}(B+1)(C1)P(A,B+1,C1,t)+$$

$$+\frac{1}{2}(B1)(C+1)P(A,B1,C+1,t)BCP(A,B,C,t)]$$
(5)
Now we introduce the “shift” operators, defined by

$$ϵ_Af(A,B,C)=f(A+1,B,C)$$
(6)
and

$$ϵ_A^1f(A,B,C)=f(A1,B,C)$$
(7)
and likewise for B and C operators. The master equation for the cyclic model now reads:

$$\frac{P(A,B,C,t)}{t}=\frac{1}{N}[(ϵ_Cϵ_A^11)AC+(ϵ_Aϵ_B^11)AB+(ϵ_Bϵ_C^11)BC]P(A,B,C,t)$$
(8)
and for the neutral drift one:

$$\frac{P(A,B,C,t)}{t}=\frac{1}{N}[(\frac{1}{2}(ϵ_Cϵ_A^1+ϵ_Aϵ_C^1)1)AC+(\frac{1}{2}(ϵ_Aϵ_B^1+ϵ_Bϵ_A^1)1)AB+$$

$$+(\frac{1}{2}(ϵ_Bϵ_C^1+ϵ_Cϵ_B^1)1)BC]P(A,B,C,t)$$
(9)
Next we transform to the intensive quantities

$$x=\frac{A}{N},y=\frac{B}{N},z=\frac{C}{N}$$
(10)
We use the system size expansion of Horsthemke and Brenig . However, the notation and style is closer to Van Kampen .

The shift operators become

$`ϵ_{A(B,C)}={\displaystyle \underset{j}{}}{\displaystyle \frac{1}{j!}}N^j{\displaystyle \frac{^j}{x^j(y^j,z^j)}}`$ (11)
$`ϵ_{A(B,C)}^1={\displaystyle \underset{j}{}}{\displaystyle \frac{1^j}{j!}}N^j{\displaystyle \frac{^j}{x^j(y^j,z^j)}}`$
Further we use the rules of transformation of random variables to define:

$$W(x,y,z,t)=P(A,B,C,t)N^3$$
(12)
Before we go ahead with the expansion, a few comments are necessary. From the rate equations, it is clear that our system does not have one single steady state or limit cycle. Instead, the mean-field solution is completely dependent on the initial conditions, and we do not have a macroscopic solution to expand about, but rather an infinity of neutrally stable cycles (for the competition case) or points (for the neutral model.) The above expansion, otherwise known as Kramers-Moyal expansion, is discussed by van Kampen . It is a risky expansion, and it does not work in most cases, mainly because higher derivatives are not small themselves. Also, with time, large fluctuations may occur, and since our system presents itself as a fluctuations–driven one, it becomes even more important to watch out for this sort of complications. In the case of neutral drift, the second order term is the leading one, so we would expect this Ansatz to work for large system sizes. In the cyclic competition model, the first order term does not drive the system out of the (neutrally) stable trajectories, and again it would be the second order term to essentially determine the fate of the system. So, we set out to investigate the applicability of the Kramers-Moyal expansion to the cyclic competition and neutral drift three-species systems.

We transform the master equation into an equation for $`W(x,y,z,t)`$ by substituting the expressions for the shift operators in it, and further by grouping together the terms of the same order in $`N`$. The terms of order $`N^1`$ cancel, and the first terms in the expansion are of order $`N^0`$; keeping only those gives us the first order equation for the cyclic case:

$$\frac{W}{t}=[(\frac{}{x}\frac{}{z})xz+(\frac{}{y}\frac{}{x})yx+(\frac{}{z}\frac{}{y})zy]W$$
(13)
For the neutral drift case the first-order term is identically zero.

The second order term is obtained by considering the terms of order $`N^1`$ in the expansion of the master equation:

$$\frac{1}{2N}[(\frac{^2}{x^2}2\frac{^2}{xz}+\frac{^2}{z^2})xz+[(\frac{^2}{x^2}2\frac{^2}{xy}+\frac{^2}{y^2})xy+$$

$$+[(\frac{^2}{y^2}2\frac{^2}{yz}+\frac{^2}{z^2})yz]W$$
(14)
and is identical for both cyclic and neutral drift cases.

## 3 The First Order Term

The concentrations $`x`$, $`y`$, $`z`$ of a three–component system are commonly represented by the distances of a point inside an equilateral triangle of unit height from its sides (Fig. 1):

For the neutral drift case the first order equation tells us that in the mean–field approximation the system will remain in the initial state forever.

After some algebra, the first order equation  for the cyclic system becomes:

$$\frac{W}{t}=\frac{}{x}(xzxy)W+\frac{}{y}(xyyz)W+\frac{}{z}(yzxz)W$$
(15)
At the centre of the triangle ($`x=y=z=1/3`$) all three expressions above are zero, so in the first order approximation the system will stay at that state forever. It can be verified that the product $`H=xyz`$ is a solution of the first-order equation, as is any function of that product. The lines $`xyz=`$ const. will then represent possible trajectories of the system, which are closed, since there is no term to drive the system out of those trajectories. In this (mean-field) approximation, the cyclic model exhibits global stability with constant system size $`N`$ . The populations of each species oscillate out of phase with one-another, and the amplitude remains constant . These are the same trajectories that are obtained when one solves the rate equations . Fig. 2 shows some of those closed trajectories.

## 4 The Second Order Term

The second order term , which represents the “diffusion term” must be added to the first order equation to give the Fokker-Planck equation for our system. This term is identical for both the cyclic (competing species) and neutral genetic drift cases. As we will see, this makes the long (evolutionary) time fate of both systems essentially the same, at least for reasonably large system sizes.

The trajectories of a given realisation of the system can easily be obtained by simulation: they initially spiral out of the centre, with the populations of each species oscillating out of phase with one–another, so that the total size of the population is conserved. The amplitude drifts until one of the populations becomes zero (species goes extinct). In other words, the slow second order term makes the system cross from a closed (mean-field) trajectory to a neighbouring one. This would suggest an adiabatic approximation. In the neutral drift model, the fluctuations drive the system from one point–solution of the mean-field equations to a neighbouring one.

The computer simulations of both systems start with an equal number of $`A`$, $`B`$, $`C`$, (start at the centre) and generate times for the next possible reaction event with exponential distribution as $`N\mathrm{ln}(rn)/AB`$, (for $`A+B2B`$ reaction, and similarly for the other two reactions,) where $`rn`$ is a random variable with uniform distribution in $`[0,1]`$ (this means the events are really independent) . The reaction which occurs first is then picked and the system is updated. The process is repeated until one of the species, and then another one, go extinct. Fig. 3 shows the variation with time of the number of $`A`$ for a realization of the system in the cyclic competition scenario, and Fig. 4 shows the evolution of the number of $`A`$ population in the neutral drift case (in both realizations the total population size is $`N=600`$.) From the time series for the population numbers it looks as if the neutral drift model is some sort of “adiabatic” approximation of the cyclic competition model. If in the neutral drift model the population number is random walking, in the cyclic competition one the amplitude of oscillations is random walking! In other words, if we average the cyclic competition model over the cyclic orbits, we get essentially the same behaviour as the neutral drift model. Then the fluctuations drive both models from a neutrally stable cycle (point) to the neighbouring ones, and these fluctuations remain quite small by virtue of the vicinity of the neighbouring cycle (point.) It is quite funny how we all tend to go around in circles in life, even though the point we end up to is still the same!

We investigated the probability distribution of first extinction times. For that we ran fifty thousands copies of the system for each (different) size of the system, and the number of survivors was plotted vs. time. In the semilog axes we get a straight line, except for the few “rare events”. The slope of this straight line is -3.01 for system size 3000, -2.98 for system size 6000, -2.99 for system size 9000. Fig. 5 shows the plots for these three different system sizes: 3000, 6000, and 9000.

The time dependence of survival probability scales approximately with $`1/N`$, but there is still a very weak dependence left for small $`N`$. This can be justified by looking at the Fokker-Planck equation for this system, which contains both a term of order $`N^0`$ and $`N^1`$.

To further investigate the behaviour of the system, we looked at the cummulative, conditional on being alive (which is equivalent to a normalization condition,) probability distribution for $`H=xyz`$. If $`P(H<h|`$alive$`)`$ increases linearly with $`H`$ for small $`H`$, the extinction rate will be non-zero, and we will get exponential decay. This corresponds to a constant distribution for $`P(H|`$alive$`)`$. Fig. 6 shows these plots for $`t=1.25`$, $`t=1.5`$, $`t=1.75`$, and $`t=2.0`$ in linear axes.

For the neutral drift case the data from the simulations give us a slope of -2.992 for system size 6000, -2.991 for system size 3000. The data collapse when time is rescaled with $`1/N`$, and that is in agreement with the fact that there is only a term of the order $`1/N`$ present in the Fokker-Planck equation. In Fig. 7 we show the number of survivors vs. time plots for cyclic (competition) and drift cases together, for system size $`N=6000`$. They clearly agree.

## 5 Probability Distribution at Long Times

In order to get an expression for $`W`$ for long times, we looked at our “experimental” data. It is possible to find a complete set of eigenfunctions of the Laplace equation, which vanish at the boundary. Having this complete set of eigenfunctions for the equilateral triangle , the probability density can be written:

$$W=\underset{m,n}{}c_{m,n}\varphi _{m,n}$$
(16)
where $`\varphi _{m,n}`$ are the abovementioned symmetric eigenfunctions. The non-normalized eigenfunctions are then obtained from the general expression:

$$\varphi _{(m,n)}(x,y)=\underset{(m,n)}{}\pm \mathrm{exp}[\frac{2\pi i}{3}(nx+my)]$$
(17)
where we are using the relation $`x+y+z=1`$ and summation over the index pair $`(m,n)`$ means summation for $`(n,mn)`$, $`(n,m)`$, $`(nm,m)`$, $`(nm,n)`$, $`(m,n)`$, $`(m,mn)`$ .

We get a symmetric eigenfunction when $`m+n`$ is a multiple of 3, $`m`$ is also a multiple of 3, but $`m2n`$, and $`n2m`$. 111In , of all the symmetric eigenfunctions, Pinsky only considers the non-degenerate ones. There are also two-fold degenerate symmetric eigenfunctions, which we verified to be orthogonal to the original eigenfunctions.

In order to calculate the expansion coefficients, we took snapshots of the system at given times, and then used the relation:

$$c_{(m,n)}(t)\frac{1}{p}\underset{l=1}{\overset{p}{}}\varphi _{(m,n)}(x_l(t),y_l(t))$$
(18)
where $`p`$ is the number of experimental points. The time-evolution of the expansion coefficients is given in Fig. 8. The doubly–degenerate functions die out quite soon, while the Lamé symmetric functions persist.

With the coefficients obtained this way we can then construct the probability density at different times. Some snapshots at the $`W(t)`$ are shown in Fig. 9.

At $`t=0`$ the $`W(0)`$ is a $`\delta `$-function. As time passes, for $`t=0.1`$ we see the $`W`$ starts to spread, and it spreads even more for $`t=0.2`$, becoming a “cake” for very long times ($`t=1.5`$,) for both the cyclic system and the neutral drift one. The “cake” obtained this way shows some Gibbs oscillations, which are due to the inclusion of only a few terms in the expansion, and the presence of the absorbing boundary. For a uniform probability distribution the expansion coefficients are proportional to $`1/m`$, where $`m`$ is the first index in the pair $`(m,n)`$. It can be seen that our “experimental” expansion coefficients approach those values.

The second–order Fokker–Planck equation for the drift case accepts solutions of the form:

$$W(x,y,z,t)=e^{\lambda t}W(x,y,z)$$
(19)
with $`\lambda `$ the smallest eigenvalue of the spatial equation. We solved the eigenvalue problem by using the Galerkin method, with the help of Maple, and obtained a value of $`\lambda =3.01`$ for the smallest eigenvalue, when we keep the first six eigenfunctions in the expansion. This is in very good agreement with our “experimental” data.

## 6 Case When Rates Are Not the Same

Now consider the cyclic (competition) case when the rates of the three equations are not the same, i. e. the rate equations read like:

$`N{\displaystyle \frac{dA}{dt}}=c_{13}ACc_{12}AB`$
$`N{\displaystyle \frac{dB}{dt}}=c_{21}ABc_{23}BC`$ (20)
$`N{\displaystyle \frac{dC}{dt}}=c_{32}BCc_{31}AC`$
We can always rescale the $`A,B,C`$ so that the matrix $`c_{ij}`$ be symmetric, and the equations read:

$`N{\displaystyle \frac{dA}{dt}}=A(\beta C\gamma B)`$
$`N{\displaystyle \frac{dB}{dt}}=B(\gamma A\alpha C)`$ (21)
$`N{\displaystyle \frac{dC}{dt}}=C(\alpha B\beta A)`$
Again, $`A+B+C=N=`$const. By manipulating the rate equations, we can find the second integral of motion to be $`A^\alpha B^\beta C^\gamma =H=`$const. The new centre will then be not at the point (1/3, 1/3, 1/3), but at $`A=\alpha ,B=\beta ,C=\gamma `$. The lines along which the quantity $`A^\alpha B^\beta C^\gamma `$ is constant are shown in Fig. 10.

We need to check whether the exponential decay behaviour is universal when the symmetry is broken this way. For this we ran fifty thousand simulations for different combinations of $`\alpha `$, $`\beta `$, and $`\gamma `$, for system size $`N=3000`$, starting from the new centre ($`A=N\alpha `$, $`B=N\beta `$, $`C=N\gamma `$.) The number of survivors vs. time plots were again obtained, and those plots show that the exponential decay behaviour is robust. The time scale is now dependent on the initial distance from the boundary, with the equal rates case having the largest distance, and so the largest time scale (and the smallest slope.) Some of those plots are shown in Fig. 11. In that figure, the times for the equal rates case (for which the sum of rates is 3) are multiplied by three, to make them comparable to the times for the unequal rates case, for which the sum of rates is 1.

## 7 Conclusions

We have considered an $`ABC`$ model in both the cyclic competition and neutral genetic drift versions, and studied the long-term behaviour of such a model. The number of the $`A,B,C`$ species in the cyclic competition case oscillates with a drifting with time amplitude, until one of the species (and then the next one in the cycle) goes extinct. In the neutral drift case the number of the $`A,B,C`$ species drifts, until one of them becomes zero. In both scenarios the number of survivors vs. time plots show an exponential decay, with the same exponent. The result is verified by writing and solving the Fokker–Planck equation for the second model. Finally, its robustness is checked against variations in the rate of the three different reactions in the system. It is very interesting to note that there is no difference in the time scale for the ensemble of species in cyclic competition and the case of neutral genetic drift.

There is growing concern about the effects of habitat fragmentation in the survival of the species . If the population of a certain species goes extinct in one patch (e.g. a herd, school, swarm) while it still survives in other patches, then what is known as “rescue effect” can prevent global extinction . Otherwise, the species is doomed to go extinct altogether. The models with individuals of “species” $`A`$, $`B`$, and $`C`$ distributed in a lattice have been studied by Szábo et al. . If we try to model habitat fragmentation as an ensemble of patches, each of those patches can be considered as one copy of our non-spatial system. However, with continuous migration in and out of the patch, the non-spatial picture described in this paper would be seriously perturbed. One aspect of this immigration and emigration has been recently discussed by Togashi and Kaneko . In a first approximation, the broad range of extinction times could relate to the persistence of species that exist in different “versions” such as certain lizards, or even some kinds of bacteria. Applied to its epidemiological scenario, it would relate to the endemicity of certain diseases, namely those with mutating pathogen, and then the outbreaks of epidemics in certain patches (areas). The many patches will then form a network . We have work in progress which places our $`ABC`$ ensemble in one of these small worlds with a scale-free topology.

## 8 Acknowledgements

We thank Nicolaas G. van Kampen for constructive criticism, and Michael Döbeli for helpful and interesting discussions and suggestions.


text:
# Using All-Sky Surveys to Find Planetary Transits

## 1 INTRODUCTION

In the past three years, great strides have been made in the detection of extrasolar planets (XSPs). To date, nearly all of the roughly 100 known XSPs have been discovered using the radial velocity (RV) technique. However, RV detections, in and of themselves, yield only a few planetary parameters, namely the period $`P`$, the eccentricity $`e`$, and $`M\mathrm{sin}(i)`$, where $`M`$ is the mass of the planet and $`i`$ is the inclination of its orbit. By contrast, if a planet transits its host star, much more information is available. First, of course, the $`M\mathrm{sin}(i)`$ degeneracy can be broken. Second, the ratio of the radii of the planet and host star can be measured. Therefore, provided that the star can be classified well enough to determine its mass and radius, then the planet’s radius and hence its density can be determined. Third, and perhaps most important, if the transits can be observed with sufficient signal-to-noise ratio (S/N), then one can probe otherwise unobservable details of a planet, such as its oblateness (Hui & Seager, 2002), atmospheric conditions (Charbonneau et al., 2002), and perhaps satellites and rings. Regardless of how a planet is initially discovered, once it is determined to transit its host star, this wealth of information can in principle be extracted by intensive follow-up observation of these transits. This fact has been amply demonstrated by the discovery and analysis of the transiting planet HD209458b (Charbonneau et al., 2000; Cody & Sasselov, 2002).

At the moment, all ongoing and proposed transit surveys are carried out in relatively narrow pencil beams. They make up for their small angular area with relatively deep exposures. These surveys fall into two basic classes: field stars (Howell et al., 2000; Brown & Charbonneau, 1999; Mallen-Ornelas et al., 2001; Udalski et al., 2002)111 http://www.psi.edu/{̃}esquerdo/asp/asp.html
http://www.hao.ucar.edu/public/research/stare/stare.html
http://bulge.princeton.edu/{̃}ogle/, and clusters (Street et al. 2000; Burke et al. 2002)222 http://star-www.st-and.ac.uk/{̃}yt2/WEB\_GROUP/top.html. These surveys are potentially capable of establishing the frequency of planets in various environments, but they are unlikely to find the kinds of transits of bright stars that would be most useful for intensive follow-up analysis. Although some of the surveys of field stars are considered “wide field”, their total survey areas are small compared to $`4\pi `$ str. One project that has the potential to cover a very large area is WASP (Street et al., 2002), which plans to employ five cameras, each with a $`9.^{}5\times 9.^{}5`$ field of view.

An alternative method is to conduct an all-sky survey. Instead of continuous observation of all targets (which is impossible from a practical standpoint for an all-sky survey), this approach would necessarily involve revisiting each target in the sky at regular, semi-regular, or random intervals throughout the course of the project. This kind of observing strategy will not yield a continuous light curve on any star, as the current transit surveys do. Rather, this plan will generate an only sporadically sampled light curve. However, the long time baseline for the survey will eventually generate just as many individual observations of a single star. Transit-like dips in the data stream will not be visually obvious, but by repeatedly phase-folding the full light curve back on itself over a range of periods, one can detect the dips from the transits. (See § 3.)

This approach is especially relevant given the fact that there are several all-sky surveys already being planned for objectives other than transit detections. It should be possible, for instance, to utilize the photometric data stream of upcoming astrometric missions for transit detection. Space-based projects such as GAIA333http://astro.estec.esa.nl/GAIA/gaia.html and DIVA444http://www.ari.uni-heidelberg.de/diva/diva.html would take hundreds of observations of millions of stars over mission lengths of years with the aim of obtaining precise astrometry. These data could equally well be analyzed for planetary transits.

There are also several existing or proposed ground-based all-sky surveys, including the Large-aperture Synoptic Survey Telescope (LSST)555http://www.lssto.org/lssto/index.htm, the Panoramic Optical Imager (POI), and the All-Sky Automated Survey (ASAS)666http://www.astrouw.edu.pl/~gp/asas/asas.html (Pojmanski, 2000). These surveys will be imaging the entire sky every few days, a qualitatively similar cadence to those of GAIA and DIVA. While LSST and POI will likely be saturated by stars of $`V12`$, which are of greatest interest for transit follow-ups, ASAS is of particular interest in the present context because of its very small aperture.

In this paper, we examine the process of analyzing photometric data streams from all-sky surveys to find transits. We calculate the sensitivity to XSP detection, the distance out to which these detections will be possible, and the number of false-positive detections due to random noise. We derive a general expression for the number of stars that can be probed by this technique as a function of the total sensitivity of the survey. We apply our analysis to the problem of telescope design and conclude that very small, 5 cm telescopes are optimal for finding transits of bright stars.

## 2 SCALING RELATIONS

In order to accurately determine the number of transits a given survey would be expected to detect, we must carefully define the set of stellar systems that can be probed for transits by the survey. This will evidently depend on the number density $`n`$, the luminosity $`L`$, and the radius $`R`$, of the stars being probed; as well as the semi-major axis $`a`$ and radius $`r`$ of their planets. A quick (and naive) formulation would then state that the total number of systems that can be probed for transits $`N_p`$ is

$$N_p=\frac{4}{3}\pi n\frac{R}{a}[d_{max}(L,R,a,r)]^3$$
(1)
for a homogeneous population of such stars, where $`d_{max}`$ is the distance out to which a transit can be detected.

However, the quantity $`d_{max}`$ is ill-defined. For fixed $`L`$, $`R`$, $`a`$, $`r`$ and distance $`d`$, the transit detection scales as $`(1x^2)^{1/2}`$ where $`x`$ is the transit impact parameter normalized to $`R`$; that is, $`0x1`$. Hence, the signal-to-noise ratio scales as $`(1x^2)^{1/4}`$. Since detections normally require a minimum S/N, $`d_{max}`$ must also be a function of the impact parameter $`x`$. Thus, for a photon-noise limited survey, $`d_{max}(1x^2)^{1/4}`$, and therefore

$$\frac{dN_p}{dx}=\frac{4}{3}\pi n\frac{R}{a}[d_{max}(L,R,a,r,x)]^3$$
(2)
where $`d_{max}(L,R,a,r,x)=d_{max}(L,R,a,r,0)(1x^2)^{1/4}`$. Here $`d_{max}(L,R,a,r,0)`$ is the distance out to which a transit can be detected for an edge-on ($`i`$ = $`90^\mathrm{o}`$) orbit. We must then integrate over all values of the impact parameter $`x`$ from 0 to 1, and so

$$N_p=\frac{dN_p}{dx}=\frac{4\pi nR}{3a}\eta [d_{max}(L,R,a,r,0)]^3$$
(3)
where,

$$\eta =_0^1(1x^2)^{3/4}𝑑x=\frac{\sqrt{\pi }}{2}\frac{\mathrm{\Gamma }(7/4)}{\mathrm{\Gamma }(9/4)}\mathrm{\hspace{0.17em}0.719}.$$
(4)
We now determine the dependence of $`N_p`$, the total number of systems probed, on the remaining parameters $`L`$, $`R`$, $`a`$, and $`r`$. To do so, we analyze the detection requirement,

$$N_t\left(\frac{\delta }{\sigma }\right)^2\mathrm{\Delta }\chi _{\mathrm{min}}^2,$$
(5)
where $`N_t`$ is the number of observations of the transit over the length of the survey, $`\delta `$ is the fractional change in the star’s brightness during the transit, $`\sigma `$ is the fractional error of an individual flux measurement, and $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ is the minimum acceptable difference in $`\chi ^2`$ between a fit that assumes a constant flux and one that takes account of a transit. As we discuss in § 3, $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ must be set sufficiently high to avoid spurious detections due to random noise.

* To determine the dependence of $`N_p`$ on $`L`$, we note that in equation (3), the only factor that depends on $`L`$ is $`d_{max}^3`$. For a particular star, the flux $`f=L/(4\pi d^2)`$, thus $`dL^{1/2}`$, and so $`N_pL^{3/2}`$.
* For the dependence of $`N_p`$ on $`R`$, we note that $`fd^2`$, and that $`\sigma f^{1/2}`$, so $`d\sigma `$. From equation (5), we see that $`\sigma \delta N_t^{1/2}`$. Since $`N_t`$ is the total number of observations of the star during transits over the length of the survey, $`N_t=N_{\mathrm{obs}}(2R)/(2\pi a)`$. Also, $`\delta =(\pi r)^2/(\pi R)^2`$. Thus, $`d_{max}R^{3/2}`$. Combining these relations with the explicit factor of $`R`$ in equation (3) itself, we arrive at $`N_pR^{7/2}`$.
* For the dependence of $`N_p`$ on $`a`$, we see that $`d\sigma `$, and $`\sigma \delta N_t^{1/2}`$. Using $`N_ta^1`$, we have $`d_{max}^3a^{3/2}`$, and so $`N_pa^{5/2}`$.
* Finally, for the dependence of $`N_p`$ on $`r`$, the only factor that depends on $`r`$ is $`\delta `$, from $`\delta =(\pi r)^2/(\pi R)^2`$. Thus $`\sigma r^2`$, and so $`N_pr^6`$.

Consolidating the dependence of $`N_p`$ on the various parameters, we finally arrive at,

$$N_p=\frac{4}{3}\pi n\eta d_0^3\left(\frac{R_0}{a_0}\right)\left(\frac{L}{L_0}\right)^{3/2}\left(\frac{R}{R_0}\right)^{7/2}\left(\frac{a}{a_0}\right)^{5/2}\left(\frac{r}{r_0}\right)^6,$$
(6)
where $`d_0`$ = $`d_{max}(L_0,R_0,a_0,r_0,0)`$. That is, $`d_0`$ is the distance out to which a planet-star system with $`i=90^{}`$ and the fiducial parameters $`L_0,R_0,a_0,r_0`$ can just barely be detected at the S/N threshold.

We now seek to simplify equation (6) by integrating over the local stellar population at fixed absolute magnitude, $`M_V`$, and so we replace the three independent variables $`(n,R,L)`$ by the single variable $`M_V`$. We consider two regimes $`M_V6`$ and $`M_V6`$ (with one overlapping bin at $`M_V=6`$, which we will use later to check for consistency).

We first treat the fainter regime. Here, the main sequence is relatively narrow. Hence, $`R`$ may be regarded as a function of $`L`$ (and so, therefore, of $`M_V`$), while $`n`$ is simply the number density of stars in a given magnitude bin. Hence, the “integration” amounts to a simple multiplication of factors.

We adapt the number density of stars $`n`$ from the empirically determined local stellar luminosity function (LF): for the range (9 $``$ $`M_V`$ $``$ 18) we use the LF reported in Zheng et al. (2001), and for the range (6 $``$ $`M_V`$ $``$ 8) we use the LF of Bessell & Stringfellow (1993).

To estimate stellar radii, we combine the linear color-magnitude relation, $`M_V=3.37(VI)+2.89`$ from Reid (1991), a color/surface-brightness relation $`\mathrm{log}(R/R_{})=0.69+\mathrm{\hspace{0.17em}0.2226}(VK)0.2M_V`$, based on the data of van Belle (1999), and $`VIK`$ color-color relations for dwarfs from Bessell & Brett (1988). We can therefore calculate the relative number of systems with a fixed $`a`$ and $`r`$ as a function of $`M_V`$, which we designate

$$F(M_V)=\left[\frac{n(M_V)}{n_0}\right]\left[\frac{L(M_V)}{L_0}\right]^{3/2}\left[\frac{R(M_V)}{R_0}\right]^{7/2}$$
(7)
where $`n_0,L_0,`$ and $`R_0`$ are the normalizations chosen below.

For the upper main sequence, $`M_V6`$, we evaluate $`F(M_V)`$ directly using the Hipparcos catalog (ESA, 1997). For example, the LF for $`M_V=4`$ would be computed by summing $`_i[(4/3)\pi D_i^3]^1`$ over all stars within the Hipparcos completeness limit, $`V<7.3`$, having $`3.5<M_V<4.5`$, and lying within 50 pc. The distance $`D_i`$ is the minimum of 50 pc and the distance at which the star would have $`V=7.3`$. Actually, we are not directly interested in the LF at $`M_V=4`$, but rather in the integral of $`L^{3/2}R^{7/2}`$ over the subpopulations that make up the $`M_V=4`$ bin of the LF. Hence

$$F(M_V=4)=\underset{3.5<M_{V,i}<4.5}{}\left(\frac{L_i}{L_0}\right)^{3/2}\left(\frac{R_i}{R_0}\right)^{7/2}\left(\frac{4}{3}\pi n_0D_i^3\right)^1.$$
(8)
The stellar radii are determined from Hipparcos/Tycho $`(B_T,V_T)`$ photometry and the color/surface-brightness relation of Gould & Morgan (2003), $`\mathrm{log}R/R_{}=0.597+0.536(B_TV_T)M_{V_T}/5`$, ultimately derived from van Belle (1999).

To normalize $`F`$, we adopt the values associated with $`M_V=5`$ stars: $`L_0=0.86L_{}`$, $`R_0=0.97R_{}`$, $`n_0=0.0025\mathrm{pc}^3`$.

The resulting function (Fig. 1) shows that the majority of stars that are probed will be F and G type ($`2M_V6.5`$). To check how this distribution depends on the size of the volume sampled, we recalculate the distribution function for the case in which the observed stars cover a much larger volume – one which would be better described by a thin disk, rather than a spherically uniform distribution. For this case we find that the scalings shown in equation (7) are replaced by $`F(M_V)nhLR^2`$, where $`h`$ is the scale height of each population. The distribution function is still dominated by F and G stars, but there are more K and early to mid M stars ($`7M_V12`$). Of course, it is common knowledge that magnitude-limited $`(FnL^{3/2})`$ samples of main-sequence stars will be dominated by F and G stars. The interesting feature of Figure 1 is that this result does not qualitatively change despite the addition of the factor $`R^{7/2}`$ in equation (7), which very strongly favors later-type stars.

## 3 RANDOM NOISE

Equation (6) describes what kinds of XSP systems can be detected by a certain survey, given a photometric detection limit. The threshold $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ is determined by taking account of the fact that the data stream from an XSP search must be analyzed for any combination of the parameters $`R`$, $`a`$, and $`r`$ within reasonable ranges. Such analysis will, however, yield a number of false-positive detections due to random noise. The threshold value of $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ must be chosen to yield a manageable number of candidate systems for follow-up observations.

To determine $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$, we generate 1000 independent simulated streams of photometric observations of a single system with a host star of one solar mass and a circular planetary orbit. We attempt to simulate a schedule that would be characteristic of an all-sky survey that re-images a given star approximately every few days with varying intervals between observations. For each data stream, we generate 1000 observations at irregular intervals over 1800 simulated days, and then phase the observations for a range of periods from 3.0 to 3.1 days and a range of transit lengths. The number of the observations and the duration of the simulated survey are arbitrarily chosen as plausible characteristics for the type of all-sky survey we envision.

The duration of the transit depends on the orientation of the system with respect to our line of sight, and we analyze the data for eight equally spaced, progressively more inclined orientations. For each orientation, we test for 16 equally spaced phases. We test the resulting data sets for transit-like dips in the light curve, which we define simply as intervals during which the local mean light curve dips significantly below the global average. This $`(8\times 16)`$ grid structure is chosen empirically: we find that the number of false positives increases linearly with grid density below this density and then flattens above it.

The result shows that for this kind of system, in order to restrict follow-up analysis to the 0.1% of the full sample most likely to yield a true transit detection, the value for $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ should be set to $``$ 36.6. We take the highest value of $`\mathrm{\Delta }\chi ^2`$ from each of the 1000 sets, and of those highest values, we then sort the 1000 sets from highest value of $`\mathrm{\Delta }\chi ^2`$ to lowest. In Fig. 2 (inset), we plot the highest value of $`\mathrm{\Delta }\chi ^2`$ for each of the 1000 data streams.

To check the robustness of this number, we run additional simulations of 25 data sets, each with different observation schedules. First, we generate a set of observations that are randomly distributed throughout a 1800 day mission, with the requirement that no two observations be less than 10 minutes apart. Then we regenerate the data set with observations that are evenly spaced throughout the project. We also conduct the analysis on data sets with evenly spaced pairs of observations, with the observations in each pair separated by 14.4 minutes; and then again with the pair-spacing at 43.2 minutes.

These different configurations test for different types of observing schedules. For instance, a ground-based all-sky survey would most likely observe a star at essentially random times throughout a project. On the other hand, a space-based mission with a slowly rotating telescope (similar to the Hipparcos satellite) would observe a star in pairs of observations separated by minutes. For various reasons there could be a certain regularity imposed on either the space-based or ground-based observations.

We find that the value for $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ does not depend strongly on the observing schedule. We determine this by taking the highest value of $`\mathrm{\Delta }\chi ^2`$ for each of the 25 independent data streams, sorting the highest values of $`\mathrm{\Delta }\chi ^2`$ from each set in the manner described above, and then plotting the results for each type of observing schedule (Fig. 2). We find that there is only a $``$ 10% variation in $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ for different schedules.

The value of $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ depends on the number of different parameter configurations that are tested. When we rerun our analysis using a period range of 3.0 to 3.5 days, we increase the size of our parameter space by a factor of 5. We expect that the value of $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ depends on the size of the parameter space, $`\mathrm{\Psi }`$, according to $`\mathrm{\Delta }\chi _{\mathrm{min}}^2(\mathrm{\Psi }_{new})\mathrm{\Delta }\chi _{\mathrm{min}}^2(\mathrm{\Psi }_{old})=2\mathrm{ln}(\mathrm{\Psi }_{new}/\mathrm{\Psi }_{old})`$. This prediction agrees with the simulations, which show that an increase in parameter space by a factor of 5 leads to an increase in $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ by $`2\mathrm{ln}(5)=3.2`$.

We then predict $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ if the size of the parameter space is expanded to search for transits with periods between 2 and 10 days (the range of periods most probable for the sort of fiducial values of the other parameters we have chosen). This would increase $`\mathrm{\Psi }`$ by a factor of 80 compared to the period range of 3.0 to 3.1 days. Therefore, $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ for the expanded range of periods is higher than 36.6 by $`2\mathrm{ln}(80)=8.8`$. So the value expected for $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ for such an analysis is about 45.

## 4 NUMBER OF SYSTEMS PROBED

The weak dependence of $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$on the observing schedule implies that the sensitivity of a project to planets essentially depends only on the total number of photons detected from each star, and not on the details of how they are collected. For a given star, this number is obviously proportional to the flux. We therefore characterize the sensitivity of the observing setup (telescope + detectors + duration + weather + etc.) by $`\gamma `$, the total number of photons that are detected from a fiducial $`V=10`$ mag star during the entire project. (Here we adopt $`V=10`$ as a reference, although the stars of interest lie in the range $`8V10`$.) We can utilize the various relations used to derive equation (6) to relate $`\gamma `$ to $`V_{\mathrm{max}}`$, the maximum apparent magnitude at which an equatorial transit can be detected,

$$V_{\mathrm{max}}=2.5\mathrm{log}\left(\frac{\mathrm{\Delta }\chi _{\mathrm{min}}^2\pi aR^3}{r^4\gamma }\right)+10.$$
(9)
Then, considering equations (5), (6), and (7), we obtain

$$N_p=730F(M_V)\left(\frac{a}{a_0}\right)^{5/2}\left(\frac{r}{r_0}\right)^6\left(\frac{\gamma }{\gamma _0}\right)^{3/2}\left(\frac{\mathrm{\Delta }\chi _{\mathrm{min}}^2}{45}\right)^{3/2},$$
(10)
where we have adopted $`\gamma _0=1.25\times 10^7`$, $`a_0=10R_{}`$, $`r_0=0.10R_{}`$, and where we have made our evaluation at $`M_V=5`$ (i.e. $`R=0.97R_{}`$, $`V_{\mathrm{max}}=10`$, $`d_0=100`$pc, and $`n=0.0025\mathrm{pc}^3`$). Note that $`\gamma _0=1.25\times 10^7`$ corresponds to approximately 625 20-second exposures with a 5 cm telescope and a broadened $`(V+R)`$ type filter for one $`V=10`$ mag fiducial star.

As mentioned above, a noteworthy feature of equation (10) is that $`N_p`$ depends on the characteristics of the survey primarily through the parameter $`\gamma `$. Moreover, since $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ depends only logarithmically on the size of the parameter space being explored, it plays a minimal role in survey design compared to the other variables in equation (10).

We envision two scenarios to which these results will be applicable. In one, a stream of photometric measurements from a space-based astrometric mission, such as GAIA, could be searched for transits. In the other, a ground-based survey could use one or more dedicated telescopes to search all bright stars in the solar neighborhood.

## 5 IMPLICATIONS FOR TELESCOPE DESIGN

We now apply the general analysis of § 2 and § 4 to the problem of optimizing telescope design for quickly locating a “large” number of bright ($`V10`$) transiting systems. Since only one such systems is now known, we define “large” as $`𝒪(10)`$. From equation (10) and the 0.75% frequency of hot jupiters measured from RV surveys, there are roughly 5 such systems to be discovered over the whole sky per magnitude bin at $`M_V=5`$ for $`V_{\mathrm{max}}=10`$. Hence, from Figure 1, of order 25 are to be discovered from all spectral types. It would, of course, be possible to discover even more by going fainter, but setting this relatively bright limit is advisable for three reasons. First, as we argued in the introduction, the brightest transits are the most interesting scientifically, and most of the transits detected in any survey will be close to the magnitude limit. Second, as we discuss below, a high dynamic range, $`\mathrm{\Delta }V=V_{\mathrm{max}}V_{\mathrm{min}}`$, can only be achieved at considerable cost to the observing efficiency. Hence, if high efficiency is to be maintained, setting $`V_{\mathrm{max}}`$ fainter means eliminating the brightest (most interesting) systems. Third, at $`V_{\mathrm{max}}=10`$, we are already reaching distances of 100 pc for G stars. Hence the number of transits observed in fainter surveys will not continue to grow as $`d_0^3`$ as in equation (6).

In previous sections, we ignored the loss of sensitivity to systems that are brighter than $`V_{\mathrm{min}}`$, which is set by saturation of the detector (or more precisely, by the flux at which detector non-linearities can no longer be accurately calibrated). This fraction is $`10^{0.6\mathrm{\Delta }V}`$, or 6% for $`\mathrm{\Delta }V=2`$, which we therefore adopt as a sensible goal. That is, we wish to optimize the telescope design for,

$$8=V_{\mathrm{min}}<V<V_{\mathrm{max}}=10.$$
(11)
(In any event, essentially all stars $`V<8`$ have already been surveyed for XSPs using RV, and the problem of determining which among the planet-bearers have transits is trivial compared to the problem of conducting an all-sky photometric variability survey.)

Optimization means maximizing the photon collection rate, $`\gamma /T`$, where $`T`$ is the duration of the experiment and $`\gamma `$ is, again, the total number of photons collected from a fiducial $`V=10`$ mag star. Explicitly,

$$\gamma =KTD^2\frac{(\mathrm{\Delta }\theta )^2}{4\pi },$$
(12)
where $`\mathrm{\Delta }\theta `$ is the angular size of the detector, $`D`$ is the diameter of the primary-optic, $``$ is the fraction of the time actually spent exposing, and $`K`$ is a constant that depends on the telescope, filter, and detector throughput. For our calculations, we assume $`K=K_040e^{}\mathrm{cm}^2\mathrm{s}^1`$, which is appropriate for a broad $`(V+R)`$ filter and the fiducial $`V=10`$ mag star. The design problems are brought into sharper relief by noting that $`\mathrm{\Delta }\theta =/D`$, where $``$ is the linear size of the detector and $``$ is the focal ratio, or $`f`$/#, of the optics. Equation (12) then becomes

$$\gamma =\frac{K^2T}{4\pi ^2}.$$
(13)
That is, almost regardless of other characteristics of the system, the camera should be made as fast as possible. We will adopt $`=1.8`$, below which it is substantially more difficult to fabricate optics. A more remarkable feature of equation (13) is that all explicit dependence on the size of the primary optic has vanished: a 1 cm telescope and an 8 m telescope would appear equally good! Actually, as we now show, there is a hidden dependence of $``$ on $`D`$, which favors small telescopes.

### 5.1 Considerations for Aperture Size

The global efficiency $``$ can be broken down into two factors, $`=_0_S`$, where $`_0`$ is the fraction of time available for observing (i.e., during which the sky is dark, the weather is good, etc.), and $`_S`$ is the fraction of this available observing time that the shutter is actually open. The first factor is not affected by telescope design and so will be ignored for the moment. The second factor should be maximized. The smaller the telescope aperture is, the longer the exposures can be before a $`V_{\mathrm{min}}=8`$ mag star saturates. Since the readout time is fixed, a smaller fraction of time is lost to read-out. We adopt as benchmarks a detector with pixel size of $`\mathrm{\Delta }x_p=9\mu `$m, well depths of $`10^5\mathrm{e}^{}`$, and a telescope diameter of 5 cm.

To make explicit calculations, $`\theta _{\mathrm{PSF}}`$, the full width at half maximum of the point spread function (PSF) must be specified. For the fast optics ($`<`$ few) we consider here, the diffraction limit is always much smaller than a pixel, regardless of aperture: $`\theta _{\mathrm{diff}}/\theta _p1.22\lambda /\mathrm{\Delta }x_p0.16`$ (for our fiducial choices). At the small apertures we will consider, the diffraction limit is larger than the seeing, so it is possible to make the PSF much smaller than a pixel, $`\theta _{\mathrm{PSF}}\theta _p`$. This would have the advantage of reducing sky noise and is a useful approach when it is possible to always center the telescope at the same field position as is the case for “point and stare” experiments. However, for an all-sky survey, which cycles through many fields, such precision repeat pointing is extremely difficult. Without it, precision photometry is impossible unless the sub-pixel response of the CCD is mapped out in detail. We therefore adopt a Nyquist-sampled PSF, for which the sky noise is approximately that falling on $`4\pi 13`$ pixels.

Our overall consideration for telescope design must take into account three factors. First, we with to maximize observing efficiency $`_S`$. Second, we wish to achieve the highest possible signal-to-noise ratio. Third, we must avoid any distortion problems with the optics. There are four effects through which aperture size can impact these factors. Two of these effects, observing efficiency and scintillation noise, will drive us to larger telescopes, while the other two effects, sky noise and focal plane distortion, will drive us to smaller telescopes. As we show below, for the observing parameters we have specified an aperture of 5 cm ensures a manageable (and unique) balance between the various effects.

#### 5.1.1 Exposure Time vs. Readout Time

Assuming Nyquist sampling, at most half the light from a point source falls within one pixel. We can directly calculate the ratio of time lost to readout $`T_{\mathrm{read}}`$ to the time spent exposing $`T_{\mathrm{exp}}`$,

$$\frac{T_{\mathrm{read}}}{T_{\mathrm{exp}}}=1\left(\frac{D}{5\mathrm{cm}}\right)^210^{0.4(V_{\mathrm{max}}10)}\left(\frac{W}{W_0}\right)^1\frac{K}{K_0}\frac{T_{\mathrm{read}}}{30\mathrm{s}}10^{0.4(\mathrm{\Delta }V2)},$$
(14)
where $`W`$ is the well depth of the detector pixels, and $`W_0`$ = $`10^5\mathrm{e}^{}`$ is a fiducial well depth. Note that the factor $`10^{0.4(V_{\mathrm{min}}10)}`$, which arises from the need to avoid saturation of the brightest stars (where $`V_{\mathrm{min}}=V_{\mathrm{max}}\mathrm{\Delta }V`$), has been broken up into two terms to permit easy comparisons of equation (14) with equations (16) and (18) below.

In order to maximize the efficiency $`_S`$, the fraction of observing time devoted to readout should be minimized, and therefore, according to equation (14), so should the aperture size. The telescope will operate reasonably efficiently so long as $`T_{\mathrm{read}}T_{\mathrm{exp}}`$.

#### 5.1.2 Scintillation

Another concern that arises for small apertures is the effect of atmospheric scintillation, which is characterized by (Young 1967; Warner 1988),

$$\frac{\mathrm{\Delta }I}{I}=S_0\left(\frac{D}{\mathrm{cm}}\right)^{2/3}X^{3/2}\mathrm{exp}\left(\frac{h}{h_0}\right)\left(\frac{2T_{\mathrm{exp}}}{\mathrm{sec}}\right)^{1/2},$$
(15)
where $`S_0`$ = 0.09, $`X`$ is the airmass, $`h`$ is the altitude of the observatory, and $`h_0=8`$ km is the scale height of the atmosphere. Using values of $`D=5`$ cm, $`X=1.5`$, $`h=2`$ km, and $`T_{\mathrm{exp}}=30`$ s, we find $`\mathrm{\Delta }I/I=0.0057`$.

Hence, for exposure times set by a saturation threshold, $`T_{\mathrm{exp}}D^2`$ (see eq. ),

$$\frac{\mathrm{Scintillation}\mathrm{Noise}}{\mathrm{Source}\mathrm{Noise}}=1\left(\frac{D}{5\mathrm{cm}}\right)^{1/3}10^{0.2(V10)}.$$
(16)
Therefore, despite the common perception that scintillation is a greater problem for smaller telescopes, for the fixed photon counts per exposure that are of interest in the present context, scintillation noise increases with increasing aperture. However, this dependence is fairly weak. For $`D=5`$ cm, the scintillation noise is just slightly smaller than the photon noise at $`V=10`$.

#### 5.1.3 Sky Noise

All the calculations in §§ 2, 3, and 4 have assumed that sky noise is negligible, i.e. $`VV_{\mathrm{sky}}`$ where $`V_{\mathrm{sky}}`$ is the light from the sky falling on 13 pixels (for Nyquist sampling). Assuming a somewhat conservative mean sky brightness of $`V=20.0\mathrm{mag}/\mathrm{arcsec}^2`$, one finds

$$V_{\mathrm{sky}}=10.6+5\mathrm{log}\left[\frac{D}{5\mathrm{cm}}\frac{}{1.8}\left(\frac{\mathrm{\Delta }x_p}{9\mu \mathrm{m}}\right)^1\right].$$
(17)
To determine how much of a problem sky noise will be, we can compare it to the amount of photon noise,

$$\frac{\mathrm{Sky}\mathrm{Noise}}{\mathrm{Source}\mathrm{Noise}}=\frac{3}{4}\left(\frac{D}{5\mathrm{cm}}\right)^110^{0.2(V10)}.$$
(18)
Thus, sky noise is less than photon noise for a 5 cm telescope, but would become a serious problem for a substantially smaller aperture.

#### 5.1.4 Focal Plane Distortion

Focal plane distortions toward the edge of the detector become difficult (and expensive) to correct when the field of view is too large. For example, for a $`4k\times 4k`$ detector with $`\mathrm{\Delta }x_p=9\mu `$m pixels,

$$\mathrm{\Delta }\theta =23^{}\left(\frac{D}{5\mathrm{cm}}\right)^1\left(\frac{}{1.8}\right)^1.$$
(19)
Again, this is manageable for $`D=5`$ cm but could would potentially be a problem for smaller apertures.

Although this field of view is not so large as to create focal plane distortions, it is important to note that for a field of view this large, the telescope must be placed on an equatorial mount. For the alternative, an alt-az mount, the rotation of the sky will cause stars at one edge of the field to move faster across the detector than stars at the opposite edge. For similar reasons, the telescope must track rather than using drift scan.

### 5.2 Optimal Telescope Design

The relationships described in § 5.1.1 through § 5.1.4 can be used to optimize the aperture $`D`$ for any survey parameterized by a given $`V_{\mathrm{max}}`$. One is driven to smaller apertures by the goals of minimizing scintillation noise and the fraction of time spent on readout, and to larger apertures by the goals of minimizing sky noise and field distortion. Given available $`=3.6`$cm $`4096\times 4096`$ detectors and reasonably fast $`=1.8`$ optics, a $`D=5`$cm telescope is optimal for an all-sky survey of $`V=10`$ stars. Among all existing transit programs of which we are aware, the WASP telescope ($`D=6`$ cm lens, $`=2.8`$ focal ratio, $`2k\times 2k`$, $`=3`$cm detector, Street et al. 2002) comes closest to meeting these design specifications.

We use equation (13) to determine the required duration of the experiment using our optimally-designed telescope, adopting $`\gamma =\gamma _0`$, $`K_0=40\mathrm{e}^{}\mathrm{cm}^2\mathrm{s}^1`$, $`_0=20\%`$, $`_S=50\%`$, $`=3.6`$cm, and $`=1.8`$. The total time required to conduct the survey assuming only source photon noise is $`T=4\mathrm{months}`$. Taking into account sky noise and scintillation increases that time by a factor of 2.5 to 10 months.

For surveys that intend to search for transits of stars fainter than the range we consider in this paper ($`8=V_{\mathrm{min}}<V<V_{\mathrm{max}}=10`$), the relations in § 5.1.1 through § 5.1.4 operate somewhat differently. For the magnitude range we are considering, equations (14) and (16) provide an upper limit on the aperture size, while equations (18) and (19) place a lower limit, in which both limits converge to $`D=5\mathrm{c}\mathrm{m}`$. At fainter magnitudes, the aperture size is restricted according to equations (14) and (18):

$$10^{0.2(V10)}\frac{D}{5\mathrm{cm}}10^{0.2(V_{\mathrm{max}}10)}.$$
(20)
This equation shows how the aperture size limits combine into a single scaling relation. What of the other two limiting factors? The issue of focal plane distortions is unimportant at fainter magnitudes, since $`\mathrm{\Delta }\theta <23^{}`$, which is outside of the regime where such distortions are a factor. Scintillation effects are also not a factor at fainter magnitudes, since the scintillation-noise restriction requires that $`(D/5\mathrm{cm})<10^{0.6(V10)}`$, which is a looser restriction than that of equation (20).

Therefore, for a survey of transits at a magnitude range $`8=V_{\mathrm{min}}<V<V_{\mathrm{max}}=10`$, the aperture size limits converge to $`D=5\mathrm{cm}`$, while at fainter magnitudes, the aperture size is given by equation (20).

### 5.3 Practical Implementation

In the previous section, we estimated the duration of observations required to achieve the minimum S/N to detect transits by hot Jupiters assuming certain fiducial parameters of a ground-based telescope, but calculated within the framework of a literal “all-sky” ($`4\pi `$) survey of randomly-timed observations that is more characteristic of satellite missions. The resulting estimate is useful for judging the viability of a given observing setup, but it glosses over a key issue in the detection of transits, namely the problem of folding the data. As discussed in §3, the number of folds (and hence the size of the search space that must be probed) scales directly as the duration of the experiment. This larger search space increases both the minimum $`\mathrm{\Delta }\chi ^2`$ for a robust detection and the amount of computing power needed to sift through the search space. The first effect is logarithmic in the size of the search space, so a 10-fold increase changes $`\mathrm{\Delta }\chi _{\mathrm{min}}^2`$ by only $`\mathrm{ln}102.3`$, which is well under 10%. However, the second effect is linear in the search space and so could easily overwhelm available computing resources if not carefully monitored. That is, there are important drivers for keeping the duration of observations of any given field to a minimum and hence for exploring the question of whether it is better to break up the “all-sky” survey into several smaller components, each of which could be completed in a shorter time. Indeed, the OGLE experiment (Udalski et al., 2002), the only transit experiment to successfully detect a transit (Konacki et al., 2003), was motivated by these considerations to concentrate its observations over a month duration so as to limit the number of foldings.

To estimate the true duration of the project required the achieve the minimum S/N, one must take account of two factors: First, as a practical matter, ground-based surveys from a single location can cover only an angular area $`\mathrm{\Omega }<4\pi `$. Second, during a year of continuous observations, any given patch of sky is observed only for about 6 months. These factors change the previous estimates from §5.2. The first factor implies that the observing time required to reach minimum S/N is actually lower than the time calculated from equation (13) by a factor of $`\mathrm{\Omega }/4\pi 0.5`$, since the project is observing about half the angular area on the sky. The impact of the second factor is more complex. In a given night, only half the accessible angular area $`\mathrm{\Omega }`$ can be observed. Therefore, each night the available area to observe is lower by an additional factor of 2, which means the rate of observations of a given point on the sky will double again. There are two scenarios at this point. The first scenario is that after accounting for the first factor, the time required for sufficient S/N is less than a year. In this case, the two factors combine to lower the estimate from §5.2 by a factor of about 4. In the second scenario, the time required for sufficient S/N is more than a year. In this case, after observing for 6 months the minimum S/N is not reached, and therefore the project will have to pick up again six more months later when the target is again visible, and so the second factor does not apply. As we see below, however, for the parameters we are considering, we are well within the regime of the first scenario.

Applying the first factor to the result in §5.2 reduces the calculated time from 10 months to 5 months. Since this latter duration is indeed smaller than 1 year, the second factor implies that the requisite S/N will actually be reached in 2.5 months. This is only slightly longer than the duration of the OGLE observations.

However, real experiments inevitably have larger errors than expected. (We mention one source of additional errors below.) If the errors prove sufficiently larger that the experiment requires more than 1 year (after application of the first factor) to achieve the minimum S/N, then one would not gain the advantage of the second factor. In this case, it would be better to break the sky up into strips by declination, and observe each strip for a year, so as to increase the S/N obtained in each strip during a single year, and so to permit the application of the second factor.

Another real-world consideration is that the photometric errors will not be Gaussian. For example, the Sloan Digital Sky Survey (SDSS) photometry errors, while Gaussian in their core, deteriorate to an exponential profile in the wings beginning at about $`3\sigma `$ (Ivezic, Z. et al., 2003). Such deterioration is likely to set in earlier for the small-aperture, wide-field, low-budget cameras that we envisage here. The non-Gaussian form of the errors has no practical impact on our analysis: there are so many data points that the central limit theorem guarantees that their combined behavior in each phase bin will be Gaussian (as we have implicitly assumed in §3. However, the non-Gaussian tails will tend to increase the $`\sigma `$ of the distribution relative to what would be inferred from the core, which of course will degrade the sensitivity of the experiment. If the errors are as well-behaved as those of SDSS, this problem can easily be resolved by the standard device of 3-$`\sigma `$ clipping. If not, then more complex strategies will be required. These are likely to be among the biggest practical problems facing the analysis, but in the absence of real data, they cannot be further analyzed here.

This work was supported in part by grant AST 02-01266 from the NSF.