Probing the atmosphere of a subJovian planet orbiting a cool dwarf
Abstract
We derive the 0.01 m binned transmission spectrum, between 0.74 and 1.0 m, of WASP80b from low resolution spectra obtained with the FORS2 instrument attached to ESO’s Very Large Telescope. The combination of the fact that WASP80 is an active star, together with instrumental and telluric factors, introduces correlated noise in the observed transit light curves, which we treat quantitatively using Gaussian Processes. Comparison of our results together with those from previous studies, to theoretically calculated models reveals an equilibrium temperature in agreement with the previously measured value of 825K, and a subsolar metallicity, as well as an atmosphere depleted of molecular species with absorption bands in the IR (). Our transmission spectrum alone shows evidence for additional absorption from the potassium core and wing, whereby its presence is detected from analysis of narrow 0.003 m bin light curves (). Further observations with visible and nearUV filters will be required to expand this spectrum and provide more indepth knowledge of the atmosphere. These detections are only made possible through an instrumentdependent baseline model and a careful analysis of systematics in the data.
keywords:
instrumentation: spectrographs – techniques: spectroscopic – planetary systems – planets and satellites: atmospheres – planets and satellites: individual: WASP80b1 Introduction
The remarkable progress made to date in detecting and characterising extrasolar planets has been made possible through dedicated space and groundbased facilities. In a mere couple of decades we have gone from discovering the first of these socalled exoplanets (Campbell et al., 1988; Wolszczan & Frail, 1992; Mayor & Queloz, 1995) to characterising their physical properties, the most fascinating of which has been the detection of their atmospheres. This is particularly an intriguing feature of these alien worlds, as it provides the means for understanding the mechanisms involved in formation and evolution of planets (Mordasini et al., 2012a, b; Dorn et al., 2015), as well as presenting the opportunity for detection of biomarkers (Kaltenegger & Traub, 2009; Snellen et al., 2013; Benneke & Seager, 2012), pointing to perhaps biological processes on those planets capable of harbouring lifeforms producing them.
One method through which the atmospheric envelope around a transiting exoplanet is detected is transmission spectroscopy, where minute, wavelengthdependent variations of the transit depth are measured from modelling the spectrophotometric transit light curves (Seager & Sasselov, 2000). This method generally probes the upper layers of the atmosphere, where lower pressure levels together with shorter optical paths lead to transmission of stellar radiation through the exoatmosphere. This process leaves distinct spectral imprints on the observed radiation, one consequence of which is wavelengthdependent relative planetary radius (obtained from measuring the transit depth).
Such studies are performed with observations either from space (Ehrenreich et al., 2007; Sing et al., 2011; Deming et al., 2013; Wakeford et al., 2013; Knutson et al., 2014; Sing et al., 2016) or with groundbased facilities (Narita et al., 2005; Gibson et al., 2012a; Stevenson et al., 2014; Sedaghati et al., 2015, 2016; Lendl et al., 2016; Mallonn & Strassmeier, 2016), both of which have their own advantages. Spacebased observations have the advantage of not being affected by atmospheric extinction, contamination and turbulence, and therefore benefit from having the entire electromagnetic range as probing domain. For instance, NASA’s James Webb Space Telescope (JWST; Gardner et al., 2006) will be able to probe exoatmospheres in the infrared wavelengths, where transmission spectroscopy from the ground is rather difficult. On the other hand, groundbased facilities benefit from telescopes with large collecting mirrors, which is essential for performing spectroscopy at high time and spectral resolution.
ESO’s FOcal Reducer and low dispersion Spectrograph (FORS2; Appenzeller et al., 1998) mounted at the Cassegrain focus of Unit Telescope 1 (UT1) of the Very Large Telescope (VLT) has been the instrument of choice to perform such observations (Bean et al., 2010; Bean et al., 2011; Sedaghati et al., 2015, 2016; Lendl et al., 2016; Nikolov et al., 2016). It offers two multiobject spectroscopic modes: (1) Mask eXchange Unit (MXU, custom designed and lasercut masks), (2) MOS (MultiObject Spectroscopy, movable slitlets via 19 pairs of arms). The option of wide slits, multiple reference star selection for telluric correction in conjunction with the great light collecting power of the 8.2m telescope, make this instrument ideal for performing differential spectrophotometric observations, required here.
Previously, it was determined that the Longitudinal
WASP80b is a gas giant transiting a cool, possibly late Ktype 11.88V magnitude dwarf, with a mass of 0.554 0.035 M and radius of 0.952 0.026 R, orbiting its host star with a 3.068 day period (Triaud et al., 2013). It is one of only a handful of gas giants orbiting a latetype dwarf host (e.g. WASP43b, Kepler45b, HATP54b, HATS6b; Hellier et al., 2011; Johnson et al., 2012; Bakos et al., 2015; Hartman et al., 2015) and has a dayside temperature within the Tdwarf range (Triaud et al., 2015). These authors also found a transmission spectrum indistinguishable from a flat line and no evidence of active region crossing by the planet (Triaud et al., 2013; Mancini et al., 2014), despite the stellar spectral analysis indicating high levels of activity (Mancini et al., 2014).
In this work, we present a red transmission spectrum of this planet from light curves that are somewhat compromised by systematic effects. We will study the role of telluric and instrumental effects in introducing such discrepancies in the light curves, before consideration of any astrophysical phenomena.
2 Data analysis
2.1 Observations
A single transit of WASP80b was observed using Antu, the 8.2m UT1 of the VLT with the FORS2 instrument on the 16th of June 2013, with the data taken as part of the programme 091.C0377(C) (PI: Gillon), before the LADC upgrade. The programme itself comprised four transit observations, two of which were not performed due to bad weather, while the light curves from the fourth transit suffered from severe systematic effects and are subsequently not used for the analysis in this work. The instrument has a 6.8 field of view and is equipped with two detectors, where the redoptimised MIT set was utilised for the observations. The instrument was used in the MXU mode, which means that a custom designed mask with 10″wide slits positioned on the target and comparison stars, was placed in the light path. This essentially acts as a blocking mask to free the CCD for recording of the simultaneous spectra. The 600z grism (with the order sorter filter OG590) was used as the dispersing element, yielding a spectral range of m, although the exact wavelength coverage for each target is dependent on the horizontal location of its slit on the CCD.
The entire science observational sequence lasted 4.87h, with the first frame taken at 05:17UT and the last taken at 10:04UT. The complete transit lasted 2.26h, with first contact at 06:17UT and the last one at 8:32UT. The LADC was left in park position during the observing sequence, with the two prisms fixed at their minimal separation distance of 30mm. The standard 100kHz readout mode was used, which together with exposure time of 25s (apart from the first seven frames where adjustments were made in order to reach the optimal value), yielded 277 exposures, 130 of which were during the transit. The conditions were clear throughout the night and the seeing varied between 0.97″ and 2.35″. The field started at airmass of 1.24, rose to 1.08 and the last frame was taken at airmass 1.45. A copy of the mask was created for the purpose of wavelength calibration, with narrow 1″ slits centred on the science slits. Bias, flatfield and arclamp (for the purpose of wavelength calibration) images were taken, as part of the routine daytime calibration sequences, before and after the observations.
2.2 Data reduction
We used our previously written reduction pipeline based on PyRAF (Sedaghati et al., 2016), which is optimised to the science goals of transmission spectroscopy. In addition to the standard steps of overscan and bias shape subtraction, spectral flatfielding and wavelength calibration, sky background subtraction and cosmic ray contamination removal, the pipeline also includes steps for optimizing the size of the extraction box used to obtain the onedimensional stellar spectra, performed independently for each star. We also experiment with the optimal extraction algorithm (Horne, 1986) and check to see if this in any way improves the quality of the extracted spectra, which was not the case. Once the wavelength calibrated spectra are obtained, the pipeline makes corrections to the dispersion solution calculated for each target. This is first done for all the spectra of each star through optimisation of the cross correlation function, and then the spectra of the stars are further corrected with respect to each other to ensure a consistent solution across all apertures. This is a necessary and important step due to the low resolution of the spectra. An example set of these wavelength calibrated and corrected spectra is shown in Figure 1.
2.3 Broadband light curves & systematic effects
To obtain the broadband light curves, we integrate the series of spectra for each star within the largest possible common domain (shown as the grey area in Figure 1), which are shown in Figure 2(a). At first glance we observe an airmassdependent trend for all the targets, as well as the clear transit signal of WASP80b. Initially we correct all the light curves for both of these factors, the results of which are shown in Figure 2(b). These values are later on used to search for possible physical variants responsible for some of the systematic trends in the final differential transit light curve. We also produce differential transit light curves of WASP80 with respect to all the observed comparison stars, shown in Figure 2(c). This shows a transit light curve that is very precise and its deviations from a transit model are mostly due to systematic trends. Most significantly, we observe an almost Vshape transit, which could be interpreted as crossing of a stellar spot by the exoplanet, as was observed for WASP52b (Kirk et al., 2016). However, this claim cannot be made until all possible sources of correlated noise in the data have been considered.
Instrumental effects
Prior to the upgrade of the FORS2 instrument, where the LADC prisms (Boffin et al., 2015) were replaced, inhomogeneities in the degraded antireflective coating of the old unit were causing differential transmission through the telescope optics. This would then manifest itself as systematic flux variations in the differential light curves. We now take a closer look at this effect.
Boffin et al. (2016) presented the improvement in the transmission of light through the optical elements of FORS2 by comparing a stack of flatfield images from before and after the LADC exchange. For our analysis of the systematic trends attributed to the old degraded LADC (used for the observations here), we use their stacked flatfield image taken with the R_SPECIAL+76 filter, shown in Figure 3. The individual twilight flatfields used to construct this frame were taken approximately a year after our observations. To this effect, the exact initial position of the field of view relative to the LADC configuration, inferred through observation of flatfield inhomogeneities, is not precisely known. Therefore, we can only try to estimate this position from possible trends in the light curves. In Figure 3 we plot the observed stars at this estimated starting position and trace the path of each star through the entire observing sequence, where every 10th exposure is shown with a dot.
To obtain the possible starting point of the field of view relative to the LADC setup, we rotate the configuration shown in Figure 3 for 360 at 1 steps to cover all possible initial positions. If the LADC deficiencies are the cause of systematic trends in the light curve, then at the true path of the field there will be a spike in the calculated correlation between the light curve of each star (the ones that are corrected for airmass and transit, ref. Fig 2(b)) and the flatfield value along the path taken by the star, shown in 3 and read in 5(a,b). This calculated correlation as a function of field rotation angle is shown in Figure 4, where the indicated peak at 86 starting angle could point to a possible relative orientation of the LADC at the beginning of the observing sequence
In order to study the impact of the LADC inhomogeneities on the measured light curves, we read the normalised flatfield value along the path of each star, for this given starting position, which is shown in Figure 5(a). We use the normalised flatfield, as its variations are a direct indication of the problems caused by the degraded LADC coating. Additionally, since we are also interested in the correction of the differential transit light curves, as a second order correction, we also calculate the ratio of these readings to the values read along the path of WASP80. This is shown in Figure 5(b). The calculated correlation between the light curve residuals and the flatfield variations due to the degraded LADC is rather weak, which means that the LADC is not a major contributor to the systematic trends. However, it must be noted that the flatfield images used to construct the composite image in Figure 3 were taken more than a year after the epoch of the observations analysed in this work. Therefore, our conclusion about the impact of the LADC deformities is rather tentative.
A further possible cause of correlated noise is the spatial stability of the instrument. We characterise the positioning of stellar spectra on the detectors along two axes, spatial and spectral ( & ). The former variation is measured by taking the mean of the midpoints of several fitted Gaussian profiles to the twodimensional images across the entire frame in the spatial direction, which are plotted relative to the initial frame in Figure 5(e). We measure the latter in a similar way, but instead fit several telluric absorption features in the onedimensional spectra, shown in Figure 5(c). Studying the relations between both of these variables indicated no significant correlation with the flux variations, which is perhaps due to the stability of the instrument and the data reduction procedure accounting for any remaining residual effects.
Telluric effects
In addition to the instrumental effects mentioned previously, telluric conditions also play a significant role in introducing timecorrelated noise into the measured timeseries data. Namely, variations in the measure of atmospheric seeing, due to turbulence, cause changes in the stellar point spread function (PSF) along both the dispersion and the physical axes. We measure the former as the FWHM of Gaussian profile fits to several telluric absorption lines in the one dimensional spectra, variations of which are shown in Figure 5(d) and denoted as “<sky> fwhm”. This is essentially the convolution of the instrumental profile and the true stellar profile, and therefore a measure of variations in the stellar psf given the instrument response is constant. The FWHM of the PSF along the physical axis is also measured through recording the characteristics of several Gaussian function fits to various columns (orthogonal to the dispersion axis) of the two dimensional spectra. These are shown in Figure 5(f), as “psf fwhm” in units of arcseconds, which is another approach to measuring the seeing variations.
We compare the variations of seeing measured separately for each star (from the psf FWHM) to the normalised flux after the correction for the airmass trend. For WASP80, additionally, the transit feature is also removed for these calculations, in order to maintain a constant degree of freedom in all correlation plots. Such relations are plotted in Figure 6, where the comparison stars 5 and 6 have not been included due to their significantly low SN. This will be the case for the remainder of our analysis.
For the measure of correlation, we again use Pearson’s parameter described briefly earlier, the calculated values for which are given in each panel. Additionally, we calculate the 1, 2 and 3 confidence intervals, shown as different shades of grey for visualization purposes. The levels of correlation between individual flux variations and atmospheric seeing are consistent across all the stars and point to a moderate association between the two parameters. The relation between each pair of parameters is calculated using an orthogonal least squares approach (Isobe et al., 1990), using the odr package from python’s scipy (Jones et al., 2001). This is preferred to the standard linear regression approach due to the present uncertainties in determining “psf FWHM”. We use the gradients of these relationships, given in Figure 6, to transform the measured seeing values to flux variations for each star, which in turn will be used as an input of our systematic model, to be described shortly. We choose this approach instead of constructing a correction model to avoid the introduction of biases in our parameter determination and error estimation.
2.4 Light curve models
We model the transit light curves with the implementation of the analytic solution of Mandel &
Agol (2002) as the mean function. Additionally, we estimate the correlated noise as a Gaussian Process (GP), which provides a modelindependent stochastic method for the inclusion of the systematic component into our model (Rasmussen &
Williams, 2006). The application of GPs to modelling exoplanet transit light curves and transmission spectroscopy was introduced by Gibson et al. (2012b). The strength of this approach in accounting for systematic trends was highlighted in the detection of hazes in the atmosphere of the exoplanet HD189733b, where auxiliary information from the HST were used to train the GP model (Gibson
et al., 2012c). Subsequently, we use their GeaPea and Infer modules
In the Bayesian framework, we write the likelihood function, , as a matrix equation in order to introduce offdiagonal elements to the covariance,
(1) 
where is the residual vector with being the transit model, is the input parameter matrix ( being the number of inputs and being the number of data points or flux measurements), and are collections of noise and transit model parameters respectively, and represents the covariance matrix. Ideally we would like to make every element of the covariance matrix a free parameter in our model, but the size of the matrix () makes this rather impractical. Hence we “model” those elements using a covariance function, more commonly referred to as the kernel. With the aid of this kernel, the covariance matrix is written as , being an input of the kernel, the Kronecker delta and the variance term, the last two of which ensure the addition of Poisson or white noise to the diagonal of the covariance matrix.
For the kernel, we use the Squared Exponential (SE)
(2) 
where is the maximum covariance and are inverse scale parameters for all the input vectors (essentially the columns of the matrix). This parameter determines the length of the “wiggles” in the function. We choose this form for the kernel as it is the defacto default form for GPs as it is infinitely differentiable and easy to integrate against most functions.
Finally, in the definition of the analytical transit function, we use the quadratic limbdarkening (Kopal, 1950) law to describe the centre to limb variations of brightness across the stellar disk. We made this choice after the comparison of the Bayesian Information Criterion (BIC; Schwarz et al., 1978) values for models of higher complexity, such as the three parameter (Sing et al., 2009) or the nonlinear Claret (2000) laws.
Broadband
We model the best two broadband differential transit light curves shown in Figure 2(c), which are produced by using the comparison stars 1 and 2 (top two), as the other light curves are either too noisy or suffer from large systematic effects. These light curves are obtained by the integration of the relevant spectra within the domain that is shaded in grey in Figure 1. In addition to the transit function, we also include a quadratic term (as a function of parallactic angle) in the description of the analytical transit model, to account for the telescope rotation dependent effects in the differential light curves.
To find the maximum posterior solution, , we need to optimise the Bayesian relation,
(3) 
which is true up to a constant term, where the final term on the right hand side of the equation is the prior probability assumed for the parameters of the complete model. We generally use uninformative priors for the transit parameters, with the exception of the linear and quadratic limb darkening coefficients (), the orbital period, , and the eccentricity, , according to:
(4) 
which fixes the period and eccentricity of the orbit to what has previously been determined by Triaud et al. (2015), and also ensures positive brightness across the stellar disk. For the noise model parameters, we set the following prior probability distributions,
(5) 
to ensure positive values for all the parameters, as well as using a gamma distribution
For the starting values of the fitted parameters, we use the NelderMead simplex algorithm (Nelder & Mead, 1965) to find an initial set of optimal solutions for the transit and noise parameters. To find the maximum posterior solution, we optimise the log posterior given in Equation 3. We obtain these parameter posterior distributions by using the MonteCarlo MarkovChains method to explore the joint posterior probability distribution of our multivariate models (Collier Cameron et al., 2007; Winn et al., 2008; Gibson et al., 2012b). We run four independent MCMC simulations with 100 walkers, of 100 000 iterations each. Once the chains are computed, we extract the marginalised posteriors for the free parameters to check for mutual convergence and for possible correlations among any pairs of parameters. The transit parameters that are fitted for are midtransit time (), scaled semimajor axis (), relative planetary radius (), impact parameter (), limbdarkening coefficients (), the three coefficients of the baseline model. This model is a second order polynomial of the parallactic angle, q, whose complexity is chosen using again the BIC selection rule. This baseline model was initially chosen as a quadratic in time, which would describe low frequency variations due to the colour difference between the target and the comparison star. However, this approach results in out of transit flux fit that is worse than using the parallactic angle, and additionally leads to an overestimation of planetary radius as compared to previous results from photometry in our wavelength region. In addition to these, we also fit for the noise model parameters ().
In order to decide what input parameters to include in the description of the GP model, as well as which of the two final broadband light curves suffer less from systematic effects, we model both of these broadband data series with our previously described model using up to three input parameters for the kernel. We did not find any evidence of correlation with any other optical state parameters, to justify their initial inclusion into our model. To compare the fitted models, we calculate their BIC values which includes a penalty term for addition of complexity to a model. These calculated values are given in Table 2.4.1, where their comparison indicates that our systematic model is best able to describe the broadband light curve relative to the comparison star 1, when time, seeing and the ladc inhomogeneity are all used to model the covariance matrix. However, it must be noted that there is some evidence against this model relative to the case that does not include the LADC variations (BIC6), but we choose to include this final parameter due to the small increase in the reduced chisquared statistic (). The correlations and posterior distributions of the parameters of this noise model are given in Figure 7.

Footnotes
 pubyear: 2017
 pagerange: Probing the atmosphere of a subJovian planet orbiting a cool dwarf–4
 Or Linear as both adjectives are used in the definition of the LADC.
 It must be noted that this value is 3 offset from what is calculated based on the orientation of the stacked flatfield image.
 Written in Python programming language and available freely from the following repository: https://github.com/nealegibson
 Also known as the Radial Basis Function
 The probability density function for a gamma distribution of a variable is given by .