A&A 391, 369-377 (2002)
DOI: 10.1051/0004-6361:20020802
G. Kovács ^{1} - S. Zucker^{2} - T. Mazeh ^{2}
1 - Konkoly Observatory, PO Box 67, 1525,
Budapest, Hungary
2 -
Wise Observatory, Tel Aviv University, Tel Aviv, 69978, Israel
Received 28 February 2002 / Accepted 4 April 2002
Abstract
We study the statistical characteristics of a
box-fitting algorithm to analyze stellar photometric time series in
the search for periodic transits by extrasolar planets. The algorithm
searches for signals characterized by a periodic alternation between
two discrete levels, with much less time spent at the lower
level. We present numerical as well as analytical results to predict
the possible detection significance at various signal parameters. It
is shown that the crucial parameter is the effective signal-to-noise
ratio - the expected depth of the transit divided by the standard
deviation of the measured photometric average within the transit.
When this parameter exceeds the value of 6 we can expect a significant
detection of the transit. We show that the box-fitting algorithm
performs better than other methods available in the astronomical
literature, especially for low signal-to-noise ratios.
Key words: methods: data analysis - stars: variables: general - stars: planetary systems - occultations
Application of variants of the PDM method in the analyses of variable star observations goes back to earlier times than that of the DFT. This is primarily because the PDM algorithm does not require the computation of trigonometric functions, which put a heavy load on the computers, especially in those early days. The most frequently cited implementation of the PDM idea is that of Stellingwerf (1978). However, earlier versions had appeared already in the '60s and early '70s, like those of Lafler & Kinman (1965, hereafter the L-K method), Jurkevich (1971) and Warner & Robinson (1972, hereafter the W-R method). Actually, it can be shown that up to a frequency- (or trial period-) independent constant, the method of Jurkevich (1971) is equivalent to that of W-R (see, Kovács 1980). Furthermore, without the additional feature of overlapping bin structure, the method of Stellingwerf (1978) is equivalent to that of Jurkevich (1971).
The study of the algorithm presented in this paper has been stimulated by the increasing interest in searching for periodic transits by extrasolar planets (e.g., Gilliland et al. 2000; Brown & Charbonneau 2000; Udalski et al. 2002), which follow the discovery of the transit of HD 209458 (Charbonneau et al. 2000; Henry et al. 2000) by its planetary companion (Mazeh et al. 2000). Due to the short duration of the transit relative to the orbital period (typically less than 5%), the signal expected is extremely non-sinusoidal. Considering the shallowness of the transit (typically less than 2% in the case of Jupiter-size planets) and the expected high noise level of the ground-based small telescopes capable of performing large-scale surveys (e.g., Borucki et al. 2001), we suggest an algorithm that utilizes the special form of the signal. The algorithm studied here (see also Gilliland et al. 2000 and Udalski et al. 2002) is based on direct Least Squares (LS) fits of step functions to the folded signal corresponding to various trial periods. We present numerical simulations as well as analytical considerations to estimate the ability of the algorithm to detect a faint signal in a noisy time series. We show that the algorithm performs significantly better and more efficiently than the published variants of PDM, DFT or some LS modification of the latter.
Let us denote the data set by , i=1,2, ..., n. Each x_{i} includes an additive zero-mean Gaussian noise with standard deviation. The noise is presented by assigning to each data point a weight w_{i}, defined as . It is further assumed that have a zero arithmetic average.
For a given trial period we consider a folded time series, which is a permutation of the original time series. This series is denoted by and the corresponding weights by . We fit a step function to the folded time series with the following parameters:
For any given (i_{1},i_{2}), we minimize the expression
(1) |
(2) |
(3) |
(4) |
(5) |
The most obvious meaning of follows immediately from Eqs. (1) and (4). At the maximum value of , of Eq. (4) is related to the average variance of the noise. By using the definition of for the transit depth, and the corresponding estimates of H and L of Eq. (2), we find that at the correct test period yields also an estimate of , i.e., .
In a practical implementation of the above procedure^{}, we suggest to divide the folded time series into m bins and evaluate by using these binned values. This approach is very efficient computationally, and yields an exact LS solution with a time resolution defined by the number of bins. Although the lower resolution affects the efficiency of the signal detection, in all interesting cases a good compromise can be made between computational constraints and the effectiveness of signal detection (see next section). One can further minimize the amount of computation by recognizing that for a given period each s(i_{1},i_{2}) and r(i_{1},i_{2}) can be obtained by adding the corresponding bin values to the already evaluated functions of smaller i_{1} and i_{2}.
Figure 1: An example of the power of the BLS method. The time series, the normalized BLS frequency spectrum and the folded time series are shown in the upper and lower two panels, respectively. The signal parameters are displayed in the header (see text for details). | |
Open with DEXTER |
In all subsequent tests the signals have a period of and span a timebase T of . The signal is sampled at times , where is a uniformly distributed random variable in [0,1]. This distribution corresponds to the timings we may obtain during a short but concentrated and continuous observational campaign from several ground-based observatories or from space. If the sampling or test periods are different from the ones used in this paper, the main results presented here still remain valid, assuming that the distribution of the data in the folded time series at any trial period is basically uniform. Our non-periodic sampling yields aliasing-free spectra in the frequency band of interest. We do not deal with aliasing effects originating from nearly periodic sampling, because this problem affects all period searching algorithms in a similar way.
In all our simulations, the search for the fractional transit length is performed in the [0.01,0.10] range. The upper limit is well above the expected maximum fractional transit length for planets (see Defaÿ et al. 2001), but is in the right range for detached binaries.
The frequency spectra are computed in the band, a range that contains frequency components of the true period from the 10th subharmonics up to the 3rd harmonics. Usually we use 4000 frequency steps, because, due to the specific signal shape, the line profiles are very narrow. Therefore, our algorithm requires a much finer sampling than DFT (see also Stellingwerf 1978).
To characterize the Signal Detection Efficiency we introduce
(see also Alcock et al. 2000):
(6) |
Figure 2: BLS spectra computed for two different transit phases for the same number of bins and signal parameters, shown in the header. | |
Open with DEXTER |
Figure 3: Dependence of the of the BLS spectra on the transit phase for noiseless signals with parameters shown in the header. | |
Open with DEXTER |
Figure 4: As in Fig. 3, but for a signals of shorter fractional transit length. | |
Open with DEXTER |
In order to have a better idea of the possible ranges of for different numbers of bins, we show in Fig. 3 the as a function of the fractional transit phase , where denotes the shift in time of the starting moment of the transit with respect to some arbitrary epoch. The transit length is chosen so that at m=50 the bin size is about half the transit length. The figure shows that even this number of bins displays relatively large variations of as a function of the transit phase, let alone the m=20 case. The high number of bins clearly yields high, more stable . However, the necessary increase in the CPU time is rather high (see later). On the other hand, the low number of bins may yield a fairly low , due to occasional partial coverage of the transit and wide bin size. The periodicities in come from the equidistant bin distribution and are given by (mq)^{-1}in units of qP_{0}.
The dependence on the transit phase is a function of the transit length. To show this function we plot in Fig. 4 the same dependence for a transit length of q=0.015, instead of 0.035 as in Fig. 3. We see now that the corresponding to m=20 is stable, although with lower values. This is so because for a small number of bins the full transit is included in one bin in most cases - yielding values independent of the transit phase. The strong dependence on the transit phase occurs when the bin size is comparable to the actual length of the transit. In both examples, a high enough number of bins yields higher, and more stable , which, in the presence of noise, means also a higher probability of signal detection (see Sect. 3.3).
One of the most important features of any search algorithm is its ability to distinguish between false and true signals. In our case, this translates to the estimation of the statistical significance of the highest peak in the BLS spectrum. To study this question we performed intensive numerical tests to derive the Probability Distribution Function (PDF) of in case the input data contain only Gaussian noise.
Figure 5: Probability distribution functions of for the BLS method at various numbers of test frequencies shown at the horizontal lines. Thick lines are for the empirical (numerical), thin lines are for the semi-theoretical results described in the text. | |
Open with DEXTER |
Our results show that the PDFs depend only on the number of trial frequencies n_{f}, and are immune (within the numerical stability of the simulations) to changing the number of data points or number of bins used by the BLS algorithm. In order to ensure reasonable numerical stability, we use in all cases 1500 realizations to estimate the corresponding PDF. The number of data points is set to 500, 1000, or 2000, while the number of bins is chosen to be 50 or 100. Figure 5 displays the empirical PDFs obtained for various numbers of trial frequencies.
One could expect the probability of finding outstanding peaks in a noisy spectrum to increase with the number of sampled frequencies, if the samples were statistically independent. We derive a theoretical PDF by using this assumption, and show that adjustment of this PDF to the empirical one is possible, with parameters different from the ones predicted by the assumption of large independent samples.
The computation of consists of two major steps: (a) selection of the largest at each trial period, (b) selection of - the largest of all the values computed for the n_{f} trial frequencies. In general we have , and therefore it is obvious that many of the tested values cannot be considered to be independent samples. Therefore, in this - admittedly not exact, but, as we shall show below, quite practical - approach, it is assumed that the distribution of can be approximated by the one obtained for independent samples of , where is to be determined by numerical simulations.
Assume, for the sake of simplicity, that all data points have the same
error, .
We take
independent values of the
random variable
,
identify the highest value
and compute .
For a large sample size we can
assume the
values of
and
are constant.
Therefore, one can use Eq. (6) to write the probability that exceeds a specified value X:
p | = | ||
= | |||
= | |||
= | (7) |
(8) |
= | |||
= | (9) |
(10) |
(11) |
Figure 6: Signal detection efficiency for noisy signals as a function of . The same bin number of 100 and fractional transit length of 0.033 are used in all simulations. | |
Open with DEXTER |
Figure 7: Differential signal detection efficiency for noisy signals as a function of . Shaded circles and filled squares show results for , and , respectively. The transit length is the same in all cases and is chosen to be relatively small, i.e., q=0.015. | |
Open with DEXTER |
Tests performed with other numbers of bins show that the above pattern remains the same, with the sole difference that for a larger number of bins reaches larger maxima, and therefore, the linear region between and 13 becomes steeper. It is important to note that the region of around 6 is critical in all cases, because of the separation between the stochastic and deterministic signal detections. As we could expect, this regime coincides with the limit of DDE (see Eq. (10)).
In the above simulations we fixed the fractional transit length at a relatively large value, in order to ensure reasonable resolution with 100 bins. With a shorter transit length, at the same number of bins, the region becomes fuzzier, with a somewhat higher average value of (at least for moderately smaller transit lengths). However, the global properties remain the same as described above.
The response of to changing the number of bins is illustrated in Fig. 7. We see that a substantial increase of can be gained in many cases by moving from 50 bins to 200 bins. Of course, for longer transit lengths, the gain is smaller.
Figure 8: As in Fig. 6, but for the W-R method with 100 bins. | |
Open with DEXTER |
Figure 9: Comparison of the BLS and W-R methods for a realization of a signal with parameters shown in the header (other parameters are standard, as given in Sect. 3). The uppermost panel shows the folded/binned time series (dots) with the period and the fit (continuous line) obtained by the BLS method. | |
Open with DEXTER |
The next example examines the L-K method. We recall that in
this method the squared differences between the bin-averages are
computed for each folded time series. For large bin numbers and
smoothly varying time series, this method should yield very small
(in the limiting case, zero) L-K statistics at the correct period.
Figure 10: As in Fig. 9, but for the L-K method. On the ordinate axis for the L-K method, SD stands for "standard deviation'', the corresponding statistics of the L-K method. Minima of SDindicate periodicities in the signal by the the L-K statistics. | |
Open with DEXTER |
As expected, DFT is also not preferable for analyzing signals with periodic transits of short duration. In the noiseless case, DFT yields very slowly decreasing power for the successive harmonics. Although this is not a good property for correct period identification, DFT shows a reasonable stability against noise. As shown in Fig. 11, some remnant power at (or close to) the original harmonics is still visible in the DFT spectrum, but a reliable identification is not possible even at this high SNR (). At an even higher SNR, corresponding to , the main component becomes the highest amplitude peak, with a simultaneous increase of the harmonics and with a still considerable contribution from other parts of the spectrum, originating from the noise.
Figure 11: As in Fig. 9, but for the DFT method. | |
Open with DEXTER |
Finally, one can attempt to use multifrequency LS Fourier fit (FLS) for better approximation of the signal shape and thereby increasing the in the case of the Fourier method. (We note that for data distribution leading to orthogonal Fourier base, the method of Defaÿ et al. 2001 is equivalent to this direct LS Fourier approach.) Because of the substantial increase in execution time for this method, we use a smaller number of data points and limit the order of the Fourier fit to 10. By plotting the standard deviation of the residuals, we obtain the result shown in Fig. 12. Although the performance of FLS could be improved by using more harmonics, FLS did, in general, a considerably poorer job at moderately high noise levels (similarly to the one shown in Fig. 12).
Figure 12: As in Fig. 9, but for the the Fourier-sum-fitting Least Squares method with 10 components. Minima in the FLS spectrum indicate periodicities in the signal. | |
Open with DEXTER |
n | m | ||
1000 | 50 | 1.0 | 6.5 |
100 | 1.6 | 4.1 | |
200 | 4.0 | 1.6 | |
500 | 19.1 | 0.3 | |
5000 | 50 | 4.3 | 7.6 |
100 | 5.0 | 6.6 | |
200 | 7.6 | 4.3 | |
500 | 22.6 | 1.5 | |
8000 | 500 | 25.0 | 2.1 |
Notes:
- 5000 test frequencies are computed. - A SUN Ultra 300 MHz machine is used with an optimized FORTRAN compiler. |
The algorithm studied here assumes only two levels of the periodic light curve. This assumption ignores all other features that are expected to appear in planetary transits. Thus, we ignore the gradual ingress and egress phases of the transit, which carry important information about the parameters of the planetary orbit (e.g., Sackett 1999). The lengths of these phases are short compared to the transit and thus they are not expected to affect significantly the results of the search. Another effect we ignore is the limb-darkening effect, which has indeed been shown to be small in the case of HD 209458 (e.g., Deeg et al. 2001). The effectiveness of the algorithm relies on the above simplifying assumption, which is justified as long as we are interested in a detection tool. After the periodicity is detected we can try to recover subtle features of the folded light curve, in order to derive the stellar and the planetary characteristics.
Our main interest is in cases where the signal-to-noise ratio is small, and one cannot identify the signal by monitoring a single transit, because the stellar drop in intensity is buried in the noise. Contrary to the search of transits by the HST in 47 Tuc (Gilliland et al. 2000), where the noise was small relative to the expected transit dip, we concentrate on cases in which the periodic signal can be detected only after many measurements are accumulated and the unknown transit is monitored many times. To be able to deal with a large number of observations, of a thousand or more, we have introduced binning into the folded data. We have shown that as long as the bin size is small compared to the expected transit length, the efficiency of the method is not affected.
One additional factor that determines the computational load of the algorithm is the range of transit length searched for. The maximum possible transit length can be estimated if we know the orbital, stellar and planetary radii. For a given stellar mass, the stellar radius can be derived by the mass-radius relation, and the orbital radius can be derived for any period. Recent theories give some estimates for the planetary radii. Therefore, for a given stellar mass we can estimate the maximum duration of the transit, which for HD 209456 is only a few percent of the period. For most ground-based and space searches for planetary transits one would have some idea of the stellar mass of all transit candidates, and therefore we can make our algorithm computationally more efficient by imposing a variable maximum duration on the transit length.
The significance of the detection depends primarily on the effective signal-to-noise ratio of the transit. The signal is the stellar brightness within the transit, relative to the brightness outside the transit, and the noise is the expected scatter of the measured average of the stellar brightness inside the transit. The scatter is composed, obviously, of the observational noise as well as of stochastic variation of the stellar intensity. It seems that the effective signal-to-noise ratio should exceed 6 in order to get a significant detection. This requirement should be taken into account when planning future searches for extrasolar planetary transits.
Acknowledgements
We acknowledge grants OTKA T-038437, T-026031 and T-030954 to G.K. and grant 00/40 from the Israeli Science Foundation.