Free Access
Volume 598, February 2017
Article Number L3
Number of page(s) 5
Section Letters
Published online 26 January 2017

© ESO, 2017

1. Introduction

Transmission spectroscopy is an essential tool for characterising the atmospheres of transiting exoplanets (see e.g. Charbonneau et al. 2002; Pont et al. 2013; Madhusudhan et al. 2014; Sing et al. 2016, and references therein). Snellen (2004) demonstrated that narrowband exoplanet features (e.g. sodium) could be probed by analysing the shape of the stellar absorption lines as a planet occults its host star, that is, by studying the chromatic Rossiter-McLaughlin (RM) effect; more recently Di Gloria et al. (2015) have shown that this effect can also be used to probe broadband signatures (e.g. Rayleigh scattering). Recently, Borsa et al. (2016, hereafter B16) presented a new technique using (a version of) line-profile tomography, with the intent of studying chromatic changes in planetary radii. However, we demonstrate in this Letter that the technique in B16 is unfortunately flawed.

In principle, the application of line-profile tomography is well motivated for exoplanet atmosphere characterisation if correctly implemented. This technique isolates the starlight behind the planet during transit (Collier Cameron et al. 2010), and the ratio of the integrated flux within the local profile (behind the planet) to the out-of-transit profile is equal to the brightness of the occulted starlight. In turn, the ratio of the occulted starlight is dependent on the planet-to-star radius ratio. As such, if one had space-based spectra then the planet radius could be recovered from the spectra alone, and doing so in various passbands would characterise the planetary radius wavelength dependency. For ground-based spectra this is not possible due to, for example, transparency variations of the Earth’s atmosphere. To remove these effects, ground-based spectra must first be continuum normalised; they can then be multiplied by a transit light curve to allow one to study the local, occulted stellar profiles directly by subtracting the in- from out-of-transit observations (see, e.g. Cegla et al. 2016a). Since the transit light curve is dependent on the planet-to-star radius ratio, it sets the planet radius that is recoverable when examining the brightness ratio between the local and out-of-transit profiles (meaning one cannot recover new radius variations following B16).

In this study, we first break down the physical implications of the B16 technique, and then simulate the planet transit of HD 189733 b to illustrate the impact of the technique’s shortcomings on the measured planet radius; we also reapply the B16 technique to HARPS data, with a more rigorous error propagation. In Sect. 2, we demonstrate how the choice of the transit light curve normalisation sets the recoverable planet radius. We present the simulated planet transits, and our reanalysis of the HARPS data in Sect. 3, and show how an underestimation of the errors can lead to spurious claims of planetary radius variations. Finally, we summarise our conclusions in Sect. 4.

2. Limitations of a brightness-based approach to chromatic exoplanet radii

2.1. The Borsa et al. 2016 approach

In B16, the authors attempted to study passband-dependent planet radius variations by averaging together the cross-correlation functions (CCFs) for subsets of HARPS spectral orders. All of these CCFs, regardless of passband, were continuum normalised and then further scaled using a Mandel & Agol (2002) transit light curve based on the system parameters determined in the full HARPS passband, except for passband dependent limb darkening coefficients.

For each passband, they created master out-of-transit CCFs by averaging together all the individual out-of-transit CCFs in a given night, and then subtracting the out- from the in-transit data to obtain CCFs for the starlight occulted by the planet. The authors then fitted Gaussian functions to the respective CCFs, presumably to act as a proxy for the integrated flux within, and argued that the ratio between the areas (determined from the Gaussian fits) of the local to the out-of-transit CCFs could serve as a measure of the brightness of the missing starlight, β, occulted during the in-transit observations.

The authors then compared this empirical β to the expected brightness ratio based on an approximation to the solution for integrating the starlight behind the planet given the particular geometry of the system (and assuming a particular function for the limb darkening). Since this solution was dependent on the planet radius, the authors argued that they were able to disentangle a planet radius measurement for each passband.

We demonstrate in Sect. 2.2 that the technique outlined by B16 can only recover the planetary radius set by the transit light curve used in the initial normalisation. However, we emphasise that such a limitation is set due to a flux-based approach and would not be present in an RV-based approach, such as that in Di Gloria et al. (2015) or that in traditional line profile tomography (i.e. following the Collier Cameron et al. 2010, formulation).

2.2. Transit light curve normalisation

B16 started by determining the area of their out-of-transit master CCFs, Aout, for each passband. Since the out-of-transit master CCFs (CCFout) were continuum normalised their area is equal to their equivalent width, EWout, that is, (1)The authors also measured the areas of the local CCFs behind the planet, Aloc. Since the in-transit CCFs (CCFin) are normalised by the transit light curve, the area of the local CCF (CCFloc) is (2)where flc is the flux from the light curve, EWloc is the equivalent width of the local CCF, and βlc is the fraction of starlight occulted by the planet under the assumptions of the transit light curve used in the normalisation.

In an ideal case, where the local stellar photospheric profiles can be represented by constant Gaussian functions (assumed both above and in B16), then the only difference between the local and disc-integrated out-of-transit CCFs is the broadening by stellar rotation present in the disc-integrated observations. Since the rotational broadening preserves the equivalent width, then (3)where the EWout and the EWloc are not only equal to each other but are also equal to the in-transit equivalent width, EWin (since it is simply a summation of local profiles of equal EW), and the ratio of the areas becomes: (4)Hence, under the above assumptions one can only recover the planet radius injected into the model transit light curve, regardless of which passband is studied. This is because the transit light curve normalisation effectively sets the area of the local profile. We note, that the transit light curve normalisation for the above follows Cegla et al. (2016a, hereafter C16), where CCFloc = CCFoutflcCCFin; whereas, the transit light curve normalisation in B16 was (5)(Borsa, priv. comm.), such that (6)where is the local CCF obtained following C16. Since the continuum is a free parameter in the Gaussian fits, the ratio of the areas is still equal to that in Eq. (4).

As such, the technique implemented by B16 represents a circular argument. Moreover, if one had the broadband photometry necessary for the correct spectral normalisation, then the planet radii could be determined directly from the light curves alone.

We stress that a solely brightness-based approach to transmission spectroscopy uses only the equivalent width, and therefore precludes any retrieval of information on exoplanet radii. It is the inclusion of the spectral dimension that is necessary to determine Rp (i.e. utilising the Doppler information, as is done in Collier Cameron et al. 2010 and Di Gloria et al. 2015).

3. Systematic effects on chromatic radius measurements

To try to understand how B16 obtained results mimicking Rayleigh scattering, we simulated the transit of HD 189733 b and applied their technique. To demonstrate that we could recover our model inputs, we also present results wherein we applied the normalisation from C16 with the correct planet radii for each passband; in doing so we discovered the approximations for β in B16 underestimated the planet radii, and thus we also explored a numerical approach for calculating this brightness ratio that was more accurate than the approximation used in B16.

3.1. Analytical brightness ratio approximation

Unfortunately, integrating the limb darkened brightness underneath a planet lying off stellar disc centre is not straightforward and even approximate analytical expressions are quite complex (especially if considering ingress and egress regions). For this reason, B16 used the β formalism presented in Collier Cameron et al. (2010), which was based on the analytical approximations of Ohta et al. (2005) for the RM effect. Therein the brightness ratio for a fully in-transit planet (i.e. no ingress or egress regions were considered) was defined, under the standard linear limb darkening law, as: (7)where R is the stellar radius, Rp is the planet radius, u1 is the linear limb darkening coefficient, and μ is the centre-to-limb planet position. We note that , where θ is the centre-to-limb angle, and xp,yp is the centre of the planet (see Collier Cameron et al. 2010; Cegla et al. 2016a, for details). We also note that Ohta et al. (2005) state that the accuracy of this approximation diminishes with increasing Rp/R, arguing that the additional terms in the analytical solution contribute to ~1% if Rp/R is 0.1 and up to few percent if Rp/R ~ 0.3. Given that Rp/R is only predicted to vary a few precent in wavelength for particular atmospheric characteristics, such as Rayleigh scattering, this approximation may inject systematic errors that could be misinterpreted as having a physical origin (even if the light curve normalisation is done correctly).

3.2. Numerical brightness ratio approximation

To investigate the impact of the accuracy of Eq. (7), we decided to calculate β numerically. For this approach, the in-transit starlight blocked by the planet is still defined as β = Floc/F, with F and Floc defined as the fluxes of the total stellar disc and the stellar disc under the planet, respectively. Note the observed brightness of the un-occulted star can be analytically determined by integrating a given limb darkening law over the projected stellar disc; for a linear limb darkening this is (8)where φ is the azimuthal angle. As previously stated, calculating the flux behind the planet analytically is not trivial. Hence, we calculated Floc numerically by constructing a square stellar grid with a width of 2Rp centred about the planet position (xp,yp), with n equal steps in the vertical and horizontal direction. Contributions from steps that did not lie beneath the planet and/or on the stellar disc were excluded. Thus, we approximated the flux behind the planet as (9)where Ixy is the limb darkened intensity at a given position in the aforementioned grid and (2Rp/n)2 is the corresponding area.

thumbnail Fig. 1

Recovered average planet radius, Rp, from the simulated data for each passband as a function of wavelength (plotted using the middle wavelength in the passband). Filled circles represent when the simulated Rp was constant, and hollow circles represent when the simulated Rp varied; in both cases these are shown in red with the error bars reported by B16 for comparison purposes only. Results from the B16 procedure are shown in black, while those from the numerical approximation herein and the C16 formulation are in blue. The lines represent linear fits to the simulated Rp shown for viewing ease.

Open with DEXTER

Our aim was to try to recover Rp as a function of wavelength, wherein the injected Rp was only used to construct the correct light curves (acting as if we had simultaneous multi-colour photometry). Hence, when trying to recover Rp, we started with the broadband planet radius and then allowed it to vary by up to ±0.005 R in steps of 0.0001 R. The recovered planet radius then corresponded to the planet radii that minimised the difference between β and Aloc/Aout.

thumbnail Fig. 2

Recovered average planet radius, Rp, for each passband as a function of wavelength for the observed HARPS data. Results from B16 are shown in green (where the wavelength error bars represent the passband wavelength region), and those from our reanalysis are shown in red when using the oversampled CCF and in black when using only every one in four points in the CCF. Subplot: the nightly recovered Rp when using only every one in four points in the CCF (night indicated by colour). The horizontal dashed lines show the mean Rp recovered by the B16 method on the simulated data (i.e. the solid black points in Fig. 1).

Open with DEXTER

3.3. Simulated star-planet system

We used the simulated stellar grid of Cegla et al. (2015, 2016b) and injected into each grid cell a Gaussian profile with a full-width half-maximum (FWHM) of 5 km s-1 (note this width is similar to the expected value for the stellar photosphere). In the simulated star we did not consider any astrophysical effects (i.e. granulation or starspots etc.) other than rigid body stellar rotation, which was set to the value obtained by C16, 3.25 km s-1. We also assumed an edge-on (i = 90°) aligned orbit.

The transit was sampled in 21 equal steps in phase from −0.02−0.02, centred about mid-transit, with an additional sample at phase = 0.03 to serve as a completely out-of-transit reference. We simulated a transit for each of the seven passbands (from 400−700 nm) used in B16, and applied a linear limb darkening using the coefficients (for each passband) these authors provided. For each of the seven passband transits we injected a planet with a constant radius equal to the value assumed by B16 for the whole HARPS passband (Rp = 0.1581 R, hereafter referred to as the broadband Rp), but varied the limb darkening accordingly.

For this set of transits, we tested the impact of the transit light curve normalisation. In the first case, we followed the procedure in B16, and in the second case we normalised the data following C16 and used the numerical approximation in Sect. 3.2 to estimate Rp. Examining the first case allowed us to examine any errors introduced using the β approximation and/or the B16 normalisation. On the other hand, the second case offered a test case to ensure we could recover the model inputs.

For a second test, we repeated the above, but varied the radius of the simulated planet; for this we selected Rp equal to the values reported by B16 for each passband. Again we tested two cases: first following the B16 procedure (where the light curve limb darkening varies in each passband, but the light curve radius remains fixed at the broadband Rp), and the second using the C16 normalisation (we note the assumed light curve has the correct Rp here) and the numerical approximation in Sect. 3.2.

We then examined the recovered Rp as a function of stellar disc position, and found only a slight dependance on disc position when following B16. However, if one point each at the ingress and egress regions were included then the dependence on disc position was strong, and including such data would systematically decrease the recovered Rp (as the β formulation is not valid in these regions).

Moreover, we found that, regardless of the B16 or C16 normalisation, using the β approximation always underestimated the limb darkening behind the planet and therefore also underestimated the true planet radius. This is because the analytical approximation assumes the limb darkening behind the planet is constant, and equal to the value behind the centre of the planet. In reality, the stellar photosphere behind the planet exhibits a range of limb darkening. This is why the numerical model in Sect. 3.2 is necessary to recover the Rp injected into the simulated data.

The Rp reported in B16 comes from averaging together the planet radii recovered across the stellar disc. If the limb darkening effects are sufficiently removed (and the stellar profile is constant), then this provides a good means to boost the signal-to-noise in the reported Rp. In Fig. 1, we present the average recovered Rp as a function of wavelength from the simulations, for both tests (when Rp was constant and when it varied). As expected from Sect. 2.2, the B16 procedure always results in nearly the same Rp, regardless of whether the true Rp varied or not.

For our numerical approach and the C16 normalisation, we demonstrate accurate recovery of Rp (regardless of whether or not we include ingress and egress data), but only if the light curve normalisation is done with the correct Rp for each passband (using the broadband Rp for all passbands meant only the broadband Rp was recovered). Hence, regardless of the normalisation (i.e. B16 or C16) or the brightness formulation (i.e. β or our numerical approximation), we could only retrieve the parameters injected into the system via the transit light curve normalisation, as expected from Sect. 2.

3.4. Reanalysis of the HARPS data

Our application of the B16 procedure on the simulated data cannot explain the wavelength-dependent planet radii reported in B16. To further investigate this aspect, and to ensure we have applied the B16 method correctly, we have reanalysed the same three transits of HD 189733 b following their technique, but using the Levenberg-Marquardt least-squares minimisation from MPFIT (Markwardt 2009, and references therein) rather than IDL’s GAUSSFIT1. The results are plotted in Fig. 2, alongside those from B162. We demonstrate we can reproduce (red points in Fig. 2) results in 1–2σ agreement with B16 (in green); hence, we are confident we have applied their technique properly (in both the simulated and observed data). However, we argue that with the correct treatment of the uncertainties the recovered planet radii (in black) are consistent with a flat line (within 1−3σ of the mean Rp recovered in the simulated data, i.e. the solid black points in Fig. 1), as expected from Sect. 2. We believe the reason B16 report a trend with wavelength, and we do not, is largely due to differences in our error analysis and Gaussian fitting techniques.

In B16, the recovered Rp for each stellar disc position and all three transits were averaged together to provide one Rp for each passband, and the reported errors came from the rms of these individual planet radii (i.e. the standard deviation divided by the square root of the total number in the passband). In our analysis, we report the weighted mean for each passband, with the weights being the inverse square of the error for each individual planet radii (where the error was calculated by propagating the errors on the CCF areas as reported from the Gaussian fits following C16, and assuming negligible error on the limb darkening and stellar disc positions). The error on the weighted mean then was simply the square root of the inverse sum of the weights squared. If the errors on individual Rp were all exactly equal to the standard deviation, then the two approaches would yield the same result.

In addition to this slight difference in error analysis, B16 also applied their Gaussian fits to the oversampled CCF grid provided by the HARPS pipeline (Borsa, priv. comm.). We caution against such an approach, as the oversampling will lead to a significant underestimation of the errors. Hence, we also fit Gaussians to data composed of every one in four points from the original CCFs (to compensate for the original sampling rate of 0.25 km s-1 for a mean pixel width of 0.82 km s-1); these results are shown in black in Fig. 2.

Table 1

Best fits to observed data.

To test the significance of a trend in Rp with wavelength, we fitted the data with both a flat line and a linear regression, and calculated the reduced chi-squared, , and the Bayesian information criterion (BIC); the results are shown in Table 1. We note that even if a wavelength-dependent Rp is found, it does not confirm the B16 technique is valid, as we have already shown it is not mathematically possible to retrieve radius variations. Rather, it would serve as a red flag that we do not fully characterise the interplay of the various complexities present in the observations. In particular, stellar activity can alter the observed stellar line shapes and their equivalent widths – which in turn could lead to spurious radius variations following Sect. 2. Since HD 189733 is a known active star this is likely scenario; and in agreement with Fig. 5 from B16, wherein the single night analysis with the most apparent slope, July 2007, is also the most magnetically active (Cegla et al. 2016a). Moreover (and as noted by B16), McCullough et al. (2014) have argued that the apparent wavelength dependency in their independent observations of this system are best explained by un-occulted starspots rather than the planet atmosphere.

When using the oversampled CCFs, both our analysis and B16’s indicate a slight improvement in fit for the model with a wavelength-dependent slope. However, we find a much worse fit to the data than that found with the B16 results. The high from our reanalysis indicates an underestimation of the uncertainty in the data, as one would expect when using the oversampled CCFs. We cannot explain the very low for the B16 wavelength-dependent fit, which indicates the model is overfitting the data. We note these tests were only performed on our reanalysis of the oversampled data for comparison with B16; for our conclusions on the best-fit, we refer the reader to the analysis on the CCFs sampled every one in four points.

For the properly sampled dataset, we found only a marginal improvement in the fit for the wavelength-dependent model, and do not deem this improvement to be statistically significant (see Table 1). Moreover, the best-fit flat model (Rp = 0.1569 ± 0.0003) lies within 3σ of the mean Rp predicted by the simulations. The slight improvement for the sloped model is also heavily influenced by only a couple data points from a single transit, in August 2007, as shown in the subplot of Fig. 2. If the best-fit model is robust, it should withstand removing the August transit; however, doing so means the data is then best-fit by a flat line (, BIC = 11.1 and , BIC = 12.6 for flat and sloped line, respectively). Consequently, we believe B16 likely report a wavelength-dependent trend in Rp due to insufficient error analysis, and that its agreement with the literature may be purely coincidental.

4. Conclusions

We outline our conclusions on a point-by-point basis below.

  • The technique presented in B16 using the ratio of the areas of the local (starlight behind the planet) to the out-of-transit CCF cannot be used to determine Rp, as Rp must be known a priori for the transit light curve normalisation required for ground-based spectra. This is shown both analytically and using a simulated star-planet system.

  • The analytical β approximation used in B16 also introduces (slight) systematic trends with planet position due to inadequately accounting for limb darkening and fractional area occultation effects, and underestimates the value of Rp.

  • We postulate that the Rp variations reported in B16 are likely due to underestimated errors (largely originating from use of oversampled CCFs), as our reanalysis of the HD 189733 b transits further demonstrates that the only Rp recoverable is that injected into the transit light curve normalisation.

  • Chromatic RM measurements from ground-based spectra are not possible without taking the Doppler information into account. Hence, for future measurements, we advise readers to either follow the works of Snellen (2004), Di Gloria et al. (2015) or to apply the line-profile tomography of Collier Cameron et al. (2010) directly on each spectral passband.


MPFIT did not produce significantly different results compared to IDL’s GAUSSFIT, but it did allow us to propagate our errors more thoroughly (see C16 for details).


We note that we used the same transit parameters as B16, but there is a slight difference in the template mask used to obtain the CCFs. B16 used the archival data available from the ESO website, where 2 nights used the G2 mask and 1 night used the K5 mask, whereas our data always used the K5 mask. However, this difference is unlikely to impact the analysis since each night had its own master CCFout.


We thank the referees, I. A. G. Snellen and S. Albrecht, for their careful reading and constructive comments, which improved the clarity of the manuscript. We also thank F. Borsa for useful discussions. Additionally, H.M.C. thanks E. de Mooij for suggesting we use simulated stars to test the work herein. H.M.C. and C.A.W. gratefully acknowledge support from the Leverhulme Trust (grant RPG-249). H.M.C., V.B., C.L. and A.W. acknowledge the financial support of the National Centre for Competence in Research “PlanetS” supported by the Swiss National Science Foundation (SNSF). C.A.W. also acknowledges support from STFC grant ST/L000709/1, and A.W. acknowledges additional financial support directly from the SNSF. This research has made use of NASA’s Astrophysics Data System Bibliographic Services.


  1. Borsa, F., Rainer, M., & Poretti, E. 2016, A&A, 590, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  2. Cegla, H. M., Watson, C. A., Shelyag, S., & Mathioudakis, M. 2015, in 18th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, eds. G. T. van Belle, & H. C. Harris, 567 [Google Scholar]
  3. Cegla, H. M., Lovis, C., Bourrier, V., et al. 2016a, A&A, 588, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Cegla, H. M., Oshagh, M., Watson, C. A., et al. 2016b, ApJ, 819, 67 [NASA ADS] [CrossRef] [Google Scholar]
  5. Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L. 2002, ApJ, 568, 377 [NASA ADS] [CrossRef] [Google Scholar]
  6. Collier Cameron, A., Bruce, V. A., Miller, G. R. M., Triaud, A. H. M. J., & Queloz, D. 2010, MNRAS, 403, 151 [NASA ADS] [CrossRef] [Google Scholar]
  7. Di Gloria, E., Snellen, I. A. G., & Albrecht, S. 2015, A&A, 580, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Madhusudhan, N., Knutson, H., Fortney, J. J., & Barman, T. 2014, Protostars and Planets VI, 739 [Google Scholar]
  9. Mandel, K., & Agol, E. 2002, ApJ, 580, L171 [NASA ADS] [CrossRef] [Google Scholar]
  10. Markwardt, C. B. 2009, in Astronomical Data Analysis Software and Systems XVIII, eds. D. A. Bohlender, D. Durand, & P. Dowler, ASP Conf. Ser., 411, 251 [Google Scholar]
  11. McCullough, P. R., Crouzet, N., Deming, D., & Madhusudhan, N. 2014, ApJ, 791, 55 [NASA ADS] [CrossRef] [Google Scholar]
  12. Ohta, Y., Taruya, A., & Suto, Y. 2005, ApJ, 622, 1118 [NASA ADS] [CrossRef] [Google Scholar]
  13. Pont, F., Sing, D. K., Gibson, N. P., et al. 2013, MNRAS, 432, 2917 [NASA ADS] [CrossRef] [Google Scholar]
  14. Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59 [NASA ADS] [CrossRef] [Google Scholar]
  15. Snellen, I. A. G. 2004, MNRAS, 353, L1 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Best fits to observed data.

All Figures

thumbnail Fig. 1

Recovered average planet radius, Rp, from the simulated data for each passband as a function of wavelength (plotted using the middle wavelength in the passband). Filled circles represent when the simulated Rp was constant, and hollow circles represent when the simulated Rp varied; in both cases these are shown in red with the error bars reported by B16 for comparison purposes only. Results from the B16 procedure are shown in black, while those from the numerical approximation herein and the C16 formulation are in blue. The lines represent linear fits to the simulated Rp shown for viewing ease.

Open with DEXTER
In the text
thumbnail Fig. 2

Recovered average planet radius, Rp, for each passband as a function of wavelength for the observed HARPS data. Results from B16 are shown in green (where the wavelength error bars represent the passband wavelength region), and those from our reanalysis are shown in red when using the oversampled CCF and in black when using only every one in four points in the CCF. Subplot: the nightly recovered Rp when using only every one in four points in the CCF (night indicated by colour). The horizontal dashed lines show the mean Rp recovered by the B16 method on the simulated data (i.e. the solid black points in Fig. 1).

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.