Assessing telluric correction methods for Na detections with high-resolution exoplanet transmission spectroscopy
Adam B. Langeveld, Nikku Madhusudhan, Samuel H. C. Cabot, Simon T. Hodgkin
MMNRAS , 1–12 (2021) Preprint 15 January 2021 Compiled using MNRAS L A TEX style file v3.0
Assessing telluric correction methods for Na detections withhigh-resolution exoplanet transmission spectroscopy
Adam B. Langeveld, ★ Nikku Madhusudhan, † Samuel H. C. Cabot , Simon T. Hodgkin Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge, CB3 0HA, UK Yale University, 52 Hillhouse, New Haven, CT 06511, USA
Accepted 2021 January 7. Received 2020 December 24; in original form 2020 August 27
ABSTRACT
Using high-resolution ground-based transmission spectroscopy to probe exoplanetary atmo-spheres is difficult due to the inherent telluric contamination from absorption in Earth’satmosphere. A variety of methods have previously been used to remove telluric features inthe optical regime and calculate the planetary transmission spectrum. In this paper we presentand compare two such methods, specifically focusing on Na detections using high-resolutionoptical transmission spectra: (a) calculating the telluric absorption empirically based on theairmass, and (b) using a model of the Earth’s transmission spectrum. We test these methods onthe transmission spectrum of the hot Jupiter HD 189733 b using archival data obtained with theHARPS spectrograph during three transits. Using models for Centre-to-Limb Variation andthe Rossiter-McLaughlin effect, spurious signals which are imprinted within the transmissionspectrum are reduced. We find that correcting tellurics with an atmospheric model of theEarth is more robust and produces consistent results when applied to data from different nightswith changing atmospheric conditions. We confirm the detection of sodium in the atmosphereof HD 189733 b, with doublet line contrasts of − . ± .
07 % (D2) and − . ± .
07 %(D1). The average line contrast corresponds to an effective photosphere in the Na line locatedaround 1 . 𝑅 𝑝 . We also confirm an overall blueshift of the line centroids corresponding to netatmospheric eastward winds with a speed of 1 . ± . − . Our study highlights the impor-tance of accurate telluric removal for consistent and reliable characterisation of exoplanetaryatmospheres using high-resolution transmission spectroscopy. Key words:
Planets and satellites: atmospheres – Atmospheric effects – Techniques: spectro-scopic – Methods: observational – Planets and satellites: individual: HD 189733 b
The era of characterising exoplanet atmospheres has acceleratedgreatly over the last two decades. With technological advancementsin instrumentation, we are now able to monitor exoplanetary sys-tems and characterise in detail both their bulk properties as well astheir atmospheres. Numerous planets have been observed with bothradial velocity and transit methods, allowing for extensive analy-ses and characterisation (e.g. Fischer et al. 2014; Sing et al. 2016;Madhusudhan 2019; Pinhas et al. 2019; Welbanks et al. 2019). Ob-servations of transiting systems can give an insight into the chemicalproperties of the planetary atmosphere. During one orbit there aretwo windows which can be used for such studies: the primary tran-sit, when the planet passes in front of its host star as viewed fromEarth, and the secondary eclipse when the planet is occulted by thestar. In both instances, the planetary atmospheric spectrum can be ★ E-mail: [email protected] † E-mail: [email protected] deduced by measuring the change in observed flux outside of andduring transit (Seager & Sasselov 2000), and absorption or emis-sion due to chemical species may be seen. Such inferences can alsobe made if phase-resolved spectra are obtained at other parts of theorbit (Stevenson et al. 2014).The first observation of an exoplanet atmosphere was madeby Charbonneau et al. (2002). Using the STIS spectrograph on the
Hubble Space Telescope (HST) with a resolution of 𝑅 = 𝑅 ∼ , 𝑅 ∼ , © a r X i v : . [ a s t r o - ph . E P ] J a n A. Langeveld et al. almost 20 different chemical species in exoplanet atmospheres todate (Madhusudhan 2019). Among the most common are detectionsof the alkali species Na and K in hot Jupiters (e.g. Charbonneauet al. 2002; Redfield et al. 2008; Sing et al. 2016; Sedaghati et al.2016; Nikolov et al. 2016; Wyttenbach et al. 2017; Chen et al. 2018;Casasayas-Barris et al. 2018; Jensen et al. 2018; Deibert et al. 2019;Seidel et al. 2019; Hoeijmakers et al. 2019).Many challenges are presented when using ground-based fa-cilities to acquire high-resolution spectra of transiting exoplanets.For methods which derive the planetary transmission spectrum bycomparing the in-transit and out-of-transit stellar spectra, observ-ing times must be chosen to ensure that an adequate out-of-transitbaseline can be obtained. An insufficient out-of-transit baseline canadversely affect the quality of the transmission spectrum (Wytten-bach et al. 2015). A significant problem arises from the contamina-tion of spectra due to absorption in the Earth’s atmosphere. Thesetelluric lines can merge with stellar lines, or appear at a similar orhigher strength as features in the planetary transmission spectrum(Redfield et al. 2008; Snellen et al. 2008). In the optical domain, thepredominant sources of absorption are telluric water and oxygen.Therefore, data must be obtained on nights with low water vapourcontent and when the airmass is low to minimise the impact of thecontamination.Alongside the telluric removal process, it is essential to makefurther corrections to account for stellar reflex motion, systemicvelocity, and the planetary radial velocity. Additionally, a planetwhich passes in front of a rotating star blocks different amounts ofblueshifted and redshifted light. This causes the shape of lines in theintegrated stellar spectrum to change as the transit progresses, andis known as the Rossiter-McLaughlin (RM) effect (Rossiter 1924;McLaughlin 1924; Queloz et al. 2000; Triaud 2018). The spectralline shape is also affected by Centre-to-Limb Variation (CLV) whichdescribes the brightness of the stellar disk as a function of limb an-gle (Czesla et al. 2015; Yan et al. 2017). Together with telluriccontamination, these effects may imprint spurious signals withinthe planetary transmission spectrum and must be corrected for toprevent false identifications of chemical species. A recent studyhas shown that absorption features in the transmission spectrum ofHD 209458 b can be explained purely by the induced CLV and RMsignals (Casasayas-Barris et al. 2020), highlighting the importanceof correcting for these effects. In the near future, newly developedhigh-resolution spectrographs will be used to target rocky planetsand super-Earths. It is therefore crucial to be able to accuratelytackle the many problems associated with ground-based observa-tions in order to distinguish spectral features of these planets. Fur-ther characterisation will be possible when these results are used incombination with upcoming James Webb Space Telescope (JWST)data.Several methods have been used in previous work to removetelluric features from high-resolution spectra. For the first detectionof the atmosphere of HD 189733 b from a ground-based facility(Redfield et al. 2008), the spectrum of a telluric standard star wasrecorded immediately after the observations. A spectrum of ab-sorption in the Earth’s atmosphere was recovered by using a stellartemplate to remove the weak stellar lines of the rapidly rotatinghot B-star. This was also employed in subsequent studies (Jensenet al. 2011, 2012). Other datasets have been corrected by assum-ing a linear relationship between airmass and telluric line strength,and deriving a telluric spectrum to remove variation throughoutthe night (Vidal-Madjar et al. 2010; Astudillo-Defru & Rojo 2013;Wyttenbach et al. 2015). More recently, studies have used modelsof Earth’s atmosphere to correct for the telluric spectrum, e.g.
TER- RASPEC (Lockwood et al. 2014), molecfit (Smette et al. 2015;Kausch et al. 2015; Allart et al. 2017), and other custom-built codes(e.g. Yan et al. 2015; Casasayas-Barris et al. 2017). Contaminationcan also be corrected by removing the time-dependent variationof flux at each wavelength element using linear regression (e.g.Snellen et al. 2010; Brogi et al. 2012), or by removing system-atic trends with singular value decompositions (SVDs) or principalcomponent analysis (PCA) (e.g. de Kok et al. 2013; Birkby et al.2013; Piskorz et al. 2016). Comparison has also shown that somemethods are better at correcting H O than O (Ulmer-Moll et al.2019).In this study we investigate how different telluric correctionmethods may affect the measurements of chemical species in exo-planet atmospheres at optical wavelengths. We see inconsistencies inmeasured results depending on the method employed and the nightlyweather variation, so it is important to quantify which techniquesare best for analysing high resolution exoplanet spectra. We chooseto compare two popular methods: one which derives the tellurictransmission solely from the data, and one which uses a model ofmolecular absorption in Earth’s atmosphere. Both of these methodsaccount for the variation of telluric line strength over the multiple-hour observing window. This gives an advantage over modelling atelluric standard star, where the resultant spectrum is insensitive toatmospheric variation during the transit and requires extra time toobserve. We focus specifically on the quality of telluric correctionsin the region around the sodium D-lines, however the analysis isapplicable to the full range of the spectrograph. Other methods dis-cussed above may be more suited for Doppler-resolved moleculardetections (particularly in the near-infrared) but were beyond thescope of this work.We focus on HD 189733 b – a tidally-locked hot Jupiter orbit-ing a bright K1V host star with magnitude V = 7.67. There havealready been a number of atmospheric detections for this planet: Na(Redfield et al. 2008; Jensen et al. 2011; Huitson et al. 2012; Wyt-tenbach et al. 2015; Khalafinejad et al. 2017), H O (Birkby et al.2013; McCullough et al. 2014; Brogi et al. 2016, 2018; Cabot et al.2019), H (Jensen et al. 2012; Lecavelier Des Etangs et al. 2010;Bourrier et al. 2013), CO (de Kok et al. 2013; Brogi et al. 2016;Cabot et al. 2019), He (Salz et al. 2018), and HCN (Cabot et al.2019). Additionally, there has been evidence of Rayleigh scatteringand high-altitude hazes (Pont et al. 2008; Lecavelier Des Etangset al. 2008; Di Gloria et al. 2015), and measurement of atmosphericwind speeds (Louden & Wheatley 2015).In section 2 we give an overview of the high-resolution HARPSspectra used for our analysis. We discuss the data reduction and twomethods for removing telluric contamination in section 3, and thecalculation of the transmission spectrum which is common for bothmethods in section 4. In section 5 we evaluate the differences in thetelluric reduction processes by measuring the absorption parame-ters of the Na doublet using Gaussian fits and a binned passbandanalysis. We also measure the net atmospheric wind speed fromthe blueshift of the line profiles and discuss further atmosphericproperties. The results are briefly summarised in section 6 togetherwith a conclusion about the most advantageous method.
In this paper we use observations of high-resolution transmissionspectra of the hot Jupiter HD 189733 b as a case study. Datawere obtained from programs 072.C-0488(E), 079.C-0828(A) (PI:Mayor) and 079.C-0127(A) (PI: Lecavelier des Etangs) using the
MNRAS , 1–12 (2021) elluric correction of exoplanet transmission spectra Night Date
Table 1.
HARPS observing log of HD 189733. Nights 1–3 are numberedto remain consistent with Wyttenbach et al. (2015). The number of in andout-of-transit exposures are calculated using the mid-exposure times and theephemeris in Table A1. -1000100-1000100-0.25 -0.20 -0.15 -0.10 -0.05 0.00 0.05 0.10 0.15 0.20 0.25
Orbital Phase ( ) -1000100 P l a n e t R a d i a l V e l o c i t y ( k m s ) Figure 1.
Planetary radial velocity curves for observations on night 1 (top),2 (middle), and 3 (bottom). The dashed line shows the modelled radialvelocity curve, and the phase coverage of in and out-of-transit obaservationsare shown with blue and orange lines respectively. The transit window (timebetween first and fourth contact) is shaded.
High Accuracy Radial velocity Planet Searcher (HARPS) echellespectrograph on the ESO 3.6 m telescope in La Silla, Chile, andmade available through the ESO archive. With a resolving powerof 𝑅 = , 𝜇 mwhich gives a 1 arcsec aperture on the sky (Mayor et al. 2003). Allobservations are reduced by the HARPS Data Reduction Software(DRS) v3.5. The pipeline performs calibrations using the dark, biasand flat-field frames taken at the beginning of the night, corrects forthe blaze, and uses Thorium-Argon (Th-Ar) reference exposures forwavelength calibration. Finally, each spectrum is remapped onto auniform grid at a resolution of 0.01 Å in the solar system barycentricrest-frame.High-resolution spectra of HD 189733 were obtained withHARPS during four primary transits and are available through theESO archive. The observing log in Table 1 shows the dates andprogram IDs of the observations. The number of in-transit spec-tra were calculated using the mid-exposure time and ephemeris ofAgol et al. (2010). A combination of these datasets was first used to measure the Rossiter-McLaughlin effect (Triaud et al. 2009).They have since been analysed by several other authors to studyother atmospheric properties, leading to important results such as:a detection of sodium in the planetary atmosphere (Wyttenbachet al. 2015), spatially resolved eastward winds (Louden & Wheat-ley 2015), a temperature gradient of ∼ . − (Heng et al.2015), a measurement of the Rayleigh scattering slope (Di Gloriaet al. 2015), and an attempt to constrain the presence of atmosphericwater vapour (Allart et al. 2017). As discussed by many of theseauthors, on 29 July 2006 there are only 12 observations in total, andthere are none during the second half of the transit because of badweather conditions. We were unable to recover a clear planetarytransmission spectrum due to the lack of data before and duringtransit, and therefore chose to ignore this night.Parameters of the HD 189733 system which were used through-out this study are shown in Table A1. The known stellar radial ve-locity semi-amplitude was used to model the stellar and planetaryradial velocity curves. Spectra with mid-exposure times between thefirst and fourth contacts of the eclipse are defined as in-transit. Asseen in Figure 1, the observations on nights 1 and 3 fully cover theperiod shortly before, during, and after the transit. The orange andblue lines indicate out of and in-transit observations respectively.The complete transit duration is shown by the shaded region. Preliminary corrections need to be made to the spectra before tel-luric contamination can be considered. We extract the data from theone-dimensional HARPS (s1d) files. The full spectral wavelengthrange is 3781–6912 Å, however we limit this to 4000–6800 Å toreduce systematics at the edges of the full spectrum where strongtelluric contamination or low throughput is evident. A small numberof pixels in each spectrum contain substantially higher flux valuesthan the median, which are unlikely to be from the observed as-trophysical source. To correct for these, a mask is created for eachspectrum by applying a median filter with a width of 9 pixels andsubtracting this from the original spectrum. Any pixels in the maskwhich are greater than 10 standard deviations from the median areflagged, and the corresponding fluxes in the original spectra arecorrected to the median of the 10 surrounding pixels. Each spec-trum is then normalised by dividing by its median. The errors in theflux values are assumed to be dominated by photon noise and wepropagate the Poisson uncertainties throughout the analysis.All ground-based observations in this wavelength range arepolluted with telluric lines predominantly due to H O and O . Theselines vary in strength with observing conditions such as airmassand precipitable water vapour. Incorrect removal of tellurics willcreate false signals in the resultant planetary transmission spectrum.We now consider the removal of telluric contamination using twomethods. We first employed a method similar to that described by Vidal-Madjar et al. (2010), Astudillo-Defru & Rojo (2013), and Wytten-bach et al. (2015) – by assuming that telluric lines have a linearvariation with airmass. The airmass at the start and end of eachobservation is measured on-site and stored in the headers of theHARPS files – we assigned an average airmass to each observationusing these two measurements. Choosing only the starting or endingvalue caused a slight difference in the telluric correction. For each
MNRAS000
MNRAS000 , 1–12 (2021)
A. Langeveld et al. T () Wavelength ( Å ) R e l a t i v e F l u x Figure 2.
The effect of re-scaling the strength of telluric lines using theempirical telluric reference spectrum.
Top panel: derived telluric spectrumwith prominent features around 5899.9 Å and 5901.3 Å.
Bottom panel: un-corrected (grey) and re-scaled data (blue). The telluric lines in each spectrumare reduced so that they appear as though observations took place with aconstant airmass. Stellar lines remain unchanged. wavelength step in the spectrum, a linear fit between the naturallogarithm of the measured fluxes and their corresponding airmasswas made, taking the form 𝑙𝑛 ( 𝐹 ( 𝜆 )) = ( 𝑁 𝑘 𝜆 ) 𝑠 + 𝑐 , (1)where 𝐹 ( 𝜆 ) is the flux at wavelength 𝜆 , 𝑁 𝑘 𝜆 is the zenithal opticaldepth, and 𝑠 is the airmass. The telluric lines vary with airmass butthe exoplanet absorption lines do not. Therefore a telluric referencespectrum, 𝑇 ( 𝜆 ) , can be built using the value of the gradient inequation 1 (Wyttenbach et al. 2015) 𝑇 ( 𝜆 ) = 𝑒 𝑁 𝑘 𝜆 . (2)It is important to only use out-of-transit spectra to build the telluricreference spectrum; using in-transit spectra would also correct forabsorption due to the planet’s atmosphere. Therefore, it is idealto obtain a number of spectra shortly before and shortly after thetransit to ensure a large enough out-of-transit baseline. However,observations on night 2 began when the planet was already in transitand there are not enough out-of-transit observations to produce agood quality telluric reference spectrum. We were unable to recovera clear transmission spectrum using the limited out-of-transit data,so chose to include the in-transit sample as well. We note that someof the planetary signals could be over-corrected and muted as aresult.Using equation 2, the telluric contamination in the observedspectra can be adjusted so that they appear to have been observed ata constant reference airmass, 𝑠 ref . The adjusted spectra, 𝐹 ref , werecalculated for all observations using 𝐹 ref ( 𝜆 ) = 𝐹 obs ( 𝜆 ) 𝑇 ( 𝜆 ) ( 𝑠 − 𝑠 ref ) , (3)where 𝐹 obs ( 𝜆 ) is the observed spectrum and 𝑠 is the airmass at thetime of the observation (Wyttenbach et al. 2015). The referenceairmass was taken to be the mean airmass of the in-transit sam-ple only. An example of the adjustment to telluric lines is shown inFigure 2. The top panel shows the calculated telluric reference spec-trum T ( 𝜆 ) , and the bottom panel shows the result of the correction with the original data in grey and the corrected data in blue. It isnoted that the telluric lines are only adjusted to a level of constantairmass but are not removed completely. Following this correction,the transmission spectrum was computed as described in section 4.Residuals of telluric lines may occasionally be seen in the trans-mission spectrum despite performing the correction derived fromairmass. These lines may be confused for planetary atmospheric ab-sorption, but are actually caused by changes in Earth’s atmosphericwater content in the air above the telescope throughout the night.Similarly to Wyttenbach et al. (2015), we make a second correctionto account for this. For the data on night 3, a linear fit between thetransmission spectrum and the T ( 𝜆 ) was made. The transmissionspectrum was divided by this fit which effectively reduces partsof the spectrum where telluric residuals are correlated with T ( 𝜆 ) .We split up the full spectrum into smaller wavelength ranges andperformed the correction on each section to prevent variations inunrelated parts of the spectrum from affecting each other. Only oneiteration of the correction was needed to reduce the residuals to thecontinuum level. This is further discussed in section 5.1. molecfit An ideal method would reduce telluric contamination down to thenoise level without changing any of the existing data outside of theseregions. The empirical method discussed in the previous sectionstruggles to do this, and gets worse with low signal-to-noise.Here we consider telluric removal by modelling out Earth’satmospheric contribution. We use molecfit v1.2.0 – an ESO toolspecifically designed to correct for telluric contamination in ground-based spectra (Smette et al. 2015; Kausch et al. 2015).
Molecfit fits a line-by-line radiative transfer model (LBLRTM) of telluricabsorption to the observed spectra, producing a unique model of at-mospheric absorption for each observation. This was first performedon HARPS spectra by Allart et al. (2017), and has since been used inother HARPS studies (Cauley et al. 2017; Seidel et al. 2019; Cabotet al. 2020), as well as with other high-resolution spectrographs(Cauley et al. 2017, 2019; Allart et al. 2018, 2019; Nortmann et al.2018; Salz et al. 2018; Casasayas-Barris et al. 2019, 2020; Hoeij-makers et al. 2019; Alonso-Floriano et al. 2019; Kirk et al. 2020;Chen et al. 2020).HARPS spectra are given in the Solar System barycentric rest-frame. We use the Barycentric Earth Radial Velocity (BERV) valuesto shift each spectrum into the rest-frame of the telescope to ensurecorrect telluric modelling with molecfit . We then inspect areaswith heavy H O and O absorption around 5950, 6300 and 6475 Åand select 15 small regions (no larger than 2 Å) which contain onlytelluric lines and flat continuum (Figure 3). It is important that nostellar lines are included as they will affect the fit to the atmosphericmodel. The chosen regions remain constant for the duration of onenight of observing, however they must be re-selected for differentnights since the location of telluric lines changes with respect to thestellar lines.We provide molecfit with parameters similar to those dis-cussed by Allart et al. (2017), together with the date, location,and atmospheric conditions (e.g. humidity, pressure, temperature)stored within the HARPS output. The atmospheric model is fittedto the selected regions. With reference to the output parameters, wethen run the calctrans tool which fits an atmospheric model tothe entire spectrum, resulting in a unique telluric profile for eachspectrum with the same resolution. Finally, the spectra are dividedby their corresponding model to reduce all telluric lines down to thecontinuum noise level, as shown in Figure 4. MNRAS , 1–12 (2021) elluric correction of exoplanet transmission spectra Wavelength ( Å ) R e l a t i v e F l u x Figure 3.
Example of parts of the spectra which are contaminated with telluric lines. The 15 shaded light blue/green regions are passed to molecfit . Toppanel: modelled telluric spectrum (red) with H O absorption in bands around ∼ ∼ absorption around ∼ Middle panel: data (black) and the chosen regions (shaded) with widths <2 Å.
Bottom panel: expanded view of the green shaded region from the middle panel: by comparisonwith the molecfit model (red), we make sure that each region only contains telluric lines with no other stellar features. N o r m a li s e d F l u x A i r m a ss Wavelength ( Å ) Figure 4.
Result of telluric corrections with molecfit . The spectra for night 3 are overlaid and coloured according to their respective airmass, and the reddashed line shows the telluric model for the first observation.
Left panel: uncorrected data with clear telluric contamination around 5901.3 Å.
Right panel: corrected data following division of the model. The telluric lines have been reduced to the noise level.
Here we discuss our calculation of the transmission spectra fromobservations. This work addresses several common effects that canimprint false signatures onto the planetary transmission spectrumat high-resolution. These effects include telluric contamination, ve-locity contributions from stellar, systemic and planetary sources,as well as spurious signals due to the Rossiter-McLaughlin (RM)effect and Centre-to-Limb Variation (CLV).
Following removal of telluric contamination (sections 3.1 and 3.2),we extract the transmission spectrum while making corrections for the stellar reflex motion, systemic velocity, and planetary radialvelocity using a similar technique to Wyttenbach et al. (2015). In-dividual spectra are identified as 𝑓 ( 𝜆, 𝑡 in ) for in-transit exposuresand 𝑓 ( 𝜆, 𝑡 out ) for out-of-transit exposures by modelling the orbitusing the parameters listed in Table A1. As discussed in section 2,it is important that the observations cover the period shortly before,during, and shortly after the transit to create a good out-of-transitbaseline.The HARPS-measured stellar radial velocities at the time ofeach exposure are extracted from the Cross Correlation Function(CCF) files. An example of these measurements for observations onnight 3 is shown in Figure 5. The dashed line is a fit for both the MNRAS000
Following removal of telluric contamination (sections 3.1 and 3.2),we extract the transmission spectrum while making corrections for the stellar reflex motion, systemic velocity, and planetary radialvelocity using a similar technique to Wyttenbach et al. (2015). In-dividual spectra are identified as 𝑓 ( 𝜆, 𝑡 in ) for in-transit exposuresand 𝑓 ( 𝜆, 𝑡 out ) for out-of-transit exposures by modelling the orbitusing the parameters listed in Table A1. As discussed in section 2,it is important that the observations cover the period shortly before,during, and shortly after the transit to create a good out-of-transitbaseline.The HARPS-measured stellar radial velocities at the time ofeach exposure are extracted from the Cross Correlation Function(CCF) files. An example of these measurements for observations onnight 3 is shown in Figure 5. The dashed line is a fit for both the MNRAS000 , 1–12 (2021)
A. Langeveld et al. -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04
Orbital Phase ( ) -2.28-2.26-2.24-2.22-2.20-2.18-2.16 S t e ll a r R a d i a l V e l o c i t y ( k m s ) Figure 5.
Stellar radial velocity measurements for night 3, as measured byHARPS. In-transit observations are shown with blue diamonds, and out-of-transit with orange circles. Deviation from the modelled radial velocity curve(dashed line) during the shaded in-transit window is due to the Rossiter-McLaughlin effect. stellar reflex motion and the systemic velocity, taking the form 𝑣 ∗ = 𝐾 ∗ sin ( 𝜋𝜙 ) + 𝑣 sys , (4)where 𝐾 ∗ is the stellar radial velocity semi-amplitude, 𝜙 is the phaseand 𝑣 sys is the systemic velocity. For the in-transit spectra, there is adeviation from this fit due to the Rossiter-McLaughlin effect (Triaudet al. 2009; Di Gloria et al. 2015). Each spectrum is Doppler shiftedto correct for its corresponding measured stellar radial velocity andlinearly interpolated back to the uniform 0.01 Å grid. This cor-rects for the apparent radial velocity change due to the RM effect,however the overall shape of the transmission line profiles will stillbe affected and could lead to spurious signals – additional correc-tions are discussed in section 4.2. We simultaneously calculate theplanetary radial velocity for all observed phases using 𝑣 𝑝 = 𝐾 𝑝 sin ( 𝜋𝜙 ) , (5)where 𝐾 𝑝 = − 𝐾 ∗ ( 𝑀 ∗ / 𝑀 𝑝 ) .The maximum change in measured stellar radial velocitythroughout one night of observing is ∼
90 ms − , correspondingto a ∼ . ±
16 kms − duringthe transit. As a result, signals in the planetary transmission willbe shifted up to ∼ .
32 Å (32 resolution elements) from the restwavelength. It is therefore important to correct for this in order torecover the full transmission spectrum.We propagate the Poisson uncertainties of the raw data and usethe inverse variance as weights when combining the spectra. Thisreduces the contribution of spectra with low signal-to-noise whichmay have been affected by poor weather conditions or other time-dependent systematics. We co-add all out-of-transit exposures tocreate a master out-of-transit spectrum, denoted 𝑓 out ( 𝜆 ) . In-transitspectra contain signals from both the star and planet, whereas theout-of-transit spectra contain just the stellar signal. Classically, thetransmission spectrum can be calculated by dividing a master in-transit spectrum, 𝑓 in ( 𝜆 ) , by 𝑓 out ( 𝜆 ) , however this does not accountfor the planetary radial velocity shift. We therefore compute indi-vidual transmission spectra as (cid:60)( 𝜆, 𝑡 in ) = 𝑓 ( 𝜆, 𝑡 in ) 𝑓 out ( 𝜆 ) . (6)Each spectrum is then continuum normalised using a third-order polynomial to remove effects from instrumental systematics orweather variations. A Doppler shift using 𝑣 𝑝 is performed to moveeach spectrum into the planet’s rest-frame, and the final combinedtransmission spectrum is thus computed: (cid:60) (cid:48) ( 𝜆 ) = ∑︁ 𝑡 in (cid:60)( 𝜆, 𝑡 in )| 𝑣 𝑝 ( 𝑡 in ) − = 𝐹 𝜆, in 𝐹 𝜆, out − , (7)where 𝐹 𝜆,𝑖𝑛 is the Doppler corrected and time averaged in-transitspectrum. We remove any remaining broadband continuum varia-tions with a median filter of width 1501 pixels (15.01 Å). CLV and RM effects can change the shape of lines in the stellarspectrum during the transit. This induces spurious signals into thetransmission spectrum which must be corrected for (Queloz et al.2000; Triaud 2018). We model the extent of these effects in theresultant transmission spectra in a similar manner as Casasayas-Barris et al. (2019). We generate synthetic spectra of HD 189733using the software
Spectroscopy Made Easy (SME) (Valenti &Piskunov 1996). The synthesis routine takes as input: 1) physicalstellar parameters, for which we used those listed in Table A1 (wedo not provide a stellar rotation at this point); 2) a list of spectralline positions and strengths, for which we use the VALD3 database(Ryabchikova et al. 2015); and 3) a list of 𝜇 -angles (Mandel & Agol2002); we choose 21 angles spanning from the limb to the centre ofthe stellar disk. For each night of observation, we simulate the transitof the planet on an 80 ×
80 pixel grid centred on the host star. Eachpixel overlapping the stellar disk is assigned a synthetic spectrum,interpolated at the appropriate 𝜇 -angle from the set of calculated SME spectra. The spectrum is also Doppler shifted according to thelocal radial velocity of the star, per parameters in Table A1. Theintegrated flux over all pixels provides our simulated out-of-transitspectrum. At each phase of the transit, corresponding to exposuremidtimes, we geometrically determine the projected position of theplanet on the stellar disk (Cegla et al. 2016), and integrate thestellar flux over all pixels not occulted by the planet. The result isour simulated in-transit spectra. The jointly-modelled CLV and RMeffect on the observed transmission spectra is obtained by dividingeach simulated in-transit spectrum by the simulated out-of-transitspectrum. We do not account for potential non-LTE effects in thecores of the Na doublet, nor variable planetary radii accounting forabsorption by the atmosphere.As shown in the middle panel of Figure 6, the modelled trans-mission spectrum contains features which arise strictly from CLVand RM effects. However, these features have amplitudes (cid:46) . − smaller net blueshift in the Na doubletwhich underscores the potential impact of CLV and RM residuals, aswell as the RM velocimetric effect. We correct for small changes inthe shape of the line profiles by dividing by the model, and producea final, observed transmission spectrum by co-adding the overallnightly CLV and RM corrected spectra. MNRAS , 1–12 (2021) elluric correction of exoplanet transmission spectra -0.020.000.02-0.020.000.025886 5888 5890 5892 5894 5896 5898 Wavelength ( Å ) O b s e r v e d P h a s e R e l a t i v e F l u x Figure 6.
Comparison of the transmission spectrum and the modelled CLVand RM effects around the Na doublet for night 1. The colour scale in panels1 and 2 represents relative flux from 0.98 (black) to 1.02 (white), and the bluedashed lines indicate the transit window.
Top panel: planetary transmissionspectrum for each observed phase.
Second panel: modelled transmissionspectrum assuming no atmosphere, showing artefacts of the CLV and RMeffects.
Third panel: binned (20x) CLV/RM uncorrected (magenta) andcorrected (black) transmission spectrum in the planetary rest-frame. Thered dotted lines in panels 2 and 3 give the position of the sodium doublet inthe planetary rest frame.
In this section we show a comparison of the calculated HD 189733 btransmission spectra using the two different telluric corrections. Weconfirm and discuss the detection of Na i in the planetary atmo-sphere, assess the shape of spectral features and the amount ofabsorption, and give context for further inferred atmospheric prop-erties. Since we are focusing on the Na doublet, we restrict thewavelength range of the transmission spectra to 5870–5916 Å.
For a clearer visualisation of the spectral features, the transmissionspectra for all nights were binned to a 0.2 Å resolution (20 points perbin). We define absorption features as parts of the spectrum whichare negative in value. Figure 7 shows the transmission spectrum fornight 3 using the airmass correction method. On this night, a sec-ondary correction was performed to reduce the effects of water col-umn variation. The middle panel shows the transmission spectrumwithout the extra correction – by comparing this with the empiricaltelluric reference spectrum (top panel), it is clear that there are stillresidual telluric lines present around 5886.0, 5887.5, 5891.5, 5900.0and 5901.5 Å. Since these are of comparable strengths to the signalswe expect in the planetary transmission spectrum, the validity of anydetections is significantly reduced. The bottom panel shows the re-sult after applying the second correction, and the residual telluricshave been reduced without affecting other parts of the spectrum.However, the sodium lines could be slightly over-corrected due toincreased noise in the telluric reference spectrum in these regions.This correction was not necessary for nights 1 and 2.To assess the quality of the telluric reduction methods, we fitGaussian profiles to the lines in the sodium doublet which providesa good approximation of the line cores (Gandhi & Madhusudhan2017). We used the
LevMarLSQFitter module from astropy , Wavelength ( Å ) -1.5-1.0-0.50.00.5 T () ' ( % ) ' ( % ) Figure 7.
Effect of telluric corrections on the transmission spectrum for night3.
Top panel: empirical telluric reference spectrum.
Middle panel: transmis-sion spectrum after airmass telluric corrections.
Bottom panel: transmissionspectrum and after applying a second correction for water column variation.The full resolution spectra are shown in grey, and binned 20x in black. Bycomparing the middle panel with the derived telluric reference spectrum(top panel), residual telluric lines can be seen in regions around 5886.0,5887.5, 5891.5, 5900.0 and 5901.5 Å. These are removed following theextra correction.Line 𝜆 (Å) D (%) FWHM (Å)D2 a . ± . − . ± .
07 0 . ± . a . ± . − . ± .
07 0 . ± . m . ± . − . ± .
07 0 . ± . m . ± . − . ± .
07 0 . ± . Table 2.
Measured properties of the Gaussian fits to the sodium D-lines inthe combined transmission spectrum of HD 189733 b. Lines are labelledwith subscripts a and m to denote which telluric removal method (airmassor molecfit ) was used. which performs a Least-Squares fit using a Lavenberg-Marquardtalgorithm and accounts for errors in the flux which have been prop-agated from the photon noise on the raw spectra. We define threeproperties of the Gaussian fit: the amplitude or line contrast, D ;centroid 𝜆 ; and the full width at half maximum, FWHM. Figure 8shows the night 3 transmission spectrum around the Na doublet forthe airmass (left) and molecfit (right) telluric removal methods.By comparison, we see that corrections using airmass gives a trans-mission spectrum with more variation, particularly in regions whereresidual telluric lines were removed (see Figure 7). The Gaussianfits to the full resolution data around the Na lines have reducedchi-square values of 𝜒 𝜈 = .
54 for the airmass correction methodand 𝜒 𝜈 = .
50 for the molecfit method.For each method, we co-added the nightly spectra to ensurethat the signal-to-noise ratio is large enough for a well characterisedatmospheric sodium detection. Weaker signals would also be moreprominent. The final CLV and RM corrected transmission spectraare shown in Figure 9. The Gaussian fits have 𝜒 𝜈 = .
58 (airmass)
MNRAS000
MNRAS000 , 1–12 (2021)
A. Langeveld et al. and 𝜒 𝜈 = .
57 ( molecfit ). Given the comparable fits with bothapproaches, no conclusion can be drawn about which method isbetter based on this statistic alone. The measured parameters of theGaussian fits to the the sodium D2 and D1 lines are shown in Table 2.Following reduction using the airmass telluric correction method,we measured line contrasts of − . ± .
07 % for the D2 line and − . ± .
07 % for the D1 line. This results in an average depthof − . ± .
05 %. For the molecfit method, the line contrastswere − . ± .
07 % (D2) and − . ± .
07 % (D1), averagingto − . ± .
05 %. Although the average line contrasts for bothmethods are not within 1 𝜎 of each other, the error ranges overlap.From this we note that telluric corrections using the molecfit method resulted in a transmission spectrum with deeper Na featuresthan those for the airmass method. We find a similar difference whenmeasuring binned absorption depths, and further context is providedin section 5.2. The measured line contrasts and differences betweenthe D2 and D1 lines are comparable to the results of Wyttenbachet al. (2015) (D2: − . ± .
07 %, D1: − . ± .
07 %) andCasasayas-Barris et al. (2017) (D2: − . ± .
05 %, D1: − . ± .
05 %).
HD 189733 b has been studied extensively in the past at differentresolutions (Redfield et al. 2008; Jensen et al. 2011; Huitson et al.2012; Wyttenbach et al. 2015; Khalafinejad et al. 2017). In orderto form a comparison to these (and studies of other planets), it isbeneficial to calculate the strength of the absorption over varying binwidths centred at the sodium lines (Snellen et al. 2008; Wyttenbachet al. 2015). This is particularly useful if the line itself is too narrowto be resolved, or if there is too much noise to produce an accuratefit. We define absorption depth (AD) as the difference betweenthe mean flux of the continuum regions and the mean flux withina passband centred on the spectral feature. Red and blue contin-uum passbands are chosen in regions around the sodium doublet:5874.89–5886.89 Å for the blue band, and 5898.89–5910.89 Å forthe red band. The mean flux in central passbands of widths 0.375,0.75, 1.5, 3, 6 and 12 Å are measured. Since the sodium featurecontains two lines, the total widths of these passbands are dividedinto two, with one half centred around each line. However, the 12 Åpassband is large enough to encapsulate both features, so this isnot split and is centred in the middle of the two lines instead. Theabsorption depth is thus calculated byAD = (cid:60) (cid:48) 𝐶 − (cid:60) (cid:48) 𝐵 + (cid:60) (cid:48) 𝑅 , (8)where (cid:60) (cid:48) is the mean flux of the transmission spectrum in the blue(B), red (R) and central (C) passbands. This quantity was evaluatedover all passbands for each nightly and combined transmission spec-trum. All averages were weighted using the inverse of the squareduncertainties (variance), which were propagated from the photonnoise of the raw data. The results are shown in Table 3 (airmass) andTable 4 ( molecfit ). The measured absorption depths of the com-bined data are shallower for the airmass method than the molecfit method. This matches the difference seen in the line contrasts insection 5.1.Narrower passbands more accurately probe the sodium linecores, so these measurements will most closely match the Gaussianline contrasts. For the molecfit reduction, the absorption depthswere − . ± .
066 % and − . ± .
046 % for the 0.375 and 0.75 Å passbands respectively. These values are consistent to within1 𝜎 of each other, and the error for the 0.75 Å passband overlaps withthe error in the average Gaussian line contrast. The same comparisonis true for the airmass method.We compare our results to those of Wyttenbach et al. (2015),who analyse the same HARPS data and compare the effects of plane-tary radial velocity corrections. In the 1.5 Å passband, they measurean absorption depth of − . ± .
031 %, which agrees to within1 𝜎 with our results of − . ± .
031 % for the airmass method and − . ± .
031 % for the molecfit method. We therefore confirmthe 10 𝜎 detection of sodium in the atmosphere of HD 189733 b.Table 3 shows slight variation to the results discussed by Wyt-tenbach et al. (2015), despite similarities in the telluric removalmethods. This is likely due to differences in the reduction pipelineand handling of the data, as well as corrections for the CLV and RMeffects which change the shape and depth of the absorption lines.A consistent trend is seen in both studies: the measured depthsfor night 2 using the airmass reduction are significantly lower thanthe other two nights. As discussed previously, the observations onthis night began when the planet was already in-transit, so a longout-of-transit baseline was not obtained. A clear telluric referencespectrum was only able to be derived when the in-transit spectrawere included. Whilst this is better for telluric correction, it couldalso lead to partial removal of signals in the transmission spectrum.We were unable to detect any significant absorption in the D1 linefor this night, which gives rise to the lower overall measurements.The molecfit reduction performs better in this case. As seen inTable 4, the measured depths are much more consistent across allnights for each passband, and there is a significant improvement inthe results for night 2. Therefore, molecfit is much better at re-moving telluric contamination in spectra, regardless of the numberof in or out-of-transit observations. It is also noted that the errorsacross single passbands are comparable for both methods, thereforethe signal-to-noise ratio is improved when using molecfit .As the passbands increase in size, the measured depths forboth telluric correction methods decrease by the same proportion.We measure an absorption depth of − . ± .
007 % in the largest(12 Å) passband using molecfit . This agrees with the space-baseddetection of − . ± .
006 % (Huitson et al. 2012), and two otherground-based detections of − . ± .
020 % (Redfield et al. 2008)and − . ± .
017 % (Jensen et al. 2011) within the same 12 Åpassband. The corresponding measurement for the airmass methodwas − . ± .
008 %, which does not agree as significantly. Thissuggests that molecfit is better at correcting tellurics in the im-mediate regions around the sodium doublet.
Sodium absorption features in the lab rest-frame should be locatedat wavelengths of 5889.951 Å for the D2 line and 5895.924 Åfor the D1 line. Despite efforts made to correct for spectral devia-tions due to systemic and planetary radial velocities, the measuredcentroids of the Gaussian fits to the full resolution data (Table 2)have a net blueshift from the rest-frame Na doublet positions (al-though still encapsulated within error). We interpret this as beingdue to winds within the atmosphere of HD 189733 b, in agree-ment with other studies (Wyttenbach et al. 2015; Louden & Wheat-ley 2015; Casasayas-Barris et al. 2017). The average blueshift was − . ± .
029 Å for the airmass method and − . ± .
025 Å forthe molecfit method. This corresponds to net wind velocities of − . ± . − and − . ± . − respectively, indicatingan eastward atmospheric movement from day-side to night-side. MNRAS , 1–12 (2021) elluric correction of exoplanet transmission spectra -1.00.0 Airmass-fit Corrected Molecfit Corrected
Wavelength ( Å ) -0.50.00.5 5880.0 5890.0 5900.0 5910.0 Wavelength ( Å ) ' ( % ) R e s i d . ( % ) Figure 8.
Na i doublet absorption in the atmosphere of HD 189733 b. The top row shows the transmission spectrum for night 3 (in the planetary rest-frame)using the airmass (left) and molecfit (right) telluric correction methods. Full resolution data are shown in grey, and binned 20x in black. The bottom rowshows the residuals to the best-fit Gaussian profiles (red) in the same resolution. The expected rest-frame positions of the D-lines are indicated with the bluedashed lines. -1.00.0
Airmass-fit Corrected Molecfit Corrected
Wavelength ( Å ) -0.50.00.5 5880.0 5890.0 5900.0 5910.0 Wavelength ( Å ) ' ( % ) R e s i d . ( % ) Figure 9.
Same as Figure 8 but for the combined data across all three nights. Stacking multiple nights of data in the rest-frame of the planet increases thesignal-to-noise ratio, as is evident by visual comparison with the transmission spectrum for night 3.Night AD . 𝑎 AD . 𝑎 AD . 𝑎 AD 𝑎 AD 𝑎 AD 𝑎 − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . Table 3.
Calculated absorption depth (%) of the Na doublet for each transmission spectrum following the airmass ( 𝑎 ) telluric removal method. The flux isaveraged across six passbands with total widths of 0.375, 0.75, 1.5, 3, 6 and 12 Å. Combined depths are measured after co-adding individual nightly spectra.Night AD . 𝑚 AD . 𝑚 AD . 𝑚 AD 𝑚 AD 𝑚 AD 𝑚 − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . − . ± . Table 4.
Same as Table 3 but measuring from the transmission spectrum after removing telluric contamination with molecfit ( 𝑚 ).MNRAS000
Same as Table 3 but measuring from the transmission spectrum after removing telluric contamination with molecfit ( 𝑚 ).MNRAS000 , 1–12 (2021) A. Langeveld et al.
Since the telluric corrections are the only parts that vary be-tween our methods, we deduce that the difference between the ve-locities may be due to inaccuracies in telluric removal – particularlyin nights 2 and 3 for the airmass reduction. Nevertheless, both re-sults are lower than the − ± − measurement discussedby Wyttenbach et al. (2015). The discrepancy can be attributed tonuances in the reduction pipeline and radial velocity corrections,and application of models to correct for CLV and RM effects. Ifthe spectra are observed uniformly over the full transit and stackedin the stellar rest frame, the RM effect can be averaged out. Asshown in Figure 6, this is not true when the transmission spectra arestacked in the planetary rest-frame – spurious signals are induced(albeit small for the slowly rotating host star) and were removed bymodelling these effects.Our measured wind velocities are consistent with the resultsof Louden & Wheatley (2015) ( − . + . − . km s − ) and Casasayas-Barris et al. (2017) ( − − ) who analysed the same data. Ad-ditionally, Louden & Wheatley (2015) spatially resolved the at-mospheric winds and measured a redshift on the leading limb ofthe planet (2 . + . − . km s − ) and a blueshift on the trailing limb( − . + . − . km s − ). This implies a net eastward motion due to acombination of tidally-locked planetary rotation and strong east-ward winds. High-resolution observations in the infrared regimehave also reported similar day-side to night-side wind speeds: − . + . − . km s − (Brogi et al. 2016) and − . + . − . km s − (Brogiet al. 2018). Additionally, our measured wind velocities are consis-tent with predictions from three-dimensional General CirculationModels (GCMs), where the atmospheric flow is dominated by asuperrotating eastward equatorial jet and the hottest region of theatmosphere is advected downwind from the substellar point (Show-man et al. 2009; Miller-Ricci Kempton & Rauscher 2012; Rauscher& Kempton 2014). For HD 189733 b, models predict that the equa-torial jet extends to latitudes of ∼ ◦ . Atmospheric winds floweastward within this region – from night-to-day at the western termi-nator (leading limb) and from day-to-night at the eastern terminator(trailing limb). At latitudes above ∼ ◦ , the winds flow preferen-tially from day-to-night at both the leading and trailing limbs ofthe terminator. Therefore, models predict a net blueshift due to thewinds (Showman et al. 2013). An exoplanet’s transit depth is wavelength-dependent; if there ismolecular or atomic absorption, the atmosphere will appear opaquewhen viewed at wavelengths within this region. The planet willhave a larger apparent radius at these wavelengths and block morestellar flux, compared to a transparent atmosphere with no absorp-tion. Therefore, we can use absorption measurements to probe theatmospheric height at which the molecules or atoms are present.We define atmospheric scale height 𝐻 sc as the increase in altitudeover which the pressure reduces by a factor of 𝑒 , given by 𝐻 sc = 𝑘 𝑏 𝑇 eq 𝜇 𝑚 𝑔 . (9)Using a mean molecular weight of 𝜇 𝑚 = .
22 u (for Jupiter) andderiving the planet’s surface gravity of 𝑔 = . − from the pa-rameters in Table A1, the scale height for HD 189733 b is 200.7 km.We refer to the white-light transit depth as the decrease in stellarflux during transit, equal to the ratio of the sky-projected areas of the planet and star: Δ = (cid:18) 𝑅 𝑝 𝑅 ∗ (cid:19) . (10)For HD 189733 b, Δ = . 𝐻 𝜆 , so a wavelength-dependent transit depth is similarly defined as Δ 𝜆 = (cid:18) 𝑅 𝑝 + 𝐻 𝜆 𝑅 ∗ (cid:19) = (cid:18) 𝑅 𝑝 𝑅 ∗ (cid:19) + 𝑅 𝑝 𝐻 𝜆 𝑅 ∗ + (cid:18) 𝐻 𝜆 𝑅 ∗ (cid:19) . (11)Since the atmosphere typically extends around 5–10 scale heights(Madhusudhan 2019), 𝐻 𝜆 << 𝑅 ∗ , so Δ 𝜆 = Δ + 𝑅 𝑝 𝐻 𝜆 𝑅 ∗ . (12)For a particular wavelength, the difference between the white-lightand wavelength-dependent transit depths is a measure of the ab-sorption, 𝛿 𝜆 = Δ 𝜆 − Δ = 𝑅 𝑝 𝐻 𝜆 𝑅 ∗ = Δ 𝐻 𝜆 𝑅 𝑝 , (13)where 𝛿 𝜆 = − 𝐹 𝜆, in 𝐹 𝜆, out = −(cid:60) (cid:48) 𝜆 . (14)Using the line contrasts and passband absorption depths as a mea-surement of (cid:60) (cid:48) 𝜆 , we can therefore calculate the height of an atmo-spheric detection by rearranging equation 13 for 𝐻 𝜆 . For the molec-fit reduction, the average Na doublet line contrast of − . ± .
05 %corresponds to an apparent planetary radius of 1 . 𝑅 𝑝 and an at-mospheric height of 𝐻 Na = ±
900 km. This is ∼
51 timesgreater than the scale height and suggests that the detection of Naprobes the upper atmosphere of HD 189733 b.
HD 189733 b is one of the most extensively studied exoplanets todate. A number of chemical species have been detected throughobservations from space and ground-based instruments, as well asmeasurements of wind speeds, atmospheric circulation, and evi-dence of Rayleigh scattering and high-altitude hazes. The natureof telluric contamination makes ground-based observations of ex-oplanet transmission spectra difficult. Several methods have beenused in the past to remove excess absorption from Earth’s atmo-sphere, including observing a telluric standard star (Redfield et al.2008; Jensen et al. 2011, 2012), computing an empirical telluricspectrum (Wyttenbach et al. 2015), and using an Earth atmosphericmodel (Yan et al. 2015; Casasayas-Barris et al. 2017; Allart et al.2017).In this study we approached the challenge of removing telluriclines in high-resolution optical spectra by comparing two popu-lar methods using archival HARPS observations of HD 189733 b.Firstly, we derived a telluric reference spectrum from the data byassuming a linear relationship between airmass and telluric linestrength. The reference spectrum was used to scale the strength oftelluric lines to a constant level as if they had been observed at thesame airmass, and are thus removed during the subsequent trans-mission spectrum calculation. Our second method involved using molecfit to generate a model of molecular H O and O absorp-tion in Earth’s atmosphere. It is then straightforward to remove thecontamination through division of unique telluric contributions for MNRAS , 1–12 (2021) elluric correction of exoplanet transmission spectra each observation. The methods can be applied to all regions of theHARPS spectra, however we focused specifically on the qualityof the telluric correction around the sodium D-lines and the im-pact on the measured line properties. Other empirical methods suchas SYSREM or principal component analysis (Tamuz et al. 2005;Mazeh et al. 2007; Birkby et al. 2013), which are normally reservedfor severe telluric contamination in the near-infrared, are beyond thescope of this study.Alongside correcting for stellar, systemic and planetary ra-dial velocity effects, we modelled the Centre-to-Limb Variation andRossiter-McLaughlin effect to correct for spurious signals which areimprinted onto the transmission spectrum. The relative strengths ofthese signals were small and at the noise level of the spectrum, how-ever it is an important factor to consider if the star is rotating morerapidly. Both methods were performed on the data from each nightseparately, then combined for an overall transmission spectrum.Absorption due to atmospheric sodium (Na i) was identifiedfrom the doublet lines at ∼ ∼ − . ± .
07 % (D2) and − . ± .
07 % (D1), and absorptionin the 0.75 Å passband of − . ± .
046 %. For the molecfit reduction, the corresponding values were − . ± .
07 % (D2 linecontrast), − . ± .
07 % (D1 line contrast), and − . ± .
046 %(0.75 Å passband absorption). In cases where weather conditionsare variable over one night, the model produced with molecfit was more robust at removing telluric contamination, and absorp-tion depth measurements were more consistent across all nights.Additionally, it performs better when the observing window doesnot allow for a long out-of-transit baseline to be obtained (whenbuilding an empirical telluric reference spectrum is difficult). Theaverage line contrast of the Na doublet for the molecfit methodwas − . ± .
05 %, corresponding to an atmospheric height of10200 ±
900 km. Finally, we measured a net atmospheric blueshiftof − . ± . − arising from a combination of tidally-lockedplanetary rotation and day-side to night-side winds.From the assessments of the results, we conclude that correct-ing telluric contamination with a model of molecular (H O and O )absorption in Earth’s atmosphere using molecfit produces a bet-ter quality transmission spectrum than if the tellurics were removedempirically, e.g. using airmass. This is evident by: (a) a detection ofatmospheric sodium with greater significance, (b) more consistentmeasurements of absorption depths in bandpasses across differentnights, and (c) the ability to extract a signal when the weather condi-tions are not ideal. Molecfit allows for targeted removal of specifictelluric lines without significantly impacting data outside of theseregions. This is a key advantage over data-driven methods such asthe airmass correction used within this work, since it is likely thatthe empirically derived model will more significantly increase thecontinuum noise. Additionally, corrections can be performed on asingle spectrum without relying on a large sample of observations,which is particularly advantageous if only a partial transit is ob-tained or if variable weather conditions affect a large number ofspectra.Detections at different spectral resolutions and across a num-ber of wavelength regimes are important for thorough atmosphericcharacterisation. With an ever-growing number of identifications ofchemical species within exoplanet atmospheres, we are able to fur-ther understand and constrain dynamical and chemical process andmake inferences into how hot Jupiters form and evolve. We havealready seen how ground-based instruments such as HARPS can bea powerful tool for such studies in the optical regime, and this isset to increase with newly developed high-resolution spectrographs such as ESPRESSO (Pepe et al. 2013), HARPS-3 (Young et al.2018) and EXPRES (Jurgenson et al. 2016). ESPRESSO covers awavelength range of 380–780 nm which probes deeper into the redwavelengths, giving opportunity to make detections of the waterband at ∼ ∼ ACKNOWLEDGEMENTS
We thank the anonymous referee for their thoughtful commentswhich helped to improve the quality of this manuscript. AL ac-knowledges support from the Science and Technology FacilitiesCouncil (STFC), UK. We thank Annelies Mortier for insightfuldiscussion about the HARPS spectrograph and data.
DATA AVAILABILITY
All the data used in this work are available through the ESO sci-ence archive facility. The data are based on observations collectedat the European Southern Observatory under ESO programmes072.C-0488(E), 079.C-0828(A) (PI: Mayor) and 079.C-0127(A)(PI: Lecavelier des Etangs).
REFERENCES
Addison B., et al., 2019, PASP, 131, 115003Agol E., Cowan N. B., Knutson H. A., Deming D., Steffen J. H., HenryG. W., Charbonneau D., 2010, ApJ, 721, 1861Allart R., Lovis C., Pino L., Wyttenbach A., Ehrenreich D., Pepe F., 2017,A&A, 606, A144Allart R., et al., 2018, Science, 362, 1384Allart R., et al., 2019, A&A, 623, A58Alonso-Floriano F. J., et al., 2019, A&A, 629, A110Astudillo-Defru N., Rojo P., 2013, A&A, 557, A56Birkby J. L., de Kok R. J., Brogi M., de Mooij E. J. W., Schwarz H., AlbrechtS., Snellen I. A. G., 2013, MNRAS, 436, L35Boisse I., et al., 2009, in Pont F., Sasselov D., Holman M. J.,eds, IAU Symposium Vol. 253, Transiting Planets. pp 462–465,doi:10.1017/S174392130802694XBonomo A. S., et al., 2017, A&A, 602, A107Bourrier V., et al., 2013, A&A, 551, A63Brogi M., Snellen I. A. G., de Kok R. J., Albrecht S., Birkby J., de Mooij E.J. W., 2012, Nature, 486, 502Brogi M., de Kok R. J., Albrecht S., Snellen I. A. G., Birkby J. L., SchwarzH., 2016, ApJ, 817, 106Brogi M., Giacobbe P., Guilluy G., de Kok R. J., Sozzetti A., Mancini L.,Bonomo A. S., 2018, A&A, 615, A16Cabot S. H. C., Madhusudhan N., Hawker G. A., Gandhi S., 2019, MNRAS,482, 4422Cabot S. H. C., Madhusudhan N., Welbanks L., Piette A., Gandhi S., 2020,arXiv e-prints, p. arXiv:2001.07196Casasayas-Barris N., Palle E., Nowak G., Yan F., Nortmann L., Murgas F.,2017, A&A, 608, A135Casasayas-Barris N., et al., 2018, A&A, 616, A151Casasayas-Barris N., et al., 2019, A&A, 628, A9Casasayas-Barris N., et al., 2020, arXiv e-prints, p. arXiv:2002.10595Cauley P. W., Redfield S., Jensen A. G., 2017, AJ, 153, 217Cauley P. W., Shkolnik E. L., Ilyin I., Strassmeier K. G., Redfield S., JensenA., 2019, AJ, 157, 69Cegla H. M., Lovis C., Bourrier V., Beeck B., Watson C. A., Pepe F., 2016,A&A, 588, A127MNRAS000
Addison B., et al., 2019, PASP, 131, 115003Agol E., Cowan N. B., Knutson H. A., Deming D., Steffen J. H., HenryG. W., Charbonneau D., 2010, ApJ, 721, 1861Allart R., Lovis C., Pino L., Wyttenbach A., Ehrenreich D., Pepe F., 2017,A&A, 606, A144Allart R., et al., 2018, Science, 362, 1384Allart R., et al., 2019, A&A, 623, A58Alonso-Floriano F. J., et al., 2019, A&A, 629, A110Astudillo-Defru N., Rojo P., 2013, A&A, 557, A56Birkby J. L., de Kok R. J., Brogi M., de Mooij E. J. W., Schwarz H., AlbrechtS., Snellen I. A. G., 2013, MNRAS, 436, L35Boisse I., et al., 2009, in Pont F., Sasselov D., Holman M. J.,eds, IAU Symposium Vol. 253, Transiting Planets. pp 462–465,doi:10.1017/S174392130802694XBonomo A. S., et al., 2017, A&A, 602, A107Bourrier V., et al., 2013, A&A, 551, A63Brogi M., Snellen I. A. G., de Kok R. J., Albrecht S., Birkby J., de Mooij E.J. W., 2012, Nature, 486, 502Brogi M., de Kok R. J., Albrecht S., Snellen I. A. G., Birkby J. L., SchwarzH., 2016, ApJ, 817, 106Brogi M., Giacobbe P., Guilluy G., de Kok R. J., Sozzetti A., Mancini L.,Bonomo A. S., 2018, A&A, 615, A16Cabot S. H. C., Madhusudhan N., Hawker G. A., Gandhi S., 2019, MNRAS,482, 4422Cabot S. H. C., Madhusudhan N., Welbanks L., Piette A., Gandhi S., 2020,arXiv e-prints, p. arXiv:2001.07196Casasayas-Barris N., Palle E., Nowak G., Yan F., Nortmann L., Murgas F.,2017, A&A, 608, A135Casasayas-Barris N., et al., 2018, A&A, 616, A151Casasayas-Barris N., et al., 2019, A&A, 628, A9Casasayas-Barris N., et al., 2020, arXiv e-prints, p. arXiv:2002.10595Cauley P. W., Redfield S., Jensen A. G., 2017, AJ, 153, 217Cauley P. W., Shkolnik E. L., Ilyin I., Strassmeier K. G., Redfield S., JensenA., 2019, AJ, 157, 69Cegla H. M., Lovis C., Bourrier V., Beeck B., Watson C. A., Pepe F., 2016,A&A, 588, A127MNRAS000 , 1–12 (2021) A. Langeveld et al.
Charbonneau D., Brown T. M., Noyes R. W., Gilliland R. L., 2002, ApJ,568, 377Chen G., et al., 2018, A&A, 616, A145Chen G., Casasayas-Barris N., Palle E., Yan F., Stangret M., Cegla H. M.,Allart R., Lovis C., 2020, arXiv e-prints, p. arXiv:2002.08379Czesla S., Klocová T., Khalafinejad S., Wolter U., Schmitt J. H. M. M.,2015, A&A, 582, A51Deibert E. K., de Mooij E. J. W., Jayawardhana R., Fortney J. J., Brogi M.,Rustamkulov Z., Tamura M., 2019, AJ, 157, 58Di Gloria E., Snellen I. A. G., Albrecht S., 2015, A&A, 580, A84Fischer D. A., Howard A. W., Laughlin G. P., Macintosh B., MahadevanS., Sahlmann J., Yee J. C., 2014, in Beuther H., Klessen R. S.,Dullemond C. P., Henning T., eds, Protostars and Planets VI. p. 715( arXiv:1505.06869 ), doi:10.2458/azu_uapress_9780816531240-ch031Gandhi S., Madhusudhan N., 2017, MNRAS, 472, 2334Heng K., Wyttenbach A., Lavie B., Sing D. K., Ehrenreich D., Lovis C.,2015, ApJ, 803, L9Hoeijmakers H. J., et al., 2019, A&A, 627, A165Huitson C. M., Sing D. K., Vidal-Madjar A., Ballester G. E., Lecavelier desEtangs A., Désert J. M., Pont F., 2012, MNRAS, 422, 2477Jensen A. G., Redfield S., Endl M., Cochran W. D., Koesterke L., BarmanT. S., 2011, ApJ, 743, 203Jensen A. G., Redfield S., Endl M., Cochran W. D., Koesterke L., BarmanT., 2012, ApJ, 751, 86Jensen A. G., Cauley P. W., Redfield S., Cochran W. D., Endl M., 2018, AJ,156, 154Jurgenson C., Fischer D., McCracken T., Sawyer D., Szymkowiak A., DavisA., Muller G., Santoro F., 2016, EXPRES: a next generation RV spec-trograph in the search for earth-like worlds. Proceedings of the SPIE, p.99086T, doi:10.1117/12.2233002Kausch W., et al., 2015, A&A, 576, A78Khalafinejad S., et al., 2017, A&A, 598, A131Kirk J., Alam M. K., López-Morales M., Zeng L., 2020, AJ, 159, 115Lecavelier Des Etangs A., Pont F., Vidal-Madjar A., Sing D., 2008, A&A,481, L83Lecavelier Des Etangs A., et al., 2010, A&A, 514, A72Lockwood A. C., Johnson J. A., Bender C. F., Carr J. S., Barman T., RichertA. J. W., Blake G. A., 2014, ApJ, 783, L29Louden T., Wheatley P. J., 2015, ApJ, 814, L24Madhusudhan N., 2019, ARA&A, 57, 617Mandel K., Agol E., 2002, ApJ, 580, L171Mayor M., et al., 2003, The Messenger, 114, 20Mazeh T., Tamuz O., Zucker S., 2007, in Afonso C., Weldrake D., HenningT., eds, Astronomical Society of the Pacific Conference Series Vol.366, Transiting Extrapolar Planets Workshop. p. 119 ( arXiv:astro-ph/0612418 )McCullough P. R., Crouzet N., Deming D., Madhusudhan N., 2014, ApJ,791, 55McLaughlin D. B., 1924, ApJ, 60, 22Miller-Ricci Kempton E., Rauscher E., 2012, ApJ, 751, 117Nikolov N., Sing D. K., Gibson N. P., Fortney J. J., Evans T. M., BarstowJ. K., Kataria T., Wilson P. A., 2016, ApJ, 832, 191Nortmann L., et al., 2018, Science, 362, 1388Pepe F., et al., 2013, The Messenger, 153, 6Pinhas A., Madhusudhan N., Gandhi S., MacDonald R., 2019, MNRAS,482, 1485Piskorz D., et al., 2016, ApJ, 832, 131Pont F., Knutson H., Gilliland R. L., Moutou C., Charbonneau D., 2008,MNRAS, 385, 109Queloz D., Eggenberger A., Mayor M., Perrier C., Beuzit J. L., Naef D.,Sivan J. P., Udry S., 2000, A&A, 359, L13Rauscher E., Kempton E. M. R., 2014, ApJ, 790, 79Redfield S., Endl M., Cochran W. D., Koesterke L., 2008, ApJ, 673, L87Rossiter R. A., 1924, ApJ, 60, 15Ryabchikova T., Piskunov N., Kurucz R. L., Stempels H. C., Heiter U.,Pakhomov Y., Barklem P. S., 2015, Phys. Scr., 90, 054005Salz M., et al., 2018, A&A, 620, A97 Seager S., Sasselov D. D., 2000, ApJ, 537, 916Sedaghati E., et al., 2016, A&A, 596, A47Seidel J. V., et al., 2019, A&A, 623, A166Showman A. P., Fortney J. J., Lian Y., Marley M. S., Freedman R. S.,Knutson H. A., Charbonneau D., 2009, ApJ, 699, 564Showman A. P., Fortney J. J., Lewis N. K., Shabram M., 2013, ApJ, 762, 24Sing D. K., et al., 2016, Nature, 529, 59Smette A., et al., 2015, A&A, 576, A77Snellen I. A. G., Albrecht S., de Mooij E. J. W., Le Poole R. S., 2008, A&A,487, 357Snellen I. A. G., de Kok R. J., de Mooij E. J. W., Albrecht S., 2010, Nature,465, 1049Stassun K. G., Collins K. A., Gaudi B. S., 2017, AJ, 153, 136Stevenson K. B., et al., 2014, Science, 346, 838Tamuz O., Mazeh T., Zucker S., 2005, MNRAS, 356, 1466Torres G., Winn J. N., Holman M. J., 2008, ApJ, 677, 1324Triaud A. H. M. J., 2018, The Rossiter-McLaughlin Effect in ExoplanetResearch. Springer International Publishing AG, p. 2, doi:10.1007/978-3-319-55333-7_2Triaud A. H. M. J., et al., 2009, A&A, 506, 377Ulmer-Moll S., Figueira P., Neal J. J., Santos N. C., Bonnefoy M., 2019,A&A, 621, A79Valenti J. A., Piskunov N., 1996, A&AS, 118, 595Vidal-Madjar A., et al., 2010, A&A, 523, A57Welbanks L., Madhusudhan N., Allard N. F., Hubeny I., Spiegelman F.,Leininger T., 2019, ApJ, 887, L20Wyttenbach A., Ehrenreich D., Lovis C., Udry S., Pepe F., 2015, A&A, 577,A62Wyttenbach A., et al., 2017, A&A, 602, A36Yan F., Fosbury R. A. E., Petr-Gotzens M. G., Zhao G., Wang W., Wang L.,Liu Y., Pallé E., 2015, International Journal of Astrobiology, 14, 255Yan F., Pallé E., Fosbury R. A. E., Petr-Gotzens M. G., Henning T., 2017,A&A, 603, A73Young J., Naylor T., Brake M., Fisher M., Seneta E., Kunst P., Pepe F., Sos-nowska D., 2018, in Proc. SPIE. p. 107072J, doi:10.1117/12.2312830de Kok R. J., Brogi M., Snellen I. A. G., Birkby J., Albrecht S., de MooijE. J. W., 2013, A&A, 554, A82
APPENDIX A: SYSTEM PARAMETERS
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS , 1–12 (2021) elluric correction of exoplanet transmission spectra Parameter Symbol Value Unit Reference
Star
Stellar Mass M ∗ + . − . M (cid:12) Triaud et al. (2009)Stellar Radius R ∗ + . − . R (cid:12) Torres et al. (2008)Stellar RV Semi-amplitude K ∗ + . − . m s − Boisse et al. (2009)Effective Temperature T eff + − K Stassun et al. (2017)Projected Rotational Velocity 𝑣 sin 𝑖 + . − . km s − Bonomo et al. (2017)Surface Gravity log 𝑔 + . − . log (cm s − ) Stassun et al. (2017)Metallicity [Fe/H] -0.03 + . − . dex Torres et al. (2008) Planet
Planetary Mass M 𝑝 + . − . M 𝐽 Triaud et al. (2009)Planetary Radius R 𝑝 + . − . R 𝐽 Torres et al. (2008)Planetary RV Semi-amplitude K 𝑝 + − km s − DerivedEquilibrium Temperature T eq + − K Addison et al. (2019)Orbital Inclination 𝑖 𝑝 + . − . deg. Agol et al. (2010)Sky Projected Obliquity 𝜆 -0.4 + . − . deg. Cegla et al. (2016) System
Period P 2.21857567 + . − . days Agol et al. (2010)Mid-transit Time T + . − . BJD
UTC
Agol et al. (2010)Transit Duration 𝜏 + . − . days Triaud et al. (2009)Semi-major axis a 0.03120 + . − . A.U. Triaud et al. (2009)Systemic Velocity 𝑣 sys + . − . km s − Boisse et al. (2009)
Table A1.
Stellar, planetary and orbital parameters of the HD 189733 system. K 𝑝 is derived through the relationship − 𝐾 ∗ ( 𝑀 ∗ / 𝑀 𝑝 ) .MNRAS000