# Modelling solar-like variability for the detection of Earth-like planetary transits

###### Key Words.:

planetary systems – methods: data analysis – techniques: photometric – stars: activity – stars: late-type^{†}

^{†}offprints: A. S. Bonomo

###### Abstract

Context:

Aims:We present a comparison of two methods of fitting solar-like variability to increase the efficiency of detection of Earth-like planetary transits across the disk of a Sun-like star. One of them is the harmonic fitting method that coupled with the BLS detection algorithm demonstrated the best performance during the first CoRoT blind test.

Methods: We apply a Monte Carlo approach by simulating a large number of light curves of duration 150 days for different values of planetary radius, orbital period, epoch of the first transit, and standard deviation of the photon shot noise. Stellar variability is assumed in all the cases to be given by the Total Solar Irradiance variations as observed close to the maximum of solar cycle 23. After fitting solar variability, transits are searched for by means of the BLS algorithm.

Results:We find that a model based on three point-like active regions is better suited than a best fit with a linear combination of 200 harmonic functions to reduce the impact of stellar microvariability provided that the standard deviation of the noise is times larger than the central depth of the transits. On the other hand, the 200-harmonic fit is better when the standard deviation of the noise is comparable to the transit depth.

Conclusions:Our results show the advantage of a model including a simple but physically motivated treatment of stellar microvariability for the detection of planetary transits when the standard deviation of the photon shot noise is greater than the transit depth and stellar variability is analogous to solar irradiance variations.

## 1 Introduction

Since the pioneering work by Mayor & Queloz (1995), about 270 extrasolar planets
have been discovered, mostly through the periodic oscillations of the radial
velocity of their parent stars^{1}^{1}1See http://exoplanet.eu/. About 35
show transits across the disk of their stars allowing us to remove
degeneracy between orbital inclination and amplitude of the radial velocity curve,
thus providing information on planet’s mass and radius. Transit searches from the ground are limited to
giant planets because their photometric precision cannot be pushed beyond the millimagnitude limit
due to the scintillation of the Earth atmosphere.
This limitation can be overcome by space-borne telescopes. The space missions
CoRoT and Kepler have been specifically designed for searching planetary transits and
can reach the photometric precision to detect transits of terrestrial planets
in front of late-type main-sequence stars (Baglin, 2003; Moutou et al., 2005; Borucki et al., 2004).

The main source of noise in the detection of terrestrial planet transits at the level of precision of space-borne photometry is the intrinsic variability of the parent star. For late-type stars having an outer convection zone, magnetic activity is the main source of variability and can severely limit the capability of detecting transits even in relatively inactive stars. In the case of the Sun, the transit of a large sunspot group across the disk produces a dip of the solar irradiance of parts per million (ppm), whereas the central transit of an earth-sized planet produces a flux decrease of only 84 ppm (Jenkins, 2002). During phases of minimum activity, even when no sizeable spot group is detected, the solar flux is modulated by faculae, with amplitudes usually exceeding 100-200 ppm.

In principle, the extraction of the transit signal from a light curve contaminated by stellar microvariability is possible thanks to the different timescales characterizing transit duration and stellar variability. Specifically, the former ranges from a few hours to about one day, in the case of terrestrial planets with orbital periods up to about one year around a solar-like star, while the latter is dominated by the rotational modulation of the largest active regions and their evolution, taking place on time scales from weeks to months in the case of the Sun (cf., e.g., Jenkins, 2002; Lanza et al., 2003, 2007, and references therein).

Several methods have been proposed to reduce the impact of solar-like variability on transit detection, starting from simple high-pass filters (Defaÿ et al., 2001) to optimal filters (Jenkins et al., 2002; Carpano et al., 2003). Those methods work well only for continuous, evenly sampled timeseries because gaps or irregular sampling introduce aliasing problems. Moreover, optimal filters need a sufficiently long timeseries to estimate the power spectrum of the intrinsic stellar variability from the data themselves in order to construct the filter. For a slowly rotating star as the Sun, this implies a time series of at least two or three months, given that most of the variability is on rotational timescales. Recently Aigrain & Irwin (2004) have proposed methods to overcome the difficulties due to irregular sampling or gaps, based on a least square fit of individual sine and cosine components or on the application of iterative non-linear filters.

The first CoRoT blind test (Moutou et al., 2005) represented a testbed for several different methods to treat solar-like variability and showed the potentiality of those methods in combination with suitable transit finding algorithms to detect transits in a quite large set of simulated light curves including also gaps and instrumental effects. The best method to remove stellar microvariability proved to be a best fit based on the use of 200 sine and cosine harmonic functions coupled with the BLS transit finder algorithm (Kovács et al., 2002), proposed by the Geneva team.

In the present paper, we propose a different method to treat stellar variability, based on our model of the Total Solar Irradiance (hereinafter TSI) variations (Lanza et al., 2003, 2007). Previous models do not exploit knowledge of the cause of the variation, that is the intrinsic evolution and the rotational modulation of the visibility of active regions, while our model includes such effects. The residuals obtained after subtracting our model have a maximum amplitude of 100-150 ppm and an almost Gaussian distribution (cf. Lanza et al., 2003, and Sect. 3.1), suggesting us to apply algorithms developed for the detection of planetary transits in the presence of a Gaussian white noise. The performance of the proposed approach will be compared with that of the 200-harmonic best-fit algorithm of the CoRoT blind test to establish advantages and limitations of the new method. Moreover, other methods will be briefly discussed in Sect. 5, deferring a detailed comparison to successive papers.

## 2 Light curve simulations

To assess the performance of the proposed method, we shall analyse an extensive set of simulated light curves. Each light curve includes flux variations produced by Sun-like activity, planetary transits and photon shot noise, respectively.

To simulate the effects of Sun-like activity, we choose a subset of the TSI variation of duration 170 days ranging from
January 1 to June 19 2001 as observed by the radiometers of the
VIRGO experiment on board of the SoHO satellite
(Fröhlich et al., 1995, 1997, see Fig. 1). Specifically,
we use the 2.0-level TSI time series with a time sampling of one value per hour, accessible
from the VIRGO web site at the Physikalisch-Meteorologisches Observatorium
Davos – World Radiation Center^{2}^{2}2http://www.pmodwrc.ch/pmod.php?topic=tsi
/virgo/proj_space_virgo#VIRGO_Radiometry. The standard deviation of the
hourly TSI data is 20 ppm.
Short quasi-periodic gaps, not exceeding a few hours, are present in the TSI time series and are retained in
our simulations as representative of the gaps expected for a space-borne
photometric experiment. Their total duration is 36 hours, i.e., 0.9 per cent of the time series
(see Fig. 1).

The considered time series falls close to the maximum of activity cycle 23 and shows flux variations up to 2600 ppm due to the transits of large active regions across the solar disk. Therefore, it is well suited to study the impact of the largest Sun-like variability on planetary transit detection. Note that the level of activity of the Sun at the maximum of the 11-yr cycle is close to the median activity level of field dwarf stars of spectral types F and G, as found by Schmitt (1997) by comparing their coronal X-ray fluxes.

We prefer to use real solar irradiance data instead of simulating stellar microvariability with the approaches proposed by, e.g., Aigrain et al. (2004) or Lanza et al. (2006). Actually, those models are based on extrapolations of the solar variability characteristics to the case of more active stars and cannot be tested observationally yet.

To simulate stellar variability for intervals of 150 days, such as those of the long runs of the space experiment CoRoT, we extract 150-d sections of the above TSI time series choosing starting points at random between 0 and 20 days. In such a way, we get different realizations of Sun-like activity reducing possible systematic effects associated with the variations contained in a fixed TSI subset (cf. Sect. 3.1).

The light curves of planetary transits are computed by means of the IDL procedure UTM (Deeg, 1999). We consider circular planetary orbits with periods of 5, 10, 25 and 50 days, and an inclination of , i.e., the planet transits across the centre of the disk of the star. The radius of the planet is varied between 1.0 and 2.0 R in steps of R. A radius of 1 R and a mass of 1 M are assumed for the star and a quadratic limb-darkening law is adopted for the stellar photosphere. Therefore, the duration of the transits ranges from 3.11 hours for an orbital period of 5 days up to 6.67 hours for a period of 50 days.

In such a way, a total of 20 different transit light curves has been simulated. Finally, we shift the initial transit to a random phase between 0.0 and 1.0, taken from a uniform distribution, and add stellar variability and photon
shot noise assuming that Poisson statistics can be well approximated by Gaussian statistics
at the photon count levels characteristic of planetary transit searches.
We consider four levels
of noise, i.e., with standard deviations of 100, 200, 300 and 1000 ppm, which approximately
correspond to the noise levels reached in one hour integration time for G2V stars of magnitudes , 13, 14 and 16 as observed by CoRoT, respectively^{3}^{3}3http://corotsol.obspm.fr/web-instrum/payload.param/index.html.
One hundred independent white Gaussian noise realizations are computed
for each value of the standard deviation and added to each of the noiseless light curves giving a total of
8,000 light curves with transits.

In a search for planetary transits it is important to estimate the number of false alarms produced by those random combinations of the noise that may approximate a transit feature (see Doyle et al., 2000, and Sect. 4). We determine the frequency of false alarms by analysing transitless light curves obtained by combining solar irradiance variations with 100 random sequences for each of the above noise standard deviations, for a total of 400 light curves.

## 3 Analysis of simulated light curves

The simulated light curves are analysed in two steps. First we fit stellar microvariability by applying our model based on three active regions (hereinafter 3-spot model; Lanza et al., 2003) and the model proposed by the Geneva team for the CoRoT blind test, based on a least-square fit with harmonic functions (hereinafter 200-harmonic fitting; Moutou et al., 2005). Secondly, we apply the Box-fitting Least Square algorithm (or BLS; Kovács et al., 2002) to search for transits in the time series of the residuals obtained with the two different variability models to see which one gives the best performance in terms of detected transits and reduced false alarms. The residuals of the transitless light curves are analysed in the same way to find the frequency of false alarms produced by noise.

### 3.1 The three-spot model of Sun-like variability

The model is described by Lanza et al. (2003), and its performance is studied in detail by Lanza et al. (2003, 2007) in the case of the Sun at different phases of the 11-yr cycle. Therefore, we limit ourselves to a brief introduction. The model makes use of three active regions, containing both cool spots and warm faculae, to fit the rotational modulation of the TSI, and a uniformly distributed background component to fit the variation of its mean level along the 11-yr cycle. The ratio of the area of the faculae to that of the sunspots in each active region is fixed.

The rotation period is assumed as a free parameter ranging from 23.0 to 33.5 days, and its value is determined, together with the other free parameters, by minimizing the (see Lanza et al., 2003). The inclination of the rotation axis is set at the solar value during the fitting process. Considering also the mean flux level, the areas and the co-ordinates of the three active regions, the number of free parameters is eleven.

The timescale of evolution of the active region pattern on the Sun is significantly shorter than the rotation period. The longest time interval that can be modelled with three stable active regions is 14 days (cf. Lanza et al., 2003). Hence, the time series is subdivided into subsets of 14 days, each of which is fitted with the model described above. The duration of a transit is at least one order of magnitude shorter than the timescale of variation of the active region configuration, ensuring that the model is not significantly affected by transit dips. Small discontinuities at the level of ppm are sometimes present at the matching points between successive 14-d best fits. Their impact on our results is minimized by the random choice of the origin of the TSI time series used to simulate stellar activity, ensuring that those features never fall at the same points of the solar variation series.

In Fig. 2, we plot with solid lines the distributions of the residuals of the best fit to the TSI from 6 February 1996 to 6 February 2007, distinguishing three time intervals corresponding to the minimum, intermediate and maximum levels of activity along solar cycle 23, respectively. The distributions are well approximated by Gaussians except in the tails that appear slightly higher than a normal distribution. Power spectra of the residuals are plotted in Fig. 3 (left panel) revealing a characteristic inverse dependence on frequency up to 0.25 hr (Jenkins, 2002) and an almost white spectrum at higher frequencies. Note the high level of power during the intermediate phase of cycle 23 owing to the appearance of large sunspot groups in the decaying phase of the cycle.

### 3.2 Modelling Sun-like variability by 200-harmonic fitting

Stellar variability is modelled as a sum of harmonics of a fundamental frequency , where is the whole duration of each time series, i.e., days. The number of considered harmonics is fixed at , so the highest frequency in the model corresponds to a period of 1.5 days, well below the frequency range where most of the transit power is concentrated. Therefore, the fitting technique is not expected to greatly affect the transit signal. The model flux at time is given by:

(1) |

Including the average level , the model has 401 free parameters that can be determined by a least-square fit. We applied Singular Value Decomposition (hereinafter SVD) to solve the system of normal equations for this least-square problem, as described by Press et al. (1992).

The distributions of the residuals are well approximated by Gaussians except in the tails that appear slightly higher than a normal distribution (see Fig. 2, dashed green plots). Their standard deviations are smaller than those of the residuals obtained with the 3-spot model. Power spectra of the residuals are plotted in Fig. 3 (right panel) and are similar to those obtained with the 3-spot model, except for the cutoff at a frequency corresponding to the period of 1.5 days.

### 3.3 Searching for planetary transits

The time series of the residuals obtained by the above fitting approaches are analysed to detect transits by means of the BLS algorithm. It assumes that transits are strictly periodic with a period and that the noiseless time series assumes only two values, i.e., outside of transits and inside transits. In other words, each transit is assumed to have a simple box-shaped profile, that is enough for detection purposes. The time spent in transit at the level is , where the fractionary transit length is a small number . The time of the first transit in the time series is indicated by . In the presence of a white Gaussian noise with a constant standard deviation , assuming that the average level of the time series is zero, the algorithm aims to find the best model by an estimation of the four free parameters , , , and that minimize the between the model and the data. We refer the reader to Kovács et al. (2002) for further details. Here we note only that a useful measure of the signal-to-noise ratio of a transit detection is given by:

(2) |

where is the number of points in the whole time series (i.e., in our 150-day time series with one-hour sampling), and is the standard deviation of the average of all measurements during transits. According to Kovács et al. (2002), a transit detection can be considered significant at the 3- level against false alarms if . In the case of our tests, we redetermine appropriate thresholds (see Sect. 4).

The implementation of the BLS algorithm is performed by means of the Fortran subroutine accessible at http://www.konkoly.hu/staff/kovacs.html. We slightly modified it to avoid negative values of that would correspond to periodic positive spikes, retaining only best fits with a positive value of that correspond to transit dips. The parameter space is searched for transits with a period between 1 and 150 days, with 30,000 values of the frequency uniformly spaced between and day; the fractionary length of the transit is varied between 0.002 and 0.08, assuming a number of 700 phase bins in the light curve folded at each trial period. Such a sampling of the parameter space ensures that none of the simulated transits will get lost because of an inappropriate choice of the free parameters in our search (for an optimal choice of the free parameter steps and ranges see, e.g., Kovács et al., 2002; Doyle et al., 2000; Jenkins et al., 2002).

The large computational load of our experiment is managed by running our analyses on a large grid infrastructure consisting of three high-performance computing nodes located in Catania, Palermo and Messina. The total number of CPUs is 560, each one being an AMD Opteron dual-core 64-bit processor with 1 or 2 GB of RAM (Becciani, 2007). The CPU time for fitting solar variability and performing transit search is on the average minutes per light curve with the 3-spot method and minutes per light curve with the 200-harmonic method, while about 2.5 days of elapsed time have been necessary to analyse our complete set of 8,000 light curves.

## 4 Results

A typical sample of our results is given in Fig. 4, where the distributions of the values of are plotted for simulations with a planet of R, the labelled values of the orbital period , and a standard deviation of the noise ppm. The vertical dash-dotted lines indicate the false-alarm thresholds for the 200-harmonic fitting and for the 3-spot model , respectively. They are derived from the analysis of the transitless light curves under the requisite that the frequency of false alarms be (i.e, not exceeding one in a hundred).

The histograms in Fig. 4 show the distributions of for the cases where the orbital period of the planet was correctly identified or incorrectly identified, respectively. For a correct identification of the period, we require that the orbital period given by the BLS routine is within days from the true period. When the period is correctly identified and , we count the detection as a positive one. Note that very few cases with and outside of the acceptable range are visible in the panels of Fig. 4. These are not counted as positive detections.

The frequencies of positive detections obtained for different values of planetary radius, orbital period and noise level are reported in Table LABEL:det_freq, omitting records for values of giving no detections, except for the cases with the smallest orbital periods. In the first column we list the planetary radius , in the second the standard deviation of the Gaussian photon shot noise , in the third the orbital period , in the fourth the frequency of positive detections with the 200-harmonic method and in the fifth the frequency of positive detections with the 3-spot method , respectively. We note that for a noise standard deviation of 100 ppm, the 200-harmonic method is better than the 3-spot method, especially for planetary radii R. On the other hand, for higher noise levels ( ppm), the 3-spot method gives a better performance than the 200-harmonic method, especially for planetary radii R.

As an illustrative application of our results, let us consider, e.g., a G2V star of observed in the CoRoT white-light exochannel with an integration time of one hour. In this case, ppm and planets with radii as small as R and orbital period days should be detectable in more than 95 per cent of the cases if stellar variability and photon shot noise are the only effects limiting photometric precision. However, a proper assessment of the performance of a given space instrument should consider also other effects (see Bordé et al., 2003; Moutou et al., 2005), which goes beyond the purpose of the present study aimed at comparing different stellar activity fitting methods.

## 5 Discussion

The results presented above prompted us to investigate the cause of the different performance of the two methods at different levels of noise and planetary radii. We found that the 200-harmonic method, having 401 degrees of freedom to fit the entire 150-d sequence, is better suited to reproduce the flux perturbations due to magnetic activity and noise than the 3-spot method, which has only 11 free parameters for each 14-d best fit (cf. Figs. 2 and 3). However, the use of a set of orthogonal functions to fit stellar variability produces a significant distortion of the planetary transits in the residual light curves. In particular, the depth of the transits is reduced, as it is clear by considering the best fits obtained in the case of a light curve containing only transits (cf. Fig. 5). This is a consequence of the well-known Gibbs’ phenomenon (e.g., Morse & Feshbach, 1954). Since the transit dip is mathematically similar to a couple of close discontinuities of the first kind, the fitting function, which is a truncated series of orthogonal trigonometric functions, overshoots the flux variation at the ingress and at the bottom of the transit leading to a reduction of its depth in the residual curve. Note that the problem is not alleviated by increasing the order of the fitting functions, because it is a consequence of the orthogonal nature of the adopted set of fitting functions (cf. Morse & Feshbach, 1954). When the standard deviation of the noise is small ( ppm), the reduction of does not greatly affect the detection efficiency because , although reduced (cf. Eq. 2), is still above the significance threshold. Hence, the 200-harmonic method performs better than the 3-spot method since it fits better the light variations. On the other hand, when the noise level increases ( ppm), the reduction of has a larger impact on the -statistics as more and more values of fall below the threshold making the performance of the 200-harmonic method worse than that of the 3-spot method (see Fig. 6). Note that the 3-spot method is based on the use of only three fitting functions that are not orthogonal, thus it is not affected by the Gibbs’ phenomenon. Moreover, the reduction of the transit depth in the residual light curve is very small because the fitting functions mainly follow the points outside transits, especially when the transit duration is short and the orbital period is long with respect to the adopted 14-d time interval during which the parameters controlling the shape of the fitting functions are kept fixed.

The Gibbs’ phenomenon affects also other filtering methods proposed to reduce the impact of stellar microvariability, such as, e.g., that of Carpano et al. (2003), or the Wiener-like discrete filters by Aigrain & Irwin (2004).

Wavelet-based methods (e.g., Jenkins, 2002) may be useful to overcome the problems related to Gibbs’ phenomenon, given the non-orthogonal nature of their basis functions, but they are negatively affected by gaps in the time series. On the other hand, the non-linear iterative filter proposed by Aigrain & Irwin (2004) is pratically insensitive to gaps or irregular sampling and does not make use of orthogonal functions. Therefore, it is worth a detailed comparison with our method and this will be the subject of a forthcoming work.

Note that the 3-spot method performs well in the presence of gaps (see Lanza et al., 2003, 2007) because the model is naturally interpolated through them to follow the observations. The main limitations of the 3-spot model are the computing time and the necessity of adjusting some of the parameters if the star is not analogous to the Sun. The computation time is really not a problem when a parallel high-performance computing system is available, as we have demonstrated through the extensive use of a grid-based system to perform our calculations. As a matter of fact, the problem of analysing the large number of light curves coming from a space experiment is a fully-parallelizable problem and the elapsed time scales with the inverse of the number of available CPUs.

## 6 Conclusions

We have performed extensive numerical experiments to compare two different methods of fitting Sun-like variability for the detection of Earth-like planetary transits. Our results show that the 3-spot model performs better than the 200-harmonic fitting when the standard deviation of the photon shot noise is times greater than the central depth of the transit. Conversely, the 200-harmonic fitting has a better performance when the noise standard deviation is comparable with the depth of the transit. This behaviour has been interpreted as a consequence of the Gibbs’ phenomenon and is expected to be common to all filtering and fitting methods based on sets of orthogonal functions. On the other hand, the assumption of some general and simple physical hypotheses about the origin of Sun-like variability allows the 3-spot method to reach a truly appealing performance, in spite of the fact that it has only eleven degrees of freedom.

Therefore, it is worth to conduct further detailed comparisons, in particular with non-linear iterative filters, such as that proposed by Aigrain & Irwin (2004). We plan to perform such tests by simultaneously applying different methods to the same, large set of simulated light curves, thus extending the Monte Carlo approach of the present paper.

###### Acknowledgements.

The authors are grateful to the Referee, Dr. S. Aigrain, for her comments on a previous version of the present paper that greatly helped to improve their work. They also gratefully acknowledge Drs. U. Becciani, A. Costa and A. Grillo and the system managers of the Trigrid and Cometa Consortia for their technical advice and kind assistance during the implementation and running of our numerical experiments on grid-based high performance computing systems. ASB and AFL are also grateful to Drs. R. Alonso, P. Barge, P. Bordé, R. Cautain, M. Deluil, A. Leger, C. Moutou, for useful discussions on several aspects of the present work. The availability of unpublished data of the VIRGO Experiment on the cooperative ESA/NASA Mission SoHO from the VIRGO Team through PMOD/WRC, Davos, Switzerland, is gratefully acknowledged. The authors gratefully acknowledge support from the Italian Space Agency (ASI) under contract ASI/INAF I/015/07/0, work package 3170. This research has made use of results produced by the PI2S2 Project managed by the Consorzio COMETA, a project co-funded by the Italian Ministry of University and Research (MUR) within the Piano Operativo Nazionale ”Ricerca Scientifica, Sviluppo Tecnologico, Alta Formazione” (PON 2000-2006). More information is available at http://www.pi2s2.it and http://www.consorzio-cometa.it. Active star research and exoplanetary studies at INAF-Catania Astrophysical Observatory and the Department of Physics and Astronomy of Catania University is funded by MUR (Ministero dell’Università e Ricerca), and by Regione Siciliana, whose financial support is gratefully acknowledged. This research has made use of the ADS-CDS databases, operated at the CDS, Strasbourg, France.## References

- Aigrain et al. (2004) Aigrain, S., Favata, F., Gilmore, G. 2004, A&A, 414, 1139
- Aigrain & Irwin (2004) Aigrain, S., Irwin, M. 2004, MNRAS, 350, 331
- Baglin (2003) Baglin, A. 2003, Adv. Sp. Res., 31, 345
- Becciani (2007) Becciani, U. 2007, in Grid-enabled Astrophysics, L. Benacchio, F. Pasian, Eds., Edizioni Polimetrica; p. 179
- Bordé et al. (2003) Bordé, P., Rouan, D., Léger, A. 2003, A&A, 405, 1137
- Borucki et al. (2004) Borucki, W., Koch, D., Boss, A., et al. 2004, in Second Eddington Workshop: Stellar structure and habitable planet finding, ed. F. Favata, S. Aigrain, & A. Wilson, ESA SP-538, 177
- Carpano et al. (2003) Carpano, S., Aigrain, S., Favata, F. 2003, A&A, 401, 743
- Deeg (1999) Deeg, H. 1999, http://www.iac.es/galeria/hdeeg/
- Defaÿ et al. (2001) Defaÿ C., Deleuil, M., Barge, P. 2001, A&A, 365, 330
- Doyle et al. (2000) Doyle, L. R., Deeg, H. J., Kozhevnikov, V. P., Oetiker, B., Martín, E. L., et al. 2000, ApJ, 535, 338
- Fröhlich et al. (1997) Fröhlich, C., Crommelynck, D., Wehrli, C., et al. 1997 Sol. Phys. 175, 267
- Fröhlich et al. (1995) Fröhlich, C., Romero, J., Roth, H., et al. 1995, Sol. Phys. 162, 101
- Jenkins (2002) Jenkins, J. M., 2002, ApJ, 575, 493
- Jenkins et al. (2002) Jenkins, J. M., Caldwell, D. A., Borucki, W. J. 2002, ApJ, 564, 495
- Kovács et al. (2002) Kovács, G., Zucker, S., Mazeh, T. 2002, A&A, 391, 369
- Lanza et al. (2003) Lanza, A. F., Rodonò, M., Pagano, I., Barge, P., Llebaria, A. 2003, A&A, 403, 1135
- Lanza et al. (2004) Lanza, A. F., Rodonò, M., Pagano, I. 2004, A&A, 425, 707
- Lanza et al. (2006) Lanza, A. F., Messina, S., Pagano, I., Rodonò, M. 2006, AN, 327, 21
- Lanza et al. (2007) Lanza, A. F., Bonomo, A. S., Rodonò, M. 2007, A&A, 464, 741
- Mayor & Queloz (1995) Mayor M., Queloz, D. 1995, Nature, 378, 355
- Morse & Feshbach (1954) Morse, P. M., Feshbach, H. 1954, Methods of Theoretical Physics, McGraw-Hill Book Co. Inc., New York
- Moutou et al. (2005) Moutou, C., Pont, F., Barge, P., Aigran, S., Auvergne, M., et al., 2005, A&A, 437, 355
- Press et al. (1992) Press, W. H., Teukolsky, S. A., Vetterling, W. T., Flannery, B. P. 1992, Numerical recipes in Fortran, 2nd ed. (Cambridge Univ. Press); § 15.4
- Schmitt (1997) Schmitt, J. H. M. M. 1997, A&A, 318, 215