A SEA BASS on the exoplanet HD209458b
Giuseppe Morello
astro-ph.EP
















Irfu, CEA, Université Paris-Saclay, 91191 Gif-sur-Yvette, France
Department of Physics & Astronomy, University College London, Gower Street, WC1E 6BT, UK
giuseppe.morello@cea.fr














We present here the first application of Stellar and Exoplanetary Atmospheres Bayesian Analysis Simultaneous Spectroscopy (SEA BASS) on real datasets.  SEA BASS is a scheme that enables the simultaneous derivation of four-coefficient stellar limb-darkening profiles, transit depths, and orbital parameters from exoplanetary transits at multiple wavelengths. It relies on the wavelength-independence of the system geometry and on the reduced limb-darkening effect in the infrared. This approach has been introduced by Morello et al. (2017) (without the SEA BASS acronym), who discuss several tests on synthetic datasets.
Here, we (1) improve on the original algorithm by using multiple Spitzer/InfraRed Array Camera passbands and a more effective set of geometric parameters, (2) demonstrate its ability with Hubble Space Telescope/Space Telescope Imaging Spectrograph datasets, by (3) measuring the HD209458 stellar limb-darkening profile over multiple passbands in the 290–570 nm range with sufficient precision to rule out some theoretical models that have been adopted previously in theliterature, and (4) simultaneously extracting the transmission spectrum of the exoplanet atmosphere. The higher photometric precision of the next-generation instruments, such as those onboard the James Webb Space Telescope, will enable modeling the star-planet systems with unprecedented detail, and increase the importance of SEA BASS for avoiding the potential biases introduced by inaccurate stellar limb-darkening models.




































§ INTRODUCTION

Over the past 15 years, transit spectroscopy has enabled the detection of atomic, ionic, and molecular species in the atmospheres of exoplanets (e.g., ), as well as the characterization of their thermal and cloud profiles (e.g., ). The diagnostic parameter is the squared planet-to-star radii ratio, p^2 = (R_p/R_*)^2, the so-called “transit depth”. The optically thick area of a transiting exoplanet is wavelength-dependent due to the differential opacity of its gaseous envelope. The expected spectral modulations are of the order of 100 parts per million (ppm) relative to the stellar flux for the hot Jupiters <cit.>, and smaller for other classes of planets.



The shape of a transit light-curve also depends on the stellar limb-darkening effect, i.e., the wavelength-dependent radial decrease in specific intensity, coupled with the star-planet geometry and timing. Accurate knowledge of the stellar intensity profile is necessary to infer the correct transit depth from observations, especially at the precision level required to study the exoplanetary atmospheres. In most cases, a set of limb-darkening coefficients is computed either from specific stellar-atmosphere models or by interpolation on a precalculated grid. In some cases, for the highest quality datasets, empirical limb-darkening coefficients are estimated from the light-curve fits by assuming a linear <cit.> or quadratic <cit.> law:

    I_λ( μ ) / I_λ( 1 )  = 1 - u (1-μ) 
     I_λ( μ ) / I_λ( 1 )  = 1 - u_1 (1-μ) - u_2 (1-μ)^2

Here μ=cosθ, where θ is the angle between the line-of-sight and the normal to the stellar surface. Several authors have found that the linear and quadratic limb-darkening laws can be inadequate to approximate the stellar intensity profile. <cit.> showed that such inadequate parameterizations can provide a perfect match to the observed transit light-curves while significantly biasing the transit depth and other parameter estimates. <cit.> proposed a four-coefficient law, hereinafter “claret-4”, which enables a more accurate approximation to the stellar model-atmosphere intensity profiles:

    I_λ( μ ) / I_λ( 1 )  = 1 - ∑_n=1^4 a_n  ( 1 - μ^n/2 )

The claret-4 coefficients are not commonly fitted to the transit light-curves, as the strong parameter degeneracies are likely to hamper the convergence of the fit.

<cit.> suggest an optimal strategy to break some degeneracies by adopting informative priors on the orbital parameters, hence enabling the convergence of model-fits with free claret-4 coefficients. Their proposal was to measure the “geometric” parameters, a/R_* and i, from infrared transit observations, then implement the results as informative priors when model-fitting at shorter wavelengths. This approach is motivated by the smaller limb-darkening effect in the infrared, which mitigates the potential biases that are due to inaccurate models, and the negligible wavelength-dependence of a/R_* and i. Its effectiveness has been demonstrated over a set of simulated observations.

Here, we apply this technique (now named Stellar and Exoplanetary Atmospheres Bayesian Analysis Simultaneous Spectroscopy (SEA BASS)), for the first time to real datasets, taken with the InfraRed Array Camera (IRAC) onboard Spitzer and the Space Telescope Imaging Spectrograph (STIS) onboard the Hubble Space Telescope (HST).
We focus on the HD209458 system, which consists of a hot Jupiter orbiting a solar analog star. Stellar and planetary parameters available in the literature are summarized in Table <ref>. HD209458b was the first transiting exoplanet to be discovered <cit.>, and also the first one on which a chemical element has been detected (Na, ). HD209458 is currently one of the most frequently studied exoplanetary systems, thanks to the relatively large amplitude of the transit signal and the apparent brightness (V=7.6) of its parent star.






 §.§ Structure of the paper

Section <ref> reports the observations analyzed in this paper. Section <ref> describes the state-of-the-art data-detrending pipelines adopted for the different instruments, which are not novel to this work, except for the “test of aperture sizes” (Section <ref>). Sections <ref> and <ref> report the results obtained for the infrared observations; more technical details are discussed in Appendices <ref>–<ref>. Section <ref> compares our parameter results over multiple epochs and passbands, and with previous results reported in the literature. Section <ref> discusses different ways of combining the geometric parameter estimates into a  Bayesian prior to adopt in the analyses of visible light-curves. Section <ref> reports the results obtained for the visible observations. Appendix <ref> discusses the impact of using two-coefficient limb-darkening laws instead of the most accurate claret-4.





§ OBSERVATIONS





 §.§ Infrared photometry



We reanalyzed all the transits of HD209458b observed with Spitzer/IRAC across the four passbands centered at 3.6, 4.5, 5.8, and 8.0 μm, i.e., the IRAC channels 1–4, respectively. Observational and detector information for the individual light-curves and their identifiers is given in Table <ref>.

In full-array mode, frames of 256×256 pixels are recorded with a cadence of 8.4 s (frame time + delay time).  Except for the 2007 December 24 observation, a detector switch occurred after each full frame, in order to sample the same transit event almost simultaneously at 3.6 and 5.8 μm or 4.5 and 8.0 μm. The telescope needed frequent repointings to keep the centroid within one pixel. This observing strategy is not available for the warm Spitzer mission, as the two longest wavelength channels are not operative. The adjectives “cold” and “warm” associated with Spitzer refer to the datasets taken before and after the cryogen exhaustion on 2009 May 15, respectively.

In sub-array mode, 64 frames from an area of 32×32 pixels are taken in succession. The delay time between consecutive groups of 64 frames is 2 s. No repointings are needed in this mode to keep the centroid within one pixel over a few hours.




 §.§ Visible spectroscopy


We also reanalyzed two transits of HD209458b that were observed spectroscopically with the HST/STIS G430L grating (wavelength range 290–570 nm) as part of the GO-9447 program. Each visit consists of five consecutive HST orbits (duration ∼95 minutes), including 50-minute gaps due to Earth occultations. Frames are 64×1024 detector images with an integration time of 22 s and a reset time of 20 s. The dispersion direction is almost parallel to the longest side of the detector, and it is stable within the same pixel column during each visit. Observational and detector information for the individual light-curves and their identifiers is given in Table <ref>.






§ DATA ANALYSIS



In this section we describe the specific data-detrending pipelines adopted for each instrument, from the extraction of raw light-curves to the final parameter fitting.

The model-fitting procedure consists of two steps, which are common to all of our data analyses:


  * least-squares minimization between the reference light-curve and the parametric model;

  * Markov Chain Monte Carlo (MCMC) sampling of the parameters posterior distribution.

Least-squares parameter estimates are used as starting parameters for the chains. Then we run an Adaptive Metropolis algorithm with delayed rejection <cit.> for 1,000,000 iterations, with a burn-in of 100,000 iterations. 
The transit parameters that are not free in the fit are fixed to the values reported in Table <ref>. For the Spitzer/IRAC light-curves, we adopt the passband-integrated limb-darkening coefficients provided by <cit.>, here reported in Table <ref>.






 §.§ Detrending Spitzer/IRAC data

We downloaded the Basic Calibrated Data (BCD) from the Spitzer Heritage Archive <cit.>. BCD are flat-fielded and flux-calibrated frames <cit.>. Each observation consists of a time series of frames taken with regular cadences (see Section <ref> for timing details).

We take square arrays of pixels with the stellar centroid at their centers as fixed photometric apertures during each observation, following the same procedure as adopted by, e.g., <cit.>, and <cit.>. The centroid position varies within the same pixel during each observation.

For the full-phase and half-phase observations (see Table <ref>), we restrict the data analysis to an interval |ϕ | < 0.05, where ϕ denotes the time from mid-transit in units of the orbital period, commonly referred to as orbital phase. This cut allows an out-of-transit baseline longer than 6 hr, ≲1.5 times longer than those associated with transit-only observations.

We flag and correct outliers in the flux time series. Outliers are those points that deviate from the median of the four nearest temporal neighbors by more than five times the noise level, estimated as the median standard deviation of the sets of five consecutive points. Outliers are replaced with the mean value of the points immediately before and after, or, in the case of two consecutive outliers, with a linear interpolation. In order to identify more than two consecutive outliers, Ø_i (i=1 … n), we considered the light-curves binned in time (a binning factor of up to 128 for the sub-array mode, and with different starting points), such that a single outlier, O, in a binned time series corresponds to the arithmetic mean of the consecutive outliers in the unbinned version. O is then replaced with the mean value of the points immediately before and after, (A+B)/2 (in the binned version). This operation is equivalent to shifting O by (A+B)/2 - O. As we wish to correct the
unbinned light-curve, all the consecutive outliers, O_i (i=1 … n), are then shifted by (A+B)/2 - O. The percentage of outliers is always smaller than 0.4%.




  §.§.§ IRAC channels 1 and 2

The two IRAC channels with the shortest wavelengths (3.6 and 4.5 μm) both use InSb detectors <cit.>. The most significant instrumental effects arise from fluctuations of the centroid position coupled with a non-uniform inter- and intra-pixel response <cit.>. The flux measurements are correlated to the centroid of the Point Spread Function (PSF) on the detector, which varies over multiple timescales <cit.>. In particular, the most prominent effect is a sawtooth modulation associated with the periodic on-off cycling of a battery heater within the spacecraft bus, although its amplitude depends on the solar angle during the observations.

We apply the pixel-Independent Component Analysis (ICA) technique to remove the instrumental systematics in the observations at 3.6 and 4.5 μm <cit.>. The key step is the ICA of the pixel time series within the chosen photometric aperture, i.e., a linear transformation of pixel time series set into maximally independent components  <cit.>. <cit.> proposed a wavelet pixel-ICA algorithm, in which the pixel time series undergo Discrete Wavelet Transform (DWT) before the ICA separation. The so-called wavelet pixel-ICA algorithm has proven to perform as well as or better than other state-of-the-art pipelines to detrend Spitzer/IRAC data <cit.>. The n-level DWT preliminarily separates the high-frequency components while down-sampling by a factor of 2^n. This trade-off was worthwhile in detrending observations with low signal-to-noise-ratio (S/N∼0.4[Here S/N is defined as the amplitude ratio between the transit/eclipse signal and the high-frequency noise.]) such as those analyzed by <cit.>. No preference between the two variants of the pixel-ICA algorithm has been found for the observations analyzed here (S/N∼3–5) hence we selected the simpler one (no wavelet).

One of the independent components has the morphology of a transit; the others are referred to as non-transit components, and may correspond to different source signals, either astrophysical or instrumental in nature, and noise. 
The raw light-curves are model-fitted as linear combinations of the non-transit components and a transit model, as in <cit.>. We set uniform priors for the MCMC. The output chains (approximately gaussian) are used to estimate the best parameter values (means) and lower limit error bars (standard deviations, σ_, 0). 

We also take into account the uncertainty associated with the independent components <cit.>:

    σ_ICA^2 = ∑_j o_j^2 ISR_j

where ISR is the so-called Interference-to-Signal-Ratio matrix <cit.>, and o_j are the coefficients of non-transit components. The ICA variance, σ_ICA^2, is added to the MCMC-sampled likelihood variance, σ_0^2, approximately equal to the variance of the residuals. The final parameter error bars are rescaled as

    σ_par = σ_par,0√(σ_0^2+σ_ICA^2/σ_0^2) .






  §.§.§ IRAC channels 3 and 4


The two longest wavelength IRAC channels (5.8 and 8.0 μm) use Si:As detectors <cit.>. The most significant instrumental effects are the so-called “detector ramps” <cit.>.

Assuming a charge-trapping model, <cit.> derived an exponential formula to model the increase in flux per pixel over time observed at 8.0 μm:

    F(t) = a_0 - a_1 e^-t/τ_1

They recommend the use of a double exponential formula to approximate the ramp on the photometric flux (e.g., integrated over an aperture):

    F(t) = a_0 - a_1 e^-t/τ_1 - a_2 e^-t/τ_2

The double exponential formula has the asymptotic behaviour predicted by the charge-trapping model, τ_1 and τ_2 are characteristic timescales associated with pixels with different illumination levels, and in the analyses of <cit.>, led to a smaller Bayesian Information Criterion (BIC, ).

The ramp observed at 5.8 μm has the opposite behavior in our datasets, i.e., the photometric flux decreases over time. While this behavior may not be explained by the same charge-trapping model <cit.>, Equations <ref> and <ref> can also be used to approximate the ramp at 5.8 μm.

A long list of alternative ramp parameterizations has appeared in the literature for both channels (e.g., ), all of which are linear combinations of three basic forms: exponential, logarithmic, and powers of time. We tested the different functions found in the literature, and compare their results.

Another frequently adopted procedure is to discard the initial part of the ramp (e.g., ), before any data detrending. We tested different lengths for the trimmed interval of each dataset: from 0 to 91 minutes, with a 3.5-minute step.

The raw light-curves are model-fitted as the product of transit and ramp models.




  §.§.§ Test of aperture sizes


In the study of exoplanetary transits, the relevant information is the fraction of stellar light that is occulted by the planet over time, while the absolute calibrated flux is not required. However, an additive constant on the flux measurements does change their ratios, therefore biasing the inferred stellar and planetary parameters (see, e.g.,  for the case of additive flux from the nightside exoplanet emission,  for the case of stellar spots).

Here we derive an analogous formula for the transit depth affected by background. 
Let us consider the simplest transit model, i.e., a dark circle (the planet) that passes in front of another circle (the star) with uniform brightness <cit.>. Before and after transit, the detector will record the unperturbed stellar flux: F_ = α F_*. Here F_* denotes the stellar flux, the α factor accounts for the fact that typically, the user does not calculate the absolute source flux, but rather a detector number that also depends on user choices (e.g., size of the photometic aperture, PSF fitting strategy). During the full transit (when the planet circle lies inside the stellar circle), the decrease in flux is proportional to the planet-to-star area ratio: Δ F = F_ - F_ = (R_p/R_*)^2 α F_*. Here F_ denotes the in-transit flux, and R_p and R_* the planetary and stellar radii, respectively. The transit depth is defined as follows:

    p^2 = Δ F/F_ =  ( R_p/R_* )^2

Now let us assume that the measured flux includes some background, B, constant over time. Note that B will affect the absolute flux values, but not the difference Δ F. Therefore, the observed transit depth becomes

    p^2(B) = Δ F/F_ = Δ F/α F_* + B = p^2/1+B/α F_*.

In other words, if the background is neglected, the inferred planet-to-star area ratio will be biased by a factor [1+B/(α F_*)]^-1. This relation holds in the more general case with stellar limb-darkening, although the more complex light-curve morphology requires a different definition of the transit depth parameter <cit.>.

Here we propose the following “test of aperture sizes” to measure the background to be subtracted from the data. The purpose of this test is not finding the best aperture to optimize the S/N on the result, but rather using the different apertures to measure the impact of background on the inferred transit depth, and correcting for it. We assume a uniform background per pixel, b, within a region that includes the source star. The background associated with a given aperture is proportional to the number of pixels enclosed, n, eventually including fractional pixels: B=nb.  By injecting this relation in Equation <ref>:

    p^2(n, b) = p^2/1+nb/α F_*

By measuring the transit depth using different aperture sizes, the background per pixel, b, can be estimated as the best-fit parameter to Equation <ref>.
Note that in the limit of low background (nb ≪α F_*) and large apertures (constant α), by using the first-order Taylor expansion for the denominator, Equation <ref> becomes linear in n:

    p^2(n, b) ≈ p^2  ( 1 - b/α F_* n  )


Traditional approaches to estimating the background consist of measuring the mean/median flux integrated on a ring around the star, or on a dark sky region (that does not contain the source). Possible sources of error with these methods are contamination from the tails of the source PSF or inhomogeneities between different sky regions. The test of aperture sizes minimizes these issues, as the background is measured from the variation in transit depth using multiple apertures containing the source. A negative background may be due to inaccurate bias/dark subtraction. The test automatically corrects for these potential offsets. The case of a background star underneath the source PSF is more subtle, as it undermines the assumption of locally uniform background in Equations <ref>–<ref>.



 §.§ Detrending HST/STIS data

We downloaded the flat-fielded science images (extension: flt.fits, ) from the Mikulski Archive for Space Telescopes (MAST, https://archive.stsci.edu). We computed white light-curves and spectral light-curves for five passbands uniformly spaced in wavelength as in <cit.>, allowing a direct comparison of the results. We tested rectangular apertures centered along the spectral trace with the cross-dispersion width ranging from 3 to 49 pixels. The available HST wavelength solution was adopted to delimit the apertures in the dispersion direction. 

The notable instrumental effects are short- and long-term ramps (see Figure <ref>). The short-term ramp is a flux variation correlated with the spacecraft orbital phase, which is highly repeatable between orbits within a visit, except for the first orbit, as a result of thermal settling of the telescope at its new pointing <cit.>. HST observations of exoplanetary transits are purposefully scheduled to have two pre-transit orbits, and discard the first one <cit.>. The long-term ramp approximates a linear trend in the mean flux value of each orbit. We adopt the  method to correct for these systematics <cit.>. The two in-transit orbits are divided by the time-weighted mean of the two out-of-transit orbits. We combine the observations into a unique series ordered by the exoplanet orbital phase, using the ephemerides reported in Table <ref>, but allowing independent phase shifts for the two visits. This approach, also followed by <cit.>, is useful for obtaining good sampling of the full transit. Finally, the -corrected, joint time series are MCMC-fitted using a transit model with transit depth, orbital parameters, phase shifts, and four limb-darkening coefficients as free parameters. Additionally, we tested different Bayesian priors on the orbital parameters (this is explained in Section <ref>), stellar limb-darkening coefficients reported in the literature, and different limb-darkening laws (see Appendix <ref>).





















§ INFRARED RESULTS




 §.§ Spitzer/IRAC channels 1 and 2 - 3.6 and 4.5 μm



Figure <ref> shows the raw light-curves, the relevant best-fit models, and residuals for the observations taken with the Spitzer/IRAC channels 1 and 2 (3.6 and 4.5 μm, respectively).
Figure <ref> shows the results of the test of aperture sizes (Section <ref>) using four arrays: 5×5, 7×7, 9×9, and 11×11 pixels. For the two cold Spitzer observations in channel 1 (left, top two panels) and the warm Spitzer observation in channel 2 (right, bottom panel; see also Table <ref>), the measured transit depths dependence on the aperture size is significantly below the uncertainty measurements. The best-fit values fall within intervals of 23 ppm (cold, channel 1) and 11 ppm (warm, channel 2). For comparison, the final error bars in transit depth for these observations are ∼140 ppm (cold, channel 1) and 44 ppm (warm, channel 2), i.e., they are at least four times larger than the relevant variations with the aperture size. The other three observations show clear trends in transit depth versus aperture size, even if the projected transit depths, with the background removed, are consistent with the uncorrected results within 1 σ. In two cases the transit depth increases with the aperture size, revealing a negative background, most likely due to imperfect dark subtraction.

Figure <ref> and Table <ref> report our final results, after background subtraction, for the transit depth and the orbital parameters. We selected the 5×5 array (both as photometric aperture and set of pixel light-curves for ICA), as it gives the lowest residuals and smallest uncertainties in the ICA components (Equation <ref>), in agreement with our previous studies <cit.>. 

The parameters estimated from warm Spitzer observations have 2–5 times smaller error bars than those obtained from cold Spitzer observations, thanks to an optimized observing strategy with higher efficiency and no repointings (see Section <ref>). The warm transit depths appear to be systematically smaller than the cold ones. The discrepancies between the different couples of measurements at the same wavelength are 2–5 σ, or 300–800 ppm between their central values. The other transit parameters also appear to have systematic differences between cold and warm observations, but they are consistent within 1 σ.

In Appendix <ref>, we show that when we inject appropriate offsets, i.e., additive constants, into the original light-curves in order to force the transit depths to assume identical values at all epochs, the corresponding best-fit orbital parameters also tend to be closer to each other. This behavior suggests that the discrepancy between cold and warm measurements is, most likely, caused by the imperfect correction of the different detector nonlinearities. If this is the case, the offset will be reduced with a careful reanalysis of the available data in collaboration with the Spitzer team. The effect should be smaller for targets fainter than HD209458. Further studies on this hypothesis are beyond the scope of the current paper.




 §.§ Spitzer/IRAC channels 3 and 4 - 5.8 and 8.0 μm



Figure <ref> shows the raw light-curves, best-fit models, and residuals for the observations taken with the Spitzer/IRAC channels 3 and 4.
For these channels we selected larger apertures: 9×9, 11×11, and 13×13 pixels. The reason is that the PSFs are larger for these channels, so that smaller (and static) apertures are more sensitive to the pointing variations. The experiment with the use of static apertures is to avoid the issues related to the determination of centroid coordinates (e.g., ). Here we focus on the results obtained by fitting a single exponential ramp with the first 49 minutes discarded from each dataset. The results obtained using different ramp parameterizations and discarded intervals are discussed in Appendix <ref>. We anticipate that the best-fit orbital parameters are insensitive to the choice of ramp parameterization and discarded interval.



Figure <ref> shows the tests of the aperture sizes for the observations taken with the Spitzer/IRAC channels 3 and 4. The tests show stronger trends in transit depth versus aperture size than those observed at the shorter wavelengths, and these trends decrease with the aperture size, revealing a positive background (see Section <ref>). The greater impact of background at the longer wavelengths is expected because of the higher zodiacal flux and the lower stellar flux <cit.>. At 5.8 μm, the projected transit depths with the background subtracted are consistent with the results without subtraction within 1 σ, because the error bars are also larger. At 8.0 μm, the background subtraction changes the best-fit transit depth by 3.4 and 1.6 σ for the two observations.

Figure <ref> and Table <ref> report our final results for the transit depth and the orbital parameters. We selected the 9×9 array, as it gives slightly smaller residuals and error bars in most cases. The results obtained using larger arrays are almost identical (see Appendix <ref>). All parameter estimates from observations at different epochs are mutually consistent to within 1 σ.









 §.§ Overview of Spitzer/IRAC results



Figure <ref> summarizes the results obtained for the four passbands, and compares them with those reported in the previous literature. The orbital parameters estimated from different epochs and passbands are consistent within 1 σ, except in two cases, in which they are consistent within 2 σ.

Our results are very similar to those reported by <cit.>, i.e., the differences are generally much smaller than 1 σ. Two exceptions are the two warm datasets at 3.6 μm, which are outliers in their analysis. Our results are also consistent within 1 σ with those reported by <cit.> for the orbital parameters, except at 8.0 μm, for which their reported values are inconsistent with those obtained in the other passbands. They also obtained larger transit depths at 5.8 and 8.0 μm, most likely due to different ramp parameterizations (see Appendix <ref>). It appears that the significantly smaller error bars reported by <cit.> do not fully account for the degeneracies between a/R_* and i. In fact, the results that they obtained with two different fitting algorithms are inconsistent in some cases, as is clearly shown in Figure <ref> (but the transit depths assume identical values). <cit.> analysed the 4.5 μm full phase-curve, obtaining transit parameters consistent with ours within 0.4 σ. 



 §.§ The Spitzer/IRAC prior information


Following the SEA BASS procedure as developed by <cit.>, we derive the information on the orbital parameters from the infrared measurements in the form of gaussian probability distributions, which may be used as Bayesian priors in subsequent analyses at shorter wavelengths.

The most common way to combine n multiple measurements of the same parameter, x_i ±σ_i, into a unique estimate, x̅±σ̅, is by computing the weighted mean:

    x̅ = ∑_i=1^n w_i x_i /∑_i=1^n w_i ,   σ̅= √(1/∑_i=1^n w_i ),

where the weights are w_i = 1/ σ_i^2. The underlying assumption when taking the weighted mean is that there are no systematic errors <cit.>. A more conservative estimate is obtained by computing the arithmetic mean and standard deviation. Unlike the weighted mean error, the standard deviation is sensitive to differential systematic errors across the measurements, and does not scale as n increases.
We calculated both the weighted and arithmetic means of a/R_* and i over the 10 Spitzer/IRAC measurements; the results are reported in Table <ref> and Figure <ref>. 

It is known that a/R_* and i are strongly correlated parameters (e.g., ). The information on their correlation is lost if we report only their mean values and relevant uncertainties, regardless of how they are calculated. A simple way to retain this information without increasing the mathematical complexity is to consider an equivalent set of less strongly correlated parameters to express in the standard form (x̅±σ̅). A possible set consists of the impact parameter, b, and the “central” transit duration, T_0 <cit.>:

    b = a/R_*cosi
    
    T_0 = P/πarcsin ( √(1-b^2)/a/R_*sini )

The impact parameter is the sky-projected planet-star separation in units of
the stellar radius at the instant of inferior conjunction; the central transit duration is the length of the interval during which the planet center covers the star, i.e., the projected separation is not larger than the stellar radius.

Figure <ref> and Table <ref> report the individual estimates, arithmetic and weighted means of b and T_0. These parameters are less widely scattered around their weighted means than a/R_* and i, as measured by the smaller reduced chi-squared values, χ_0^2.

Figure <ref> shows the scatter plots for the two sets of parameters (a/R_*,  i) and (b,  T_0). In the former set, the correlation is striking, while in the latter, it is not evident by eye. The overplotted regions determined using the new set of parameter priors (b and T_0) are also clearly more representative of the data distribution than the analogous definitions with a/R_* and i.




§ VISIBLE RESULTS









 §.§ HST/STIS G430L - white light-curve

Figure <ref> and <ref> show the two white light-curves integrated over the STIS/G430L passband before and after data detrending, the best-fit transit model and residuals. The ∼1.7×10^-4 rms amplitude of residuals is equal to the photon noise value. Figure <ref> reports the best parameter results for a/R_*, i, b and T_0 when adopting different priors (see Section <ref>). All parameter estimates are mutually consistent within 2 σ. The parameters obtained with gaussian priors on b and T_0 are mutually consistent well within 1 σ, and also have smaller error bars than those obtained with gaussian priors on a/R_* and i. Overall, we note that the results are closer to the weighted rather than arithmetic mean values of the Spitzer/IRAC results. In particular, when the arithmetic means of b and T_0 are taken as prior (cf., third points of each panel in Figure <ref>), the mean values of the posterior distribution are driven more than halfway toward the corresponding weighted means, while the widths are reduced by a factor of up to 6. Figure <ref> reports the analogous results in transit depth and for different apertures (widths of 17, 31, 39, and 49 pixels in the cross-dispersion direction). They are all consistent within ∼1 σ, and the smallest error bars are ∼50 ppm. From now on, we refer to the results obtained with the 17-pixel aperture and Spitzer/IRAC weighted mean priors on b and T_0.















The empirical limb-darkening profile is shown in Figure <ref>, with two theoretical profiles overplotted, derived from <cit.>. The theoretical profiles were computed ad hoc using 3D hydrodynamical stellar-atmosphere models (, ) or interpolated from a precalculated grid of 1D plane-parallel models (MARCS, ). The intensities predicted by the 3D model are consistent  within 1 σ with the empirical values inferred from our light-curve fit for μ>0.2, i.e., up to 98% of the projected stellar radius, and significantly overestimated near the edge of the projected disk. In particular, the limb-to-center intensity ratios are ≳0.15 and 0 for the theoretical and fitted models, respectively. This discrepancy has a small effect on the light-curve residuals, mainly at the ingress and egress, as shown in Figure <ref>. The impact on the inferred transit depth is negligible, as can be noted in Figure <ref>. The intensities predicted by the 1D model deviate significantly from the empirical results, leading to much higher systematic residuals and a greater difference in transit depth (just above 1 σ, however), as also shown in Figures <ref> and <ref>.

The discrepancy between 3D and 1D models for HD209458 has been described before in the literature <cit.>. The novelty in this paper is the precise measurement of the stellar limb-darkening profile from the data, without using any prior information about the star, stellar-atmosphere models, or precalculated tables of coefficients. The trade-off for this admission of ignorance is an ∼1.8 times larger error bar in transit depth.




 §.§ HST/STIS G430L - transmission spectrum

We fixed the impact parameter and transit duration to their Spitzer/IRAC weighted mean values (also supported by the white light-curve analysis), and the phase shift for each HST visit to those estimated from the white light-curve fit. The only free parameters in the spectral fits are the transit depth and the four stellar limb-darkening coefficients. Fixing the values of the wavelength-independent parameters is common practice in exoplanet spectroscopy, where the focus is on the differences in transit depth over multiple passbands, rather than on their individual transit depth values (e.g., ). The error bars do not account for potential offsets, uniform in wavelength, that may be caused by the adopted parameters.

Figure <ref> shows the detrended spectral light-curves, best-fit transit models, residuals, and differences between these models and those obtained with fixed limb-darkening coefficients reported by <cit.>. The empirical fits are able to eliminate time-correlated noise in the residuals, which appear if the limb-darkening coefficients are fixed to those adopted by <cit.> (see also Figure <ref>). The improvements are statistically significant for the three passbands with the highest S/Ns, leading to smaller rms residuals by 16–34% (see also Appendix <ref>).

Figure <ref> compares the empirical limb-darkening profiles with the reference models, showing, as expected, significant discrepancies over the three passbands with the highest S/Ns. Note that the reference limb-darkening models for the spectral bins are 1D, and the discrepancies are similar to those observed for the white light-curve analysis.
Figure <ref> reports the final transmission spectrum and those calculated under different assumptions. Numerical results for the selected configuration are reported in Table <ref>.  The spectrum does not change if the impact parameter, transit duration, and two-visit phase shifts are free to vary individually in all passbands, rather than fixed to the same values, but the error bars are typically twice larger. This test shows no evidence of systematic errors introduced by our choice of fixed parameters in the spectral fits. Some biases of tens of ppm are expected when using the fixed limb-darkening coefficients, but they are overcome by the error bars. Similar discrepancies will be statistically significant in the near future, thanks to the smaller error bars that are expected with future instruments onboard the James Webb Space Telescope (JWST; ) and other missions.








§ CONCLUSIONS

We performed a self-consistent reanalysis of infrared-to-visible transits of the exoplanet HD209458b, applying an upgrade of the Bayesian scheme proposed by <cit.>, here named SEA BASS. The SEA BASS procedure enables simultaneous investigation of the stellar and exoplanetary atmospheres by spectroscopic measurements of the stellar limb-darkening profile and the optically thick area of the transiting exoplanet. This paper reports the first application of the SEA BASS approach on real datasets. The upgrade on the previous version, tested on simulated datasets by <cit.>, consists of using more than one observation over multiple infrared passbands to assess the geometric parameters with higher confidence, and the use of a set of less degenerate parameters (impact parameter and transit duration, rather than semimajor axis and inclination). 

The improved data-detrending techniques lead to consistent parameter results for all the Spitzer/IRAC observations, including those that appear as outliers in the literature.

We measured the stellar limb-darkening profiles of the HD209458 star over multiple passbands in the wavelength range 290–570 nm. The retrieved profiles are in good agreement with the 3D  model computed by <cit.> for the whole passband, and present non-negligible deviations from the corresponding 1D MARCS model and the 1D ATLAS models reported by <cit.> for the spectral bins. The potential biases in transit depth that are due to the inaccurate limb-darkening models are not larger than a few tens of ppm, hence they are comparable to the 1 σ error bars.

The analysis presented in this paper confirms the reliability of the SEA BASS approach with the current datasets. Future missions, such as the JWST, are expected to provide higher quality datasets that will lead to smaller error bars, which in turn will increase the importance of SEA BASS for avoiding the potential biases of theoretical models. We also note that HD209458 is relatively well known, because it is a Sun-like star. The uncertainties in the theoretical models are larger for other stellar types, especially for cooler stars.


G.M. thanks G. Tinetti, A. Tsiaras, I. P. Waldmann, I. D. Howarth, D. Dicken, C. Danielski, M. Martin-Lagarde, and P.-O. Lagage for useful discussions.
This work was supported by the LabEx P2IO, the French ANR contract 05-BLAN-NT09-573739 and the ERC project ExoLights (617119).
This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.
Some of the data presented in this paper were obtained from the Mikulski Archive for Space Telescopes (MAST). STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.



























































§ SPITZER/IRAC BACKGROUND





Figure <ref> compares the background estimates obtained with the test of aperture sizes with classical estimates from two detector areas with no apparent astrophysical sources. The first area consists of four 5×5 arrays in the corners of the sub-array area, excluding the edges. For the full-array observations, this area was defined to have the same position relative to the centroid. Analogously, the second area consists of four 5×5 arrays at the edges of the full-array area (not applicable to sub-array observations). The three measurements are correlated over the multiple datasets, but the background measured within the sub-array area is systematically higher than measured at the corners of the full array, denoting some light contamination from the source. Owing to the limited size of the sub-array, this small contamination cannot be eliminated. The results of the test of aperture sizes are generally closer to the full-array estimates, in particular for the observations at 5.8 μm and the first two observations at 3.6 μm during the cold Spitzer era. In a few cases, the test suggests significantly lower background than the classical estimates.






Figures <ref> and <ref> report the transit parameters obtained using different aperture sizes with and without background corrections. Except for the earlier observation at 8.0 μm, all parameter results are consistent within 1 σ, as expected for observations that are not background limited. The correction at 8.0 μm is necessary to obtain consistent transit depth values over the two observations, especially when using larger apertures. The orbital parameters, a/R_* and i, are less strongly affected by the choice of the array and background correction.





§ THE COLD-WARM SPITZER OFFSET



We previously showed that warm Spitzer transit depth measurements are systematically smaller than the cold ones in the same passband, with discrepancies at the 2–5 σ level (see Section <ref>). <cit.> observed a similar trend, but they attributed the discrepancy to the presence of high-frequency noise of unknown origin during part of the warm Spitzer observations at 3.6 μm. The same explanation cannot be applied to the observations at 4.5 μm as they themselves stated, although the discrepancy is also less significant (≲2 σ). It is reasonable to assume that there is an offset in the calibration of the two Spitzer eras. <cit.> obtained inconsistent results for the orbital parameters in the 3.6 μm warm measurements. This problem is minimized by our pixel-ICA reanalysis, i.e., the new orbital parameter estimates are mutually consistent within 1 σ (see Section <ref>); the only inconsistency remains between the warm and cold transit depth measurements. We also noted that if we add appropriate offsets to the flux measurements in order to obtain identical transit depths over the cold and warm eras, the orbital parameters also become more similar (but they are always within 1 σ, see Figure <ref>).


Such offsets are more than eight times greater than the highest background estimate in the 3.6 and 4.5 μm channels. Any potential offset equal for all pixels can be discarded, as it would be revealed by the test of aperture sizes (see Section <ref>). We show that a possible explanation is a differential error in the nonlinearity correction between the cold and warm Spitzer eras, which can be mimicked to first order by an offset in the total flux.

If the instrument response is linear, the flux integrated over the photometric aperture is proportional to the source flux:

    F_0(t) = α F_* f(t) = α F_* (1-λ(t)) ,

where F_* is the unperturbed stellar flux, and λ(t) is the occulting factor as a function of time (λ(t)=0 out-of-transit, 0<λ(t) ≪ 1 during the transit). We refer to f(t) = 1-λ(t) as the “transit model” (cf., ). Note that in the simple case defined by Equation <ref>, the transit model is coincident with the observed flux normalized to the out-of-transit value:

    F_0,(t) = F_0(t)/α F_* = f(t) .

If the instrument response is nonlinear, the proportionality constant in Equation <ref>, α, must be replaced with a non-trivial function of the source flux. We assume a first-order nonlinearity:

    F_1(t) = α (1-β F_* f(t))  F_* f(t) .

The F_1 time series normalized to the out-of-transit value is

    F_1,(t) = F_1(t)/α F_* (1-β F_*) = 1-β F_* f(t)/1-β F_* f(t) .

In the limit β F_* ≪ 1, using the second-order Taylor expansion for the denominator, and expliciting f(t) = 1-λ (t):

    F_1,(t) ≈ (1+β F_* + (β F_*)^2) × (1-β F_* (1-λ (t))) × (1-λ (t)) .

By expanding the products, and retaining up to second-order terms in β F_* and λ:

    F_1,(t) ≈ 1 - λ (t) + β F_*  λ (t) .



Now, we consider the case with linear instrument response and an additive offset, C:

    F_(t) = α F_* (1-λ(t)) + C,

The F_ time series normalized to the out-of-transit value is:

    F_,(t) = F_(t)/α F_* + C =  1-λ(t) + C/α F_*/1 + C/α F_* .

In the limit C / (α F_*) ≪ 1, using the second-order Taylor expansion for the denominator:

    F_,(t) ≈ ( 1 - C/α F_* +  ( C/α F_* )^2  ) × ( 1-λ(t) + C/α F_* ) .

By expanding the products, and retaining up to second-order terms in C / (α F_*) and λ:

    F_,(t) ≈ 1 - λ (t) + C/α F_*λ (t)



Note that Equations <ref> and <ref> have the same structure, and are equivalent if

    β F_* = C/α F_*≡ x .

The third terms in Equation <ref> and <ref>, with λ (t) = (R_p/R_*)^2, approximate the bias in transit depth due to the nonlinearity or the additive offset, respectively. As the observed cold-warm discrepancy is ∼300 ppm in transit depth (see Section <ref>), using (R_p/R_*)^2 ∼ 1.5 × 10^-2, we obtain

    x ≈ 2 × 10^-2 .

This result is consistent with the ∼3% error in the absolute flux calibration reported by <cit.>.





§ DIFFERENT RAMP MODELS FOR SPITZER/IRAC CHANNELS 3 AND 4


In this Appendix, we compare the performances of the following ramp models in statistical terms:


  * single exponential,

    F(t) = a_0 - a_1 e^-t/τ_1 ,


  * double exponential,

    F(t) = a_0 - a_1 e^-t/τ_1 - a_2 e^-t/τ_2 ,


  * linear + quadratic log,

    F(t) = a_0 + a_1 t + a_2 log(t-t_0) + a_3 [ log(t-t_0) ]^2  ,


  * linear + logarithmic,

    F(t) = a_0 + a_1 t + a_2 log(t-t_0) ,


  * quadratic log,

    F(t) = a_0 + a_1 log(t-t_0) + a_2 [ log(t-t_0) ]^2  ,


  * fifth-order polynomial,

    F(t) = a_0 + a_1 t + a_2 t^2 + a_3 t^3 + a_4 t^4 + a_5 t^5  .


Following the prescription by <cit.>, we set t_0 to be 30 s before the observations begin to prevent the divergence of the logarithm at t=0.

We focus on the two observations at 5.8 and 8.0 μm  taken in 2008 July. 

Regardless of the choice of the ramp parameterization, the root mean square (rms) of residuals decreases monotonically while discarding more data points from the beginning until random fluctuations become dominant.  The best-fit parameters also follow a trend before stabilizing.
Figure <ref> shows the best transit parameters obtained with six different ramp models (Equations <ref>–<ref>), discarding different time intervals from the beginning. When the single exponential formula (Equation <ref>) is used, the estimated transit depth becomes stable after the first 10–20 minutes are discarded, i.e., it varies by less than 1 σ if more data are discarded. The orbital parameters are even more robust, as the amplitudes of the observed trends across 0–91 minute cuts are much smaller than 1 σ. The results obtained with the other ramp parameterizations show similar or larger variations as a function of the discarded interval, the worst case is the fifth-order polynomial (Equation <ref>).  Overall, the orbital parameters, a/R_* and i, are very stable with respect to the choice of the ramp model and the length of the discarded interval. Only for the case of the polynomial ramp at 8.0 μm and discarded interval shorter than 30 minutes do they fall outside the accepted ranges. The linear + quadratic log and the double exponential formulas lead to systematically larger transit depths at 8.0 μm than those obtained with the single exponential. <cit.> adopted the linear + quadratic log formula to correct for the ramp at 8.0 μm. We conclude that the most likely cause of discrepancy with our analysis in this passband is the choice of the ramp parameterization. The statistical tests discussed in this Appendix are moderately in favor of the single exponential form.

We consider two quantitative model selection techniques: the Bayesian Information Criterion (BIC, ) and the Akaike Information Criterion (AIC/AIC_c, ). They have been calculated under the assumptions of uniform uncertainties per data point and gaussian distributed residuals [It is equivalent to rescaling the Spitzer-provided uncertainties, a common, well-motivated practice in the analysis of transits <cit.>.]:

    BIC = n logσ_res^2 + k logn


    AIC_c = ( n logσ_res^2 + 2k ) + 2k(k+1)/n-k-1 = AIC + 2k(k+1)/n-k-1

Here σ_res is the rms of residuals, n is the number of data points, and k is the number of ramp parameters. The first term in Equations <ref> and <ref> is a measure of the goodness of fit, while the second terms are penalties proportional to the number of free parameters, and the third term in Equation <ref> is a correction for small sample sizes (negligible for n ≫ k). In our cases, the differences between AIC and AIC_c  are ∼0.01 (negligible), therefore we can here equivalently refer to AIC instead of AIC_c. The best models according to these selection techniques are those with the lowest BIC or AIC <cit.>, respectively. BIC tends to select simpler models than AIC as a result of the different penalties for the number of parameters (greater for BIC). The relative strength of evidence between two models, i and j, is defined by their probability ratio <cit.>:

    ( p_i/p_j )_BIC = exp ( - 1/2 (BIC_i - BIC_j)  )  = exp ( - Δ BIC/2 ) 
     ( p_i/p_j )_AIC = exp ( - 1/2 (AIC_i - AIC_j)  )  = exp ( - Δ AIC/2 )

According to the interpretation provided by <cit.>, absolute differences smaller than 2 in BIC/AIC[<cit.> only refers to ΔBIC, but Equations <ref>–<ref> allow to extend the same interpretation to ΔAIC.] are not statistically significant and only differences greater than 6 provide “strong evidence” against the model with higher BIC/AIC.

Table <ref> reports the differences in BIC and AIC obtained for different ramp models relative to the single exponential. Positive differences denote a preference for the single exponential model, negative differences show a preference for the alternative model. The single exponential is the favorite model for the ramp at 5.8 μm according to both selection criteria, but the linear + logarithmic and the quadratic log models have nearly the same likelihoods. The situation is reversed for the ramp at 8.0 μm with a slight preference for the linear + logarithmic model, but there is no strong evidence against the single exponential. The double exponential formula at 8.0 μm is strongly disfavored by BIC, but only marginally disfavored by AIC. The fifth-order polynomial is by far the worst parameterization for all ramps (the same is valid for lower-order polynomials).


We also checked for correlated noise in the residuals, but all the autocorrelation coefficients turned out to be smaller than 0.07, regardless of the formula adopted to remove the ramp. Figure <ref> shows that the residuals scale with the bin size according to the expected behavior for gaussian noise, except when the polynomial ramp is adopted.








§ ACCURACY OF TWO-COEFFICIENT LIMB-DARKENING LAWS


<cit.> claimed that two-coefficient limb-darkening laws are inaccurate in the 290-570 nm wavelength range covered by HST/STIS G430L, even if the corresponding light-curve fit residuals appear to be limited by photon noise. This result was based on simulated data using PHOENIX stellar-atmosphere models <cit.>.
Figures <ref> and <ref> confirm the same behavior for real HD209458 datasets. In particular, we obtained nearly identical light-curve residuals when using free quadratic <cit.>, power-2 <cit.>, or claret-4 <cit.> coefficients, despite that the corresponding stellar intensity profiles are, in some cases, significantly different. The bias in transit depth can be ∼100 ppm when adopting two-coefficient limb-darkening laws, and in particular, it is slightly above 1 σ for the 403–457 and 458–512 nm passbands for the quadratic, power-2, or square-root law. We note that for the case analyzed in this paper, even the 1D claret-4 theoretical models outperform the empirical two-coefficient parameterizations (see also Figure <ref>). 




































[Agol et al.(2010)]agol10 Agol, E., Cowan, N. B., Knutson, H. A., et al.  2010, , 721, 1861
[Akaike(1974)]akaike74 Akaike, H.  1974, IEEE Transactions on Automatic Control, 19, 716
[Allard et al.(2012)]allard12 Allard, F., Homeier, D., & Freytag, B.  2012, RSPTA, 370, 2765
[Ballerini et al.(2012)]ballerini12 Ballerini, P., Micela, G., Lanza, A. F., & Pagano, I.  2012, , 539, A140
[Barman(2007)]barman07 Barman, T.  2007, , 661, L191
[Beaulieu et al.(2010)]beaulieu10 Beaulieu, J. P., Kipping, D. M., Batista, V., et al.  2010, , 409, 963
[Beichman et al.(2014)]beichman14 Beichman, C., Benneke, B., Knutson, H., et al.  2014, , 126, 1134
[Berta et al.(2012)]berta12 Berta, Z. K., Charbonneau, D., Désert, J.-M., et al.  2012, , 747, 35
[Blecic et al.(2013)]blecic13 Blecic, J., Harrington, J., Madhusudhan, N., et al.  2013, , 779, 5
[Brown(2001)]brown01 Brown, T. M.  2001, , 553, 1006
[Brown et al.(2001)]brown01b Brown, T. M., Charbonneau, D., Gilliland, R. L., Noyes, R. W., & Burrows, A.  2001, , 552, 699
[Burnham & Anderson(2002)]burnham02 Burnham, K. P., & Anderson, D. R.  2002, Model Selection and Multimodel Inference: A practical information-theoretic approach (2nd ed.), Springer-Verlag, ISBN 0-387-95364-7
[Carey et al.(2012)]carey12 Carey, S. J., Ingalls, J. G., Hora, J. L., et al.  2012, , 8442, 84421z
[Carter et al.(2008)]carter08 Carter, J. A., Yee, J. C., Eastman, J., Gaudi, B. S., & Winn, J. N.  2008, , 689, 499
[Cavanaugh(1997)]cavanaugh97 Cavanaugh, J. E.  1997, Statistics & Probability Letters, 33, 201
[Charbonneau et al.(2000)]charbonneau00 Charbonneau, D., Brown, T. M., Latham, D. W., & Mayor, M.  2000, , 529, L45
[Charbonneau et al.(2002)]charbonneau02 Charbonneau, D., Brown, T. M., Noyes, R. W., & Gilliland, R. L.  2002, , 568, 377
[Charbonneau et al.(2008)]charbonneau08 Charbonneau, D., Knutson, H. A., Barman, T., et al.  2008, , 686, 1341
[Claret(2000)]claret00 Claret, A.  2000, , 363, 1081
[Csizmadia et al.(2013)]csizmadia13 Csizmadia, S., Pasternacki, T., Dreyer, C., et al.  2013, , 549, A9
[Désert et al.(2009)]desert09 Désert, J.-M., Lecavelier des Etangs, A., Hébrard, G., et al.  2009, , 699, 478
[Fazio et al.(2004)]fazio04 Fazio, G. G., Hora, J. :L., Allen, L. E., et al.  2004, , 154, 10
[Evans et al.(2015)]evans15 Evans, T. M., Aigrain, S., Gibson, N., et al.  2015, , 451, 680
[Fischer et al.(2016)]fischer16 Fischer, P. D., Knutson, H. A., Sing, D. K., et al.  2016, , 827, 19
[Grillmair et al.(2012)]grillmair12 Grillmair, C. J., Carey, S. J., Stauffer, J. R., et al.  2012, Proc. SPIE Operations: Strategies, Processes, and Systems IV, 8448, 84481I
[Gustafsson et al.(2008)]gustafsson08 Gustafsson, B., Edvardsson, B., Eriksson, K., et al.  2008, , 486, 951
[Haario et al.(2006)]haario06 Haario, H., Laine, M., & Saksman, E.  2006, Stat. Comput., 16, 339
[Harrington et al.(2007)]harrington07 Harrington, J., Luszcz, S., Seager, S., Deming, D., & Richardson, L. J.  2007, , 447, 691
[Hayek et al.(2012)]hayek12 Hayek, W., Sing, D., Pont, F., & Asplund, M.  2012, , 539, A102
[Henry et al.(2000)]henry00 Henry, G. W., Marcy, G. W., Butler, R. P., & Vogt, S. S.  2000, , 529, L41
[Hyvärinen & Oja(2000)]hyvarinen00 Hyvärinen, A., & Oja, E.  2000, Neural Networks, 13, 411
[Hora et al.(2014)]hora14 Hora, J. L., Witzel, G., Ashby, M. L. N., et al.  2014, , 793, 120
[Ingalls et al.(2016)]ingalls16 Ingalls, J. G., Krick, J. E., Carey, S. J., et al.  2016, , 152, 44
[IRAC Instrument Handbook(v. 2.1.2)]irachand IRAC Instrument and Instrument Support Teams  2015, IRAC Instrument Handbook, v. 2.1.2
[Kipping & Tinetti(2010)]kipping10 Kipping, D. M., & Tinetti, G.  2010, , 407, 2589
[Knutson et al.(2008)]knutson08 Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T.  2008, , 673, 526
[Knutson et al.(2007)]knutson07 Knutson, H. A., Charbonneau, D., Allen, L. E., et al.  2007, , 447, 183
[Knutson et al.(2009)]knutson09 Knutson, H. A., Charbonneau, D., Cowan, N. B., et al.  2009, , 703, 769
[Kopal(1950)]kopal50 Kopal, Z.  1950, HarCi, 454, 1
[Krick et al.(2012)]krick12 Krick, J. E., Glaccum, W. J., Carey, S. J., et al.  2012, , 754, 53
[Krick et al.(2016)]krick16 Krick, J. E., Ingalls, J., Carey, S., et al.  2016, , 824, 27
[Lust et al.(2014)]lust14 Lust, N. B., Britt, D., Harrington, J., et al.  2014, , 126, 1092
[Machalek et al.(2008)]machalek08 Machalek, P., McCullough, P. R., Burke, C. J., et al.  2008, , 684, 1427
[Mandel & Agol(2002)]mandel02 Mandel, K., & Agol, E.  2002, , 580, L171
[Morello(2015)]morello15b Morello, G.  2015, , 808, 56
[Morello et al.(2016)]morello16 Morello, G., Waldmann, I. P., & Tinetti, G.  2016, , 820, 86
[Morello et al.(2017)]morello17 Morello, G., Tsiaras, A., Howarth, I. D., & Homeier, D.   2017, , 154, 111
[Morello et al.(2014)]morello14 Morello, G., Waldmann, I. P., Tinetti, G., et al.  2014, , 786, 22
[Morello et al.(2015)]morello15 Morello, G., Waldmann, I. P., Tinetti, G., et al.  2015, , 802, 117
[Nordlund & Galsgaard(1995)]nordlund95 Nordlund, Å., & Galsgaard, K.  1995, A 3D MHD Code for Parallel Computers, Tech. rep., Astronomical Observatory, Copenhagen University
[Raftery(1995)]raftery95 Raftery, A. E.  1995, Sociological Methodology, 25, 111
[Schwarz(1978)]schwarz78 Schwarz, G.  1978, The Annals of Statistics, 6, 461
[Schwarzschild(1906)]schwarzschild1906 Schwarzschild, K.  1906, Nachr. Königl. Ges. Wiss. Göttingen. Math.-Phys. Kl., p. 43
[Seager & Mallén-Ornelas(2003)]seager03 Seager, S., & Mallén-Ornelas, G.  2003, , 585, 1038
[Seager & Deming(2010)]seager10 Seager, S., & Deming, D. 2010, , 48, 631
[Sing et al.(2011)]sing11 Sing, D. K., Pont, F., Aigrain, S., et al.  2011, , 416, 1443
[Sing et al.(2016)]sing16 Sing, D. K., Fortney, J. J., Nikolov, N., et al.  2016, , 529, 59
[Stevenson et al.(2010)]stevenson10 Stevenson, K. B., Harrington, J., Nymeyer, S., et al.  2010, , 464, 1161
[STIS Instrument Handbook(v. 16.0)]stishand Riley, A.  2017, STIS Instrument Hanbook for Cycle 25, Version 16.0
[Taylor(1996)]taylor96 Taylor, J. R.  1996, An Introduction to Error Analysis (2nd ed.; University Science Books)
[Tichavský et al.(2008)]tichavsky08 Tichavský, P., Koldovský, Z., Yeredor, A., Gómez-Herrero, G., & Doron, E.  2008, ITNN, 19, 421
[Tinetti et al.(2007)]tinetti07 Tinetti, G., Vidal-Madjar, A., Liang, M.-C., et al.  2007, , 448, 169
[Torres et al.(2008)]torres08 Torres, G., Winn, J. N., & Holman, M. J.  2008, , 677, 1324
[Tsiaras et al.(2017)]tsiaras17 Tsiaras, A., Waldmann, I. P., Zingales, T., et al.  2017, arXiv:1704.05413
[Vidal-Madjar et al.(2003)]vidal-madjar03 Vidal-Madjar, A., Lecavelier des Etangs, A., Désert, J.-M., et al.  2003, , 422, 143
[Vidal-Madjar et al.(2004)]vidal-madjar04 Vidal-Madjar, A., Désert, J.-M., Lecavelier des Etangs, A., et al.  2004, , 604, L69
[Vidal-Madjar et al.(2013)]vidal-madjar13 Vidal-Madjar, A., Huitson, C. M., Bourrier, V., et al.  2013, , 560, A54
[Waldmann(2012)]waldmann12 Waldmann, I. P.  2012, , 747, 12
[Wu et al.(2010)]wu10 Wu, X., Roby, T., & Ly, L.  2010, procspie, 7737, 773716
[Zellem et al.(2014)]zellem14 Zellem, R. T., Lewis, N. K., Knutson, H. A., et al.  2014, , 790, 5
































































































