A&A 450, 1231-1237 (2006)
Research and Scientific Support Department, European Space Agency, ESTEC, Keplerlaan 1, 2200 AG Noordwijk ZH, The Netherlands
Received 31 October 2005 / Accepted 29 December 2005
Easy to use analytical formulae are presented for the computation the of light curves of extra-solar planetary transits. The equations are a function of the fractional radii of the planet and the parent star, the inclination of the orbit, and the limb-darkening coefficients of the star. Light curves can be solved for these parameters depending on the precision of the available observations. When the radial velocity curve is also available, as is normally the case to ensure the nature of the system, the masses, radii, and average density of both the star and the planet can be determined. The equations are valid for any degree of limb darkening, as well as for any type of transit. The cases of eccentric orbits, third light, or a non-zero relative luminosity of the planet can be easily taken into account. The basic assumption is that the projections of both the star and the planet on the plane of the sky are well represented by circular discs. The effects in case this assumption is not valid are also discussed. Practical applications are shown, beginning with the light curve of the photometrically discovered planet OGLE-TR-113, obtained with a ground-based telescope. As a second example, results are shown from the study of the light curve obtained for the transit of the giant planet in HD 209458 with the Hubble Space Telescope. Procedures to get the best fit parameters are briefly discussed.
Key words: techniques: photometric - stars: binaries: eclipsing - stars: planetary systems
During the past decade the study of extra-solar planets has been a very active field of research, and the method to find new cases based on planetary transits has been found to be very promising. Despite the dominant role of the radial velocity technique in the discovery of new planets, the possibility of detecting the small light dips photometrically due to planets moving in front of their parent stars has been considered for a long time. In fact the photometric method is expected to be less limited than the radial velocity techniques in the detection of Earth-like planets. For a review of the methods of detecting extra-solar planets, see e.g. Perryman (2000).
After the measurement by Charbonneau et al. (2000) with a rather small telescope, of the light variations caused by the giant planet around HD 209485, the method of planetary transits was confirmed and found to be a very rewarding source of information. This planet had been previously discovered using radial velocity changes but the evidence of the transit encouraged new searches using photometric techniques. New observations could of course benefit from the infrastructure developed for wide-field photometric surveys like those searching for microlensing events. In fact, existing instrumentation was immediately applied for such a purpose, and the success of it was soon shown by the discovery of several planets, e.g. within the OGLE-III survey (Udalski et al. 2003). The possibility of measuring accurate relative radii and orbital inclinations, leading to the actual masses, radii, and densities of the star and the planet, became a new and very interesting area of research that was not possible with only the radial velocities. The photometric precision required for these studies was found to be achievable with medium sized telescopes.
In addition to these efforts using ground-based telescopes, an excellent light curve was obtained with the Hubble Space Telescope for the first known case of HD 209458 (Brown et al. 2001). The photometric quality of the light curve required equally precise equations for its analysis. Moreover, space missions specifically designed to obtain very accurate light curves of planetary transits will also be available soon. The French-led space mission Corot, to be launched in 2006, will search for transits of extra-solar planets with the necessary photometric precision for detecting big Earths while Kepler, to be launched two years later, will get even more precise light curves, eventually leading to the detection of truly Earth-like planets. Both missions make use of the transit method to detect new planets and will provide a large number of light curves to be analyzed.
The study of the light curves of extra-solar planets has been generally based on relatively simple geometric algorithms for computing the occulted area of the stellar disc as caused by an opaque moving circle representing the planet. Deeg et al. (2001) adopted such a geometrical approach to put boundary constraints to the system elements and used detailed stellar models by Claret et al. (1995) to interpret the observational data. A similar procedure was applied by Seager & Mallén-Ornelas (2003) searching for a unique solution of planet and star parameters. In these studies it became clear that a correct treatment of the limb darkening of the star is very important in the computation of the system parameters. A more elaborated approach was presented by Mandel & Agol (2002) allowing the use of non-linear laws for limb darkening and a precise geometrical solution for the transits. Unfortunately the equations given by these authors are different for each type of transit and limb-darkening law, making the analysis and interpretation of the observations somewhat complicated.
Only rarely available synthetic light curve codes, developed for the analysis of eclipsing binaries, have been used to study planetary transits. Among these codes, EBOP (Etzel 1993) is probably the best choice since it is based on a geometrical solution for almost spherical stars as described by Nelson & Davis (1972). More sophisticated models, like the one developed by Wilson & Devinney (1971) for the analysis of close binaries, are clearly not well-suited for spherical components and the expected small value of the ratio of radii. An exception could be the ELC code by Orosz & Hauschildt (2000) due to the possibility of increasing the effective resolution in the surface elements through a Monte Carlo sampling, as indicated by Wittenmyer et al. (2005). The use of Roche geometry and detailed model atmospheres for highly-distorted stars is nevertheless not required in well-detached systems like those of planetary transits, contrary to the case of cool giants in close binaries originally targeted for the development of the ELC code.
In this paper we present analytical solutions for the transit problem adapted from those derived initially for the study of the light curves of eclipsing binaries by Kopal (1977). The original equations were developed in a long series of papers finally collected in a book (Kopal 1979) and are based on the study of the eclipses through the cross-correlation of optical apertures. This approach is very rigorous mathematically and provides a consistent method to derive the light-curve elements. It was intended to avoid the obtention of a "right combination of wrong elements'', which used to be a problem for synthetic light curve methods.
Unfortunately, the new equations did not get very much attention due to the fact that they were difficult to extend to realistic cases of close binaries, with highly-distorted components. On the other hand, the synthetic light curve codes benefited from a rapid increase in computing capability, leading to better and more effective parametric searches and thus avoiding the intrinsic problem of finding a unique solution. Another source of confusion with Kopal's equations certainly was the different misprints in several of the published results, the complex procedure to check them, and the difficulty of reproducing observational data. But the situation of essentially spherical components is actually what we were looking for in the case of extra-solar transiting planets and the equations can be effectively brought back into use to open new possibilities for the analysis of light curves as described in this paper.
The total luminosity l, as seen from Earth, is given at any time t by , where and represent the luminosity of the star and the planet, respectively, while denotes the relative loss of light due to the partial occultation of the star by the planet during the transit. Using the orbital phase , as reckoned from the primary conjunction t0, to measure the time - i.e. , with P representing the orbital period - and assuming the luminosities to be constant and normalized so that , we can re-write the light-curve equation as .
The main condition for the validity of the equations presented in this paper is that both the star and the planet be treated as circular discs moving in front of each other. This is indeed a realistic approximation since the mutual distortions caused by rotation and tides can be neglected for most extra-solar planetary systems. The assumption, which is of course also linked to the constant value of the luminosities, can be easily checked with the observations outside of the transit within the level of observational precision of the data. For the phases when there is no transit, we should obviously have and . Some comments about the effects of this assumption not being valid are included in Sect. 4.
The conjunction opposite to the transit, at , can moreover be used to directly measure the relative luminosity of the planet, . At this phase we have a total eclipse, and . If this expectedly very small light variation can be detected, the corresponding value of the luminosity can be adopted in the analysis of the transit. In most cases, the relative luminosity of the planet with respect to the star will be negligible, therefore can be safely used. It should be noted that, very recently, the luminosity of the planet in HD 209458 was detected using this method (Deming et al. 2005) in the infrared domain at m. The obtained value of is quite small but significant. The corresponding effect in the optical domain is actually consistent with . This has been confirmed with measurements using the small Canadian space mission MOST, in the optical, and in the near infrared with K-band ground-based observations (Snellen 2005).
Another interesting point is how to treat the light curves a with certain amount of third light. This is the contribution of background sources, which is found to be quite significant in the new searches for planetary transits in very crowded fields. Then, the normalization of the luminosities becomes , where L3 denotes the light of the additional contributions, and even if the depth of the transit will also be a function of the amount of third light, since . The more background light that is included, the shallower the changes in the total luminosity for the same geometrical parameters, making the process of solving it for the elements more difficult.
The whole problem of analyzing the light curve of a transiting planet is reduced to the
obtention of an analytical expression for ,
the fractional loss of light, which
will be given at any phase of the orbit by
where the integral is extended over the entire area S of the star occulted by the planet, I stands for the distribution of brightness over the star, of surface element d, and where is the so-called angle of foreshortening, i.e. the angle between the surface normal and the line of sight.
In order to proceed with the evaluation of Eq. (1) we need to have a good knowledge of the
brightness distribution over the surface of the star. Fortunately, the effect of limb
darkening is a well-known element in the theory of stellar atmospheres. The simplest
solution is to adopt a linear law of the form,
All the involved limb-darkening coefficients have been extensively tabulated (Diaz-Cordovés et al. 1995; Claret et al. 1995; Claret 2000). When analyzing light curves of planetary transits, it is important to remember that the best coefficients are those fitted by least squares to the model atmosphere rather than methods based on the conservation of the total flux. Claret (2000) has published theoretically derived values of the coefficients for 12 photometric bands in the optical and near-infrared (uvby, UBV, and RIJHK) for a wide range of temperatures, surface gravities, metal abundances, and micro-turbulence velocities, using available model atmospheres (ATLAS and PHOENIX).
In this paper, a general law of limb darkening is adopted as given by,
Using Eq. (4) for ,
it can be shown (see Kopal 1979) that associated
functions can be defined by means of
It is now possible to evaluate the
functions by integrating the
occulted area of the star. Figure 1 shows the basic parameters of the problem. The
radii of the star and the planet are represented, respectively, by
as fractional values of their separation, while the orbital inclination i is measured with respect to the plane of the sky. The apparent separation between
the centers of the projected discs of the star and the planet is given by ,
also expressed relative to the radius of the orbit. With these parameters,
we can establish the limits of integration for Eq. (1).
|Figure 1: Geometry and parameters of a planetary transit from the perspective of the observer.|
|Open with DEXTER|
In the case of circular orbits, we obviously have
Deeg et al. (2001) indeed produced an approximation to Eq. (1) using linear limb-darkening, circular orbits, and a small ratio of radii, . In this approximation, the surface brightness distribution can be considered to be constant within the area occulted by the planet and an analytical equation can be obtained. This was found to be very useful for ground-based observations obtained with relatively small telescopes, but certainly not for the more accurate light curves obtained from space. For this latter case, we mentioned the work by Mandel & Agol (2002) that extended the equations of Deeg et al. to second-order limb-darkening laws and solved Eq. (1) by means of elliptical integrals, under different assumptions for the law of limb darkening and cases of transit (involving different limits of integration).
A solution for the functions can also be obtained through numerical integration. This is essentially what is done in the analysis of the light curves of eclipsing binaries by synthetic models. It has already been mentioned that the code by Etzel (1993) could be used for this purpose. In fact, we find that EBOP provides good results and converges rapidly towards the right combination of elements. Because the original version of the code only considers a linear law of limb darkening, a modified version has been used as described by Giménez & Diaz-Cordovés (1993), where the possibility of computing light curves with non-linear limb-darkening laws was introduced. The main possible limitations of the method are the precision of the computed luminosities, which is a function of the adopted size of the integration rings, and the optimization of the elements by means of differential corrections. Due to the consideration of the surface brightness as constant over each of the rings and the relative size of the transiting planets, the finest grid available of 1 degree has to be used. This produces precisions in luminosities that are more than enough for ground-based observations, as well as for high-quality space-based light curves. Unfortunately, the differential corrections method inherent to the code was found to be non-applicable in most transits observed with ground-based telescopes since convergence is not effectively achieved for observational precisions poorer than 1/10 the depth of the transit.
In this paper, the goal is to derive a more general equation, and one that is as precise as
to deal with all kinds of transit, limb-darkening law, and
eccentricity. For this purpose, we use the approach introduced by Kopal (1977), who
reformulated the problem of evaluating the fractional loss of light as a cross-correlation
of two apertures; one representing the star undergoing eclipse, and the other the eclipsing
disc (in our case, the transiting planet). This method links the fractional loss of light
to the diffraction patterns of the two apertures as described in physical optics, allowing
the use of well known mathematical tools, like diffraction integrals. With this approach,
it is found that the
functions can be expressed, using a Hankel transform, as
After some more or less complicated algebra, Eq. (10) can be
expressed, without any loss of generality, in terms of Jacobi polynomials of the type
Gn(p,q;x) with the result
Equation (11) is the most general expansion for the associated functions in algebraic form, valid for any degree n of limb darkening, eccentricity, and type of eclipse, whether total, annular, a transit, or an occultation. The right-hand-side of the equation is a convergent series of terms and the precision that can be achieved, under the adopted assumptions, is more than enough for the needs of existing or planned photometric surveys. This precision is, of course, obtained as a function of the degree adopted for the law of limb darkening, as well as that for the series of terms in the Jacobi polynomials. We have found no practical problem evaluating for any high value of j, even well above 1000, but in most cases a good geometric precision was obtained, better than , with values of j up to only 20, or below , with values of . The application of Eq. (11) in the analysis of light curves is not only simple but also very fast. The computer code needed for the evaluation of the total luminosity at any orbital phase, for a given set of parameters, requires only 25 lines of FORTRAN making use of available mathematical libraries and, using a standard laptop, predicted light curves could be computed faster than the results could be printed.
In order to check the validity of Eq. (11) for any type of transit and any orbital phase, we compared the predicted results with those using the set of equations by Mandel & Agol (2002). Adopting different types of limb-darkening, linear, quadratic or even no limb-darkening, a set of light curves were computed for realistic system elements, reproducing known cases of extra-solar planets with both sets of equations. They were found to agree within their internal precisions to any degree in luminosity, as pointed out earlier, depending on the adopted maximum value of j. The largest differences, as could be expected, are found at phases of internal tangency.
Of course the precisions mentioned concern the integration of Eq. (1) with respect to the geometry of the problem, which is no longer a limitation. An incorrect treatment of the stellar properties (limb-darkening, surface inhomogeneities, or a non-spherical stellar shape) may nevertheless limit the accuracy of the obtained system elements as discussed in Sect. 4.
One problem is to get analytical equations for the study of planetary transits and another to solve the light curve for the system elements. As already shown, the relations between the elements of the transit and their observed characteristics are transcendental, and not capable of any type of analytic inversion. In order to find the best combination of elements reproducing observed light curves, five parameters are to be considered: the relative radius of the star , the ratio of radii , the orbital inclination i, and the limb-darkening coefficients ua and ub, in the case of a quadratic law.
The phase of the transit ingress ,
is preferred to the radius of the star as
an independent parameter since it is easily and directly measured on the light curve. On
the other hand,
is defined by the sum of the relative radii of the two components
as a function of the latitude over the star as given by the orbital inclination. At
and, according to Eq. (8),
The best combination of the five elements, , k, i, u+, and u-, should be obtained by optimization of the O-C (observed-computed) light curve as expressed by the minimum value of the parameter. There are many ways to fit the elements of a non-linear function (see e.g. Press et al. 1992), such as that given by , but two especially simple methods have been chosen in order to illustrate the use of the equations.
A parametric search in the space domain of the possible values of the first three
k, and i, was adopted for light curves with the precision typical
of ground-based observations (0.003 in magnitude), adopting values for the
limb-darkening coefficients as derived from model atmospheres. This direct method requires
boundary values for the search for the best combination of elements. The value of
can be easily constrained observationally. For k, an upper value
is given by the solution for no limb-darkening at the middle of the transit, i.e.
while a lower value in k is given by,
The second method can only be applied to more precise light curves (better than 0.001 in magnitude), like those available from space-based observations. In this case the limb-darkening coefficients can be fitted and good results are obtained with a least-squares differential correction method for three independent parameters, k, i and u+, and a parametric search for the two other, and u-. For the differential corrections a linearized Taylor expansion is taken around approximate values of the parameters, xn, for each observation, i, and least squares are used to optimize the fit towards the true values by means of:
A good example of a light curve for an extra-solar planet obtained with a ground-based telescope is that of OGLE-TR-113 (Udalski et al. 2002). The quality of the light curve is easily evaluated with the observations outside of the transit, during the nearby constant phases, showing a standard deviation of in luminosity.
Under the assumption of a circular orbit, no third light, and a dark planet relative to the star, as supported by the observational evidence, the corresponding value of for the O-C light curve can be calculated using Eq. (11) and estimated elements. The photometric precision does not allow for an independent determination of the limb darkening coefficients of the star, and they were adopted for the quadratic law (3) from Claret (2000) using the ATLAS model atmospheres for an effective temperature of K, as given by Konacki et al. (2004), in the photometric I band. The adopted values were and . It is important to note that, even though the quality of the light curve may not be good enough to distinguish between a linear and a quadratic limb-darkening law, the adoption of the second-order equation is needed to avoid systematic effects in the fitted parameters.
In order to proceed with the search for the best values of , k, and i in the parameter space, realistic ranges were chosen consistent with the shape of the light curve as described earlier in this section. Values for were used within 0.0275-0.0295, while k was kept within 0.129-0.158, and i between 80 and 90 degrees. A full analysis of the grid was then performed with an adequate resolution (0.0001, 0.001, and 0.1 for , k and i, respectively). The best elements were found to be , and , i.e. . The uncertainties are indicative of the internal errors, as well as those in the adopted value of the limb-darkening coefficients. Clearly, the orbital inclination could not be well determined, as was also found to be the case by Konacki et al. (2004) or Bouchy et al. (2004).
The absolute radii of the star and the planet can now be obtained using the size of the orbit as given by the radial velocity measurements. Both Konacki et al. (2004) and Bouchy et al. (2004) provide similar results and an average value of the semi-amplitude of the variations, m/s, was adopted. With this value and an estimated mass for the star, also derived from the spectroscopic information, the radius of the orbit was obtained and thus the absolute radii: in solar units, for the star, and in units of the radius of Jupiter, for the planet.
The light curve of the transit of the planet in HD 209458 was accurately measured using HST observations by Brown et al. (2001). The quality of the light curve, as indicated by the standard deviation of the measurements outside of the transit ( in luminosity) is clearly much better than in the previous case but it should be remembered that a slightly lower quality was achieved for the observations in one of the four series of measurements, the first one. This is visible in the dispersion of individual points for that orbit. In any case, the precision of the light curve is some two orders of magnitude better than the depth of the transit, and the limb-darkening coefficients can therefore be obtained while the differential corrections approach becomes applicable.
As explained above, a simultaneous search for the minimum , using and u- in the parameter space and optimizing the remaining k, i, and u+ with the differential corrections Eq. (14) was adopted. A linear limb-darkening solution was attempted first and a coefficient was obtained. Nevertheless, a much better fit to the observational data could be achieved with a quadratic limb-darkening law. In this case, the resulting parameters were , , and . The obtained limb-darkening coefficients were and in good agreement with model atmospheres. The corresponding relative radii of the star and the planet are thus determined to better than to be 0.1136 and 0.0137, respectively.
Taking the orbital parameters given by Mazeh et al. (2000) into account, the absolute radii of the star and the planet are found to be solar radii and Jupiter's radii, respectively. These values are quite close to those published by Brown et al. (2001) but with the precision essentially dominated by the uncertainty in the size of the orbit rather than the geometry of the system. Results also show excellent agreement with the most recent re-analysis by Winn et al. (2005) or Wittenmyer et al. (2005).
It is interesting to note that a simple, but not fully calibrated, photometric estimation by Giménez (2000) already indicated an effective temperature for the star of 6060 K with a metal content of , both agreeing with values given by the detailed spectroscopic analysis by Mazeh et al. (2000). On the other hand, a calibrated infrared method by Ribas et al. (2003), much more appropriate for this kind of study, not only provided an effective temperature of K but also the absolute radius of the star as , in excellent agreement with the value obtained from the light curve.
Analytical equations for the study of extra-solar planetary transits are given, together with simple approaches to derive accurate geometrical elements. Two representative cases have been analyzed with these equations leading to improved results. FORTRAN subroutines for the calculation of the relevant equations are available from the author on request. Of course, the method and equations described in this paper can also be applied to eclipsing binaries with small components so that the same tools can be used for the study of all possible photometric candidates derived from large surveys, independent of their nature as extra-solar planets or not. The distinction between a planet and a very low-mass star will only be possible with the radial velocity curve.
In this paper, the adopted equations are valid for any type of eclipse and degree of limb darkening but assume spherical components and homogeneous distribution of light over the stellar surface. Let us discuss these two assumptions further. Star spots are the usual source of surface inhomogeneities in solar-type stars, and their presence increases with stellar activity driven by fast rotation. Fortunately, spots can be detected through the photometric behavior of the star outside of the transit and then disentangled from the geometric properties. Moreover, in the event of the planet transiting over a spotted region, the shape of the light curve would be modified in a distinctive manner (Hui & Seager 2002). The tools presented in this paper allow for the geometrical aspects of the problem not to be a limitation in the understanding of the observations.
On the other hand, the spherical shape of the star may also be questioned as soon as stellar rotation and mutual tides are considered. The elongation effect of tidal interactions between the star and the planet should be a function of phase and therefore detectable, if at all significant (given the small mass ratio of the planet with respect to the star). In addition, during the transit, the position of the elongation with respect to the plane of the sky makes the projected shape of the star to be closest to circular. In the case of the planet, where the effect of tides may be significant, the observed relative radius can be easily corrected to obtain the equivalent to a sphere value using the ratio of masses and size of the orbit. The situation is different for the rotational effects. There is no phase dependence and therefore no photometric footprint. Fortunately the effect is expected to be so small that it does not make any difference in the analysis of the light curve. Most stars hosting planets are slow rotators, like our Sun. This may be due (see e.g. Barnes 2001) to an observational bias (so that spectral lines used for the measurement of very small radial velocities are not rotationally broadened) or to real evolutionary aspects (age and angular momentum transfer). The fact is that with small relative radii and rotational velocities, no deviation from the spherical shape can be detected with the level of precision of the available light curves. In the two cases of planetary transits analyzed in this paper, the rotational velocities of the parent stars are well known, and they correspond to a deviation from the spherical shape at the equator of the order of 10-5 in units of the relative radius. This is negligible when compared to the observational uncertainties of the derived parameters by more than 2 orders of magnitude.