TiO WASP-33b draft

High-Resolution Spectroscopic Detection of TiO and Stratosphere in the Day-side of WASP-33b


We report high-resolution spectroscopic detection of TiO molecular signature in the day-side spectra of WASP-33 b, the second hottest known hot Jupiter. We used High-Dispersion Spectrograph (HDS; R 165,000) in the wavelength range of 0.62 – 0.88 m with the Subaru telescope to obtain the day-side spectra of WASP-33 b. We suppress and correct the systematic effects of the instrument, the telluric and stellar lines by using SYSREM algorithm after the selection of good orders based on Barnard star and other M-type stars. We detect a 4.8- signal at an orbital velocity of = +237.5 km s and systemic velocity = -1.5 km s, which agree with the derived values from the previous analysis of primary transit. Our detection with the temperature inversion model implies the existence of stratosphere in its atmosphere, however, we were unable to constrain the volume-mixing ratio of the detected TiO. We also measure the stellar radial velocity and use it to obtain a more stringent constraint on the orbital velocity, km s. Our results demonstrate that high-dispersion spectroscopy is a powerful tool to characterize the atmosphere of an exoplanet, even in the optical wavelength range, and show a promising potential in using and developing similar techniques with high-dispersion spectrograph on current 10m-class and future extremely large telescopes.

techniques: spectroscopic - planets and satellites: atmospheres - planets and satellites: composition - planets and satellites: individual (WASP-33 b)

Stevanus K. Nugroho


Indonesia Endowment Fund for Education Scholar \affiliationDepartment of Earth and Planetary Science, The University of Tokyo, Tokyo 113-0033, Japan


Department of Earth and Planetary Science, The University of Tokyo, Tokyo 113-0033, Japan \affiliationResearch Center for the Early Universe, School of Science, The University of Tokyo, Tokyo 113-0033, Japan


Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, U.S.A \affiliationNASA Sagan Fellow


Department of Earth and Planetary Sciences, Tokyo Institute of Technology, Tokyo 152-8551, Japan


National Astronomical Observatory of Japan, Tokyo 181-8588, Japan \affiliationAstrobiology Center, National Institutes of Natural Sciences, Tokyo, Japan


Subaru Telescope, 650 N. A’ohoku Place, Hilo, HI 96720, U.S.A

1 Introduction

Thermal inversion in exoplanetary atmosphere is still a problem not completely understood in the theory of exoplanetary atmosphere. This inversion layer was predicted by Hubeny et al. (2003) and Fortney et al. (2008) in a highly irradiated planet, which is caused by high-temperature absorber molecules such as TiO or VO that absorb UV and visible radiation from the incoming stellar radiation, and heat up the upper atmosphere. By measuring shallower eclipse depths at 4.5 m and 5.6 m, caused mainly by HO and CO, than the thermal continuum at the nearby band (3.6 and 9 m) using the Spitzer Space Telescope (Spitzer), the first evidence of inversion layer was reported by Knutson et al. (2008) in the atmosphere of HD 209458 b. Although there have been similar observations of secondary eclipse depth by the same telescope suggesting the detection of an inversion layer, it was shown by Hansen et al. (2014) that several previous single-transit observations using Spitzer have significantly higher uncertainties than previously known, which made the measured emission-like features doubtful. Using the CRyogenic high-resolution InfraRed Echelle SPectrograph (CRIRES Kaeufl et al., 2004) in the Very Large Telescope (VLT), Schwarz et al. (2015) observed the day-side of HD 209458 b and reported the non-detection of CO at 2.3 m and constrained the nonexistence of the inversion layer in the pressure range of 10–10 bar, which is consistent with the results by Zellem et al. (2014) and Diamond-Lowe et al. (2014). Despite the non-detection of CO in the day-side of HD 209458 b, Snellen et al. (2010) detected a CO absorption feature at 5.6- using a similar instrument in the transmission spectrum. One of the possible causes of the non-detection is that the average day-side atmosphere of HD 209458 b is near-isothermal in the pressure range that they probe, which makes the day-side spectrum almost featureless (Schwarz et al., 2015).

Evans et al. (2017) reported the detection of HO emission features in a super-hot Jupiter, WASP-121b (T 2500 K), using the combination of Hubble Space Telescope (HST) and Spitzer. The presence of water was resolved and detected at 5- confidence level, strengthening the previous conclusion by Evans et al. (2016), which made WASP-121b the first exoplanet with resolved emission features and detected stratosphere layer. Their findings also include the possible detection of VO in its atmosphere at the level of 1000x solar abundance, and a thermal inversion-like atmospheric structure, even though they only used a 1D atmospheric model structure and did not take the non-equilibrium chemistry into account. Sedaghati et al. (2017) observed WASP-19b during transit using a low-dispersion spectrograph, FORS2, on VLT. They confirmed the presence of water (7.9-) and simultaneously revealed the presence of TiO (7.7-), strongly scattering haze (7.4-) and sodium (3.4-). They also constrained the relative abundance of those molecules, however, the presence of a thermal inversion remains unproved because a transmission spectrum has less information on temperature structure of the atmosphere. Other evidence of thermal inversion has also been reported in the atmosphere of WASP-33 b (T 2700 K), which is the second hottest hot Jupiter, by Haynes et al. (2015), who observed the day-side of the exoplanet using HST, and revealed a thermal excess at about 1.2 m, which is consistent with the TiO spectral feature. Other than these, there has been no significant evidence, nor direct detection of TiO and/or VO in the atmosphere of hot Jupiters.

The non-detection of TiO or VO can be caused by several factors. Hubeny et al. (2003) and Spiegel et al. (2009) suggested that gravitational settling could drag TiO/VO from the upper atmosphere to the colder layer in the deeper atmosphere. Meanwhile, owing to high-speed winds of the tidally locked hot Jupiter, condensed molecules can also be brought into the colder night side of the planet. If the temperature of the atmosphere is below the TiO/VO condensation level, it can condensate and if the vertical mixing rate, which is related to the temperature of the planet atmosphere is not high enough, it cannot be redistributed to the upper atmosphere. This effect is called the cold-trap effect. Knutson et al. (2010) found a possible connection between the UV chromospheric stellar activity and the existence of an inversion layer in the atmosphere of known hot Jupiters. Hot Jupiters orbiting active stars tend to have no inversion layer and those orbiting quiet stars showed evidences of the existence of an inversion layer in their atmosphere. Knutson et al. (2010) suggested that the increased UV intensity can probably destroy the compounds that are responsible of creating an inversion layer. Meanwhile, using high-dispersion spectroscopy, Hoeijmakers et al. (2015) reported the inaccuracy of the TiO line list that was used at wavelengths shorter than 6300 . The accuracy level, however, tends to increase at longer wavelengths, which was shown by the cross-correlation result between the spectrum of Barnard’s Star and the model spectrum created by using the corresponding TiO line list (see Figure 9 of Hoeijmakers et al., 2015).

Recently, direct detection of molecular signature in exoplanet atmosphere using high-dispersion spectroscopy is one of the most widespread approach to the attempt of exoplanet characterization (e.g Snellen et al., 2010; Crossfield et al., 2011; Birkby et al., 2013; de Kok et al., 2013; Hoeijmakers et al., 2015; Schwarz et al., 2015; Birkby et al., 2017; Esteves et al., 2017). Unlike low-dispersion spectroscopy, high-dispersion spectroscopy can resolve molecular bands into individual absorption lines. The variation of Doppler shifts during observation caused by the orbital movement enables to distinguish absorption lines in the exoplanet spectrum from telluric lines, and ensures the unambiguous detection of specific molecules. Owing to these resolved individual lines it is also possible to investigate several physical parameters of the exoplanet such as the axial tilt (Kawahara, 2012), the projected equatorial rotational velocity (Snellen et al., 2014), wind speed (Brogi et al., 2016), and the thermal inversion layer (Schwarz et al., 2015) of the planet. Obtaining the exoplanet spectrum can be done by cross-correlating the data with the exoplanet atmosphere spectrum model, after removing telluric and stellar lines using a specific method.

We observed WASP-33 b (Smith et al., 2011), which is the second hottest known hot Jupiter (wavelength dependent brightness temperature of 3620 K), orbiting quiet Scuti stars. It has a retrograde orbit with a period of 1.22 days. It is an ideal choice for transmission spectroscopy measurements: owing to its unusually large radius (Collier Cameron et al., 2010) and its high temperature making it not only suitable for secondary eclipse spectroscopy measurements, but also as the main target to find TiO/VO in its atmosphere, as the cold-trap effect is unlikely at this temperature level.

In this paper, we report a direct detection of TiO molecules and a stratosphere layer of WASP-33b, based on ground-based observation of the day-side emission spectrum in the visible wavelength range (6170-8817). The observation, data reduction, and systematic effect removal (including the correction of the blaze function variation, common wavelength grid, and telluric and stellar line removal) are described in §2. We also examine our methods to create the model spectrum and to cross-correlate the data with the model spectrum. Here, we also explain the method to confirm the radial velocity of WASP-33 and the accuracy of TiO line list used by us. In §3, the possibility of signal detection of TiO and the order-by-order (order-based) optimization of the SYSREM algorithm are explored. Here, we also show the final result after the optimization and the statistical test. This is followed by the discussion and the conclusions our findings in §4.

2 Observation and Data Reduction

2.1 Subaru Observation of WASP 33

We observed WASP-33 on UT October 26 2015 (see Table 1 for stellar and exoplanet physical and dynamical parameters) using High Dispersion Spectrograph (HDS Noguchi et al., 2002) at f/12.71 optical Nasmyth focus of the the Subaru 8.2 m telescope (proposal ID: S15B-090, PI: H. Kawahara). The observation was conducted in a standard NIRc set up (without Iodine cell) with Messia5, 2x1 binning setting, and without image rotator. Image slicer 3 (Tajitsu et al., 2012), each with slit width= 0.2 was used, resulting in the highest spectral resolution of R= 165,000 (1.8 km s resolution), which were sampled by two detectors (blue and red CCDs) with 4100 2048 pixels (0.9 km s per pixel) containing 18 orders covering 6170 – 7402 , and 12 orders covering 7537 – 8817 , respectively. We obtained 52 spectra of WASP 33, each with an exposure time of 600 s, from air mass 2.97 to 1.96 on the other meridian covering 0.23 – 0.56 exoplanet orbital phases (see Table 1 for ephemeris and derived the orbital period, the orbital coverage of our observation is shown by bold red line in Figure 1) with a typical seeing of 0.6 – 0.7.

Figure 1: The coverage of WASP-33 b orbital phase during our observation (showed by bold red line). The vertical dotted line shows the ingress and egress phase
Parameter Value
() 1.5611
1.512 0.04 2
1.495 0.031 3
() 1.509
() 6.17 0.43
Spectral type A5
(K) 7400 2004
4.3 0.2
Fe/H 0.1 0.25
M/H 0.1 0.26
Centre-of-mass velocity (km s) -2.19 0.097
-3.69 0.098
-2.11 0.059
sin (km s) 86.6310
() 117 2
-2450000 (BJD) 6934.77146 0.0005911
(days) 1.219870912
(days) 0.1143 0.000213
(days) 0.0124 0.000214
3.69 0.0115
i () 88.69516
2015 () 3.266 0.72617
() 1.67918
CGS 3.4619
(km s) 231.11 20
228.67 21
227.81 22
Table 1: System parameters of WASP-33

2.2 M-dwarfs Spectra

Hoeijmakers et al. (2015) showed that the TiO line list used by them is not accurate at short wavelengths ( 6000 ) however, it tends to be more accurate at longer wavelengths. For a robustness of analysis, we checked the accuracy of the TiO line list we use in this paper, by comparing it with the spectra of Barnard’s Star (M4V) and HD 95735 (M1.5V), as TiO is commonly found in the atmosphere of M-type stars and brown dwarfs (Burrows & Sharp, 1999; Burrows et al., 2001). In 2016, we observed Gl752A (M3V), and HD173739 (M3V) with the Subaru telescope using the same instrumental configuration as in 2015 (proposal ID: S16A-107, PI: H. Kawahara). To check the robustness of our cross-correlation we also downloaded and used the calibrated spectrum of Proxima Centauri from the ESO Science Archive Facility. The spectrum was obtained by XSHOOTER at VLT between 5336.6 and 10200 with a resolution of R= 18,000 (proposal ID: 092.D-0300(A), PI: Neves).


f2a.png0.48(a) \figf2b.png0.48(b)

Figure 2: Relative shifts of telluric lines compared to those in the first frame in the blue CCD (1) and red CCD (2). Noticeably the general trend of the shift follows the trend of temperature variation inside the dome (solid red line)

2.3 Standard Reduction

The data were reduced using IRAF tools23 and a custom-built application written in Python 2.7. We corrected the over-scan and non-linearity using CL scripts that can be obtained from the HDS website24 before debiasing the frames. During the analysis of Narita et al. (2005) data, Snellen et al. (2008) noticed the non-linearity effect in HDS. After fixing this problem empirically, the analysis gave the detection of sodium in the atmosphere of HD 209458b. The CCD response (gain) for high signal-to-noise (S/N) spectra is higher than those in a lower S/N spectra, which would make the correction of telluric and stellar lines difficult. This problem was solved by Tajitsu et al. (2010), who provided a CL script to correct this effect after applying the over-scan correction. The scattered light in the inter-order area was fitted with a cubic spline function along the slit and dispersion direction individually, using apscatter for the arc lamp frame, smoothed and subtracted from all science frames. A median flat frame was calculated from 71 dome-flat frames to correct the pixel-to-pixel sensitivity variation, and normalized using the apnormalize task in IRAF in order to conserve the fringe pattern along the slit direction. Using the CL script hdsis_ecf.cl 25, we flat-fielded and extracted 1-dimensional spectra of total 30 orders. Using hdsis_ecf.cl the science spectra were sliced along the slit from 12 to 8 pixels and from 11 to 8 pixels relative to the center of the order tracer of each order for the blue and red CCD respectively, divided them by the slice of the normalized median flat frame and sum combined. If a larger width was used, a high-frequency noise would appear along the wavelength in the extracted spectra, caused by the edge of the aperture, which was defined earlier in apnormalize by using a median flat frame. We derived the wavelength solution by identifying 133 emission lines of Thorium-Argon arc lamp frames taken at the beginning and the end of the observation, and fitted with a fifth and third order Chebyshev function corresponding to the dispersion and slit directions respectively, using ecidentify. The pixel RMS value of the residual fitting was 0.0012 for both CCDs. Then, the spectra were assigned by interpolating between the preceding and the following Thorium-Argon arc lamp frames in respect of the observation time, using refspec with the relative distance of the spectra to the reference spectra as the weight of the interpolation. Then, the wavelength solution was applied using dispcor. These steps were applied to the non-normalized median flat frame to create the corrector for the blaze function. All science frames were then divided by the blaze function corrector to create the final reduced frames. Then the next step is correcting the variation of blaze function along the time which is given in detailed in Appendix A.

2.4 Common Wavelength Grid

We measured the shifts of the spectrum during the observation, by measuring the relative shifts of strong telluric lines to the reference spectrum. First, we detected telluric lines using peakutils26 (15 – 20 lines on average per selected order) in the order that contain well-spread telluric lines in the wavelength direction. We then fitted a Gaussian function to determine the precise centroid of each line and compared it with the centroid of the same lines in the first frame. The relative shift of each order was calculated by taking the median combination of the shifts in that order. There are order dependent shifts of about 0.05 pixel, but because the information for the order that does not has strong telluric lines is not available and the value is too small to affect the result (0.05 pixels corresponding to 0.04 km s), we decided to ignore it. The results of all selected order were combined to estimate the shift of the corresponding frame, by calculating its weighted median combined.

The amplitude of relative shifts is about 0.13 and 0.08 pixel (0.12 and 0.07 km s) for blue and red CCD, respectively. Noticeably the general trend of the shifts are correlated with the shifts of the temperature inside the dome (T 1.75 K) during the observation (see Figure 2 a b). This wavelength-temperature shift relation trend is relatively weaker than those for the CRIRES in VLT reported by Brogi et al. (2013) (1.5 K change in temperatures corresponding to 1.5 pixels shifts), which gave a non-detection of the 51 Peg b signal due to the odd-even effect for detector 4.

Then, the spectra were shifted with spline interpolation, based on the weighted median combined shift of each frame, and then re-sampled into a common wavelength bin. The spectra of each order were then stacked in a matrix with wavelength bin values in its column, and spectrum numbers (or time/orbital phases) in the rows. 5- clipping was performed for each wavelength bin to identify any bad pixels/cosmic rays, which are replaced by the mean value of the bin. We excluded these bad region for the rest of the analysis, which is 5.30% and 7.19% of all pixels in blue and red CCDs, respectively, mostly because they are on the edge of the CCD.


f3a.png1.(a) \gridline\figf3b.png1.(b)

Figure 3: Reduction process for each order of both blue (a) and red CCDs (b). The first row ([1]) shows an example of the normalized 1D spectrum of WASP-33 following the correction of the blaze function variation and the common wavelength grid iteration. The next row ([2]) shows the 2D spectrum with the wavelength as the horizontal axis, each label representing the median wavelength of the order, while the vertical axis is the frame number. The third row ([3]) shows the final reduced spectra after the correction of the blaze function in common wavelength grid. The variation of brightness along the frames for all orders is due to the blaze function variation. The fourth row ([4]) shows the mean subtracted spectra as the input to SYSREM. The fifth and sixth rows ([5] [6]) show the residual spectra after running SYSREM with 1 and 4 iteration and through the double high-pass filter, at the latter stage almost all telluric lines have been removed. The masked bad regions are shown by grey areas in all rows except the first one.

2.5 Removal of Telluric and Stellar Absorption Lines

In the studied wavelength range, telluric and stellar absorption lines are dominating the spectrum, while the exoplanet to stellar flux contrast is expected at the level of 10 (when extrapolated from the best-fitted spectrum in Haynes et al., 2015). Removal of telluric and stellar lines is critical in order to detect the TiO signature by cross-correlation, as the exoplanet spectrum is buried under sharp noises due to those lines. The the spectrum of WASP-33b is expected to be Doppler-shifted from +230 km s (at its maximum) to 54 km s (at the end of the observation), while telluric and stellar lines remain relatively stationary during the course of the observation. Telluric lines, which are dominated mostly by water and oxygen lines, are varied in strength due to the changes of geometric air mass and water vapor column of the atmosphere.

The quasi-static telluric and stellar absorption lines were removed by implementing the SYSREM algorithm to remove systematic effects (variation of atmospheric condition, the changing of CCD efficiency, variation of point spread function, etc.) without any priors in a large set of light curves (Tamuz et al., 2005; Mazeh et al., 2007). It has been used either in the detrending/systematic effects removal of transit surveys (e.g SuperWASP, CoRoT light curves Collier Cameron et al., 2006; Ofir et al., 2010), or in the removal of quasi-static telluric lines in the similar analysis to this work (e.g. Birkby et al., 2013, 2017; Esteves et al., 2017). Each wavelength bin (4100 wavelength bin per order on average) is treated as a “light-curve” consisting of 52 frames (including the frames when WASP-33 b in the secondary eclipse phase). Then each “light-curve” was subtracted from their mean value before applying the SYSREM algorithm order-by-order (fourth row in Figure 3). The detailed description of the SYSREM algorithm is given in Appendix B.

The step-by-step removal of telluric can be seen in Figure 3. It can be noticed that there are curve shaped features, which can be seen in the blue CCD spectra matrix, and the degree of the curvature decreases as increasing wavelength. These features can be caused by the imperfection of the correction of the blaze function variation. As the residual was cross-correlated with the high-frequency TiO features, this feature does not affect the results and these were roughly removed after a double high-pass filter was applied (see Figure 3).

Any linear systematic variation along all wavelength bins in each order can be found and removed from the spectrum starting with the most significant one, such as air mass variation. Ideally, the expected residual is the Doppler-shifted exoplanet spectrum only, with additional noise. We perform two cases of the SYSREM reduction: One is that we remove the first ten systematics (henceforth SYSREM iteration, ) for all of the orders (SYSREM with common iteration number in §3.1). This simple procedure gives a conservative estimate of the signal detection. Another attempt is that we determine the optimum number of subsequent SYSREM iterations for each order that gives the optimum systematic removal (Order-based SYSREM optimization in §3.2). The latter is based on the fact that SYSREM tends to remove the exoplanet spectrum after the Doppler-shifted variation dominate the systematic change in the spectrum.

Before performing cross-correlation for each residual, we applied a double high-pass filter similarly to the case of M-dwarfs star spectrum, using a smoothing function with a 25-pixel width for the first filter and 51-pixel width for the second filter, to remove any low-frequency variations along the spectrum, before cross-correlating with a grid of the Doppler-shifted WASP-33b model spectrum. Then, the values of each wavelength bin were weighted by their noise, which is defined by the standard deviation of the bin as a function of time. The results of this removal process are 10 sets of smoothed weighted SYSREM residuals for each CCD.

2.6 Spectral Template

To constrain the existence of TiO molecules in WASP-33 b emission spectra by cross-correlation, model spectra was generated with various profiles and volume mixing ratio (VMR) of TiO, assuming several atmospheric models as described in Table 2. Three different temperature–pressure (T/P) profiles were adopted, which described the average vertical temperature structure of the planetary day-side atmosphere: full inversion (FI), non-inversion (NI) as shown in the right panels of Figure 4, and T/P profile from Haynes et al. (2015), known as the H-model, as shown in the right panel of Figure 5. For the NI model, it was assumed that the temperature decreases with a constant lapse rate from P= 10 bar with T= 3700 K to P= 10 bar with T= 2700 K. For the FI model, the temperature increases from P= 10 with T= 2700 K to P= 10 with T= 3700 K. For both models, a constant TiO concentration at solar abundance level was assumed, that is, VMR= and no other molecule in the atmosphere. For the H-model, seven constant VMRs of TiO from sub-solar to super-solar abundance level were assumed, VMR= [4, 5, 6, 7, 8, 9, 10]. The atmosphere was divided into 50 layers, which were evenly spaced in log pressure between 10 to 10 bar, and the altitudes were calculated by assuming hydro-static H dominant atmosphere and using the derived mass and a radius of WASP-33b from the reference (Table 1). We also modeled the non-inverted atmosphere (henceforth M-dwarf model) with T= 2600 K at P= 0.01 bar and T= 4000 K at P= 1 bar assuming a constant rate of temperature change with pressure, resulting in TiO absorption features only for VMR= 7.

Name T/P Profile log VMR
FI Full inversion 7
NI No inversion 7
H Realistic inversion (Haynes et al. (2015)) 4 – 10
Table 2: Models of planetary atmosphere

The cross section of molecules was calculated using Python scripts for Computational Atmospheric Spectroscopy (Py4CATS27). The procedure of computing the cross section is described in Appendix C. Then, the cross-sections were combined into continuum absorption coefficients and integrated along the line-of-sight through the atmosphere, resulting delta optical depths per layer. The Schwarzschild equation was solved by integrating the Planck function versus the monochromatic transmission, along the line-of-sight and convolved with a Gaussian kernel to match with the HDS resolution. In total, we produced 9 WASP-33b mock spectra and 1 non-inverted model for TiO line list accuracy analysis; 1 full inversion (FI-spec), 1 no inversion (NI-spec), 7 spectra for Haynes (H-spec), and 1 M-dwarf model. See Figure 4 for FI-spec and NI-spec, Figure 5 for H-spec, and Figure 6 for 5 M-dwarfs spectra1 M-dwarf model spectrum.

In this analysis, the systemic velocity (V) is one of the important parameters to confirm the detection (if there is any); thus by measuring radial-velocity (RV) of WASP-33 and comparing it with the previous results, the confidence and the robustness of our analysis are improved. Instead of creating the WASP-33 comparison spectra by our self, we used the stellar model spectrum from Coelho (2014) for T= 7500 K, g= 4.5, Fe/H= 0.2 and /Fe= 0. The stellar model spectrum was convolved to the HDS resolution and rotationally broadened to resemble the Doppler broadening caused by its fast rotation (henceforth Coelho model spectrum).

Figure 4: WASP-33b spectrum model using full inversion (FI) and no inversion (NI) of the TP profile. The left panel shows the WASP-33b model spectrum for both FI-spec and NI-spec. The color shade represents the wavelength range of the observed order, blue, and orange color represent the wavelength range in the blue and red CCDs, respectively. The vertical axis shows the planet to star flux contrast in the level of 10. The middle panel shows the weight function along the wavelength for the pressure range considered. The right panel shows our adopted temperature profile of the WASP-33b atmosphere.
Figure 5: WASP-33b spectrum model using TP profile of Haynes et al. (2015). The left panel shows the WASP-33b model spectrum of various VMR. The vertical axis shows the planet to star flux contrast in the level of 10. The right panel shows our adopted temperature profile of the WASP-33b atmosphere.
Figure 6: Observed spectrum of five M-dwarf stars that were used for the TiO line list accuracy comparison. The bumps on the 4 first spectra are caused by the uncorrected blaze function. Before order-by-order cross-correlation, the spectrum was normalized by its continuum profile. The color shade represents the wavelength range of the observed order, the blue color is in the blue CCD and orange color is in the red CCD.

f7a.png0.45(a) \figf7b.png0.45(b)

Figure 7: (a) The colored line is an example of blue CCD spectrum of WASP-33 after standard reduction. Different colors mean different order. The black line is Coelho model spectrum. (b) The barycentric corrected RV of WASP-33 during the observation (red star), as a result of combining the RV of each order that has the highest cross-correlation signal. The blue line is the weighted median combine of all frames that was taken as the RV of WASP-33.The error bar is the 1- scatter across the observation time.

2.7 System Velocity from Stellar Spectra

To measure the radial velocity of the WASP-33 system, we analyzed the standard reduced WASP-33 spectra and masked the telluric lines of the selected order that contain significant stellar absorption lines. Then, the spectra were cross-correlated order-by-order with the Doppler-shifted Coelho model spectrum from 100 km s to 100 km s with 0.1 km s intervals. In this cross-correlation, only the spectra in blue CCD was used, because there are a lot of telluric lines and less stellar absorption lines in the red CCD spectrum. For each selected order, the RV value with the highest cross-correlation signal was extracted, and median combined to calculate the RV of WASP-33 for the corresponding frame. Then, the weighted median combine of the RV of all frames was calculated, following the correction of the barycentric radial velocity difference to have the final WASP-33 RV value.

The results can be seen in Figure 7b, and the median of the measured RV is 3.02 0.42 km s. Although this value is different from the value 9.2 2.8 km s (Gontcharov, 2006) in the SIMBAD database, our result is consistent with Collier Cameron et al. (2010), which used Doppler tomography technique to confirm the existence of the planet and measured the center-of-mass velocity (see Table 1). Therefore we used this value as the radial velocity of the WASP-33 system, which is taken into consideration when evaluating the result of the WASP-33b vs TiO model cross-correlation analysis.

Line List Accuracy and Excluding Bad Orders

Figure 8: Cross-correlation results between M-dwarf model and five M-dwarf spectra, Barnard’s Star (black line), HD 95735 (blue line), Gl752A (green line), HD173739 (red line), and Proxima Centauri (brown line). The color of the label of the wavelength range also represents the order to which the CCD belongs. Note that the y-axis scale is not uniform, thus the orders cannot be compared. The black dash line is the rest-frame RV of each M-dwarf stars, while the vertical blue dashed line in the first three orders of the blue CCD is the position of CCF peaks on those orders. The blue and red shaded panel show the bad shape CCF that are masked for the rest of the analysis (bad). The grey shaded panel shows the masked order in the later analysis due to heavily contaminated by strong telluric lines (ht).

To check the TiO line list accuracy, five standard reduced M-dwarf spectra were cross-correlated order-by-order with the M-dwarf model. The model spectra were Doppler-shifted from 80 km s to 80 km s relative to the RV of the expected target, from SIMBAD database with 1 km s intervals. The object and model spectrum was divided by their continuum profile, calculated by applying a double high-pass filter with 501 pixels and 1001 pixels of smoothing function before the cross-correlation, in order to maintain the cross-correlation scale (from 1 to 1).

The cross-correlation results of five M-dwarf spectra versus the M-dwarf model is shown in Figure 8. As Hoeijmakers et al. (2015) expected, the accuracy of the TiO line list improved in the longer wavelength, although there are cross-correlation functions (CCFs) of several order, which blue-shifted from their expected radial velocity, and several others do not have any significant peak. For the other of the order, the peaks are located at their expected radial-velocity. The shifts are most likely caused by the inaccuracy of the TiO line list itself, because the first three orders of the blue CCD show similar blue-shifted CCF for the five different M-dwarf spectra. The shift is about 21 km s (shown by blue dashed line in Figure 8). There are several orders that have no significant peak, which can be caused by the imperfection of our simple atmosphere modeling or the inaccuracy of the line list. Note that all orders that have no significant CCF peak, have a large bad region except the last order of the red CCD (see Figure 3 and/or Figure 8). These bad orders are excluded for the rest of the analysis including order 2 in red CCD due to heavily contaminated by strong telluric lines.

3 TiO Signal Detection

3.1 Results from SYSREM with Common Iteration Number

As described in §2.5, we first analyze the spectra with a common SYSREM iteration number for all orders. A grid of Doppler-shifted WASP-33b model spectrum was cross-correlated with weighted SYSREM residuals. The Doppler-shifted model spectrum covers planet radial velocity (RV) between 169.69 km s RV 393.30 km s with 0.5 km s intervals, corresponding to half of the HDS sampling resolution. For each detector, the CCFs of every good orders (all orders excluding the bad orders explained in §2.7.1) were summed. This value is listed in the CCF matrix with a dimension of 1127 (RV as column) 52 (orbital phase as row). Then, the CCF matrix of blue and red CCDs was summed to calculate the final CCF matrix. These steps were done for all SYSREM residuals.

We then calculate the CCF map in the plane for all spectrum models by doing the following steps. The CCF of the frames (40 frames in total, excluding the frames when WASP-33b was expected in the secondary eclipse phase) was integrated along the expected RV curve:




where is the barycentric correction, is the planet orbital phase, is the mid observation time in BJD, is the ephemeris, and is the orbital period of the planet, is the semi-amplitude of the radial velocity of the planet, and is the systemic velocity 28. We consider the radial velocity of the planet semi-amplitude between 150 km s 310 km s, and systemic velocity between 80 km s 80 km s with 0.5 km s steps. The result is a 321 ( as column) 320 ( as row) CCF matrix. Then this matrix was divided by its the standard deviation to make S/N map.

Figure 9 shows the S/N map for (H-spec) and 10 (H-, FI-, and NI-spec models). In all the cases except for NI-spec, we detected positive peaks with 4-, while the map for NI-spec exhibits a negative peak at the same place. The maps for H-spec and FI-spec exhibits the strong peak in almost all SYSREM iteration number at = 237.0 km s and = 1.5 km s (henceforth peak A). Among all H-spec, the strongest signal was found in log VMR= 8. This peak is also the strongest one in the c.c. maps of both FI-spec and NI-spec for most of the SYSREM iteration number, although for NI-spec the value is negative. The positive value of the peak A in H-spec and FI-spec, and the negative value in NI-spec suggested that non-inversion atmosphere is unlikely for WASP-33b if peak A is the real signal. The second strongest peak was detected at = 192.0 km s and = 19.5 km s (henceforth peak B) in the SYSREM iteration number= 2 only.


f9a.png0.48(a) \figf9b.png0.48(b) \gridline\figf9c.png0.48(c) \figf9d.png0.48(d)

Figure 9: (a) and (b) are the cross-correlation map of H-spec (log VMR= -8) for SYSREM iteration number= 4 and 10 respectively. (c) and (d) are the cross-correlation map of FI-spec and NI-spec for SYSREM iteration number= 10 respectively. The white dashed line is the expected Kp and Vsys from the previous studies. The black line is the maximum S/N peak in the map. The top and right panels show the 1-dimensional cross section of the CCF peak along the and respectively. A and B are the peak A and peak B respectively.

Figure 10 shows the S/N of peaks A and B as a function of the SYSREM iteration, although peak B is unlikely to be a real signal, due to its physically non-realistic and values compared with the expected value from previous studies, it was chosen as a representative of the noise/false-positive signal. The S/N of peak B decreases as increasing of the SYSREM iteration number. The S/N of peak A increases from SYSREM iteration= 1 until SYSREM iteration= 4 then decreases until SYSREM iteration= 6 before begin to increase again until SYSREM iteration= 10. The increasing of S/N of peak A after SYSREM iteration= 6 is most likely due to the difference level of telluric and/or stellar lines removal in each order. It can be seen in the Figure 9 a for SYSREM iteration= 4 where multiple peaks, 5 peaks, within 1 value from the strongest peak can be found (white color). While in the Figure 9 b there are only 2 peaks, peak A and the other one at about 160.0 km s and 20 km s. From the curve of Figure 10, it is natural to adopt as a fiducial iteration number. Because the common iteration number for all orders does not fully optimize the procedure of the systematics removal, 4.3- of the significance level for H-spec (VMR= 8) we obtained should be regarded as a conservative estimate of detection level.

Figure 10: The S/N of peak A and B along the SYSREM iteration number, , for H-spec (log VMR= 8). The blue diamond-lined is the S/N of peak A and the red circle-dotted line is the S/N of peak B.

3.2 Order-based SYSREM Optimization

The level of systematics (telluric lines) should be different for each order, thus to find the optimized SYSREM iteration number of each order, following steps were also performed. A scaled artificial signal was injected at the detected RV (= 237 km s and = 1.5 km s) in the spectra before telluric and stellar lines removal. The scaling was performed according to the equation


where is the scaling constant, is the scaled artificial signal, is the planet model spectrum (H-spec with VMR= 8) from the integration of Planck function versus monochromatic transmission along the line-of-sight (see §2.6), is the black body flux of T= 7400 K representing the continuum level of the WASP-33 flux, and are the planet and star radius, respectively.

To check the dependence on the strength of the injected signal, we adopt 5 different injected signals with [0.2, 0.4, 0.6, 0.8, 1.0]. The injected signals are broadened by a rotation kernel, using fastRotBroad from PyAstronomy with = 0.4 times the expected projected velocity of WASP-33b (calculated by assuming a tidally locked condition with the parameter in the Table 1 referring to the typical broadening width of tidally locked exoplanet as shown in Kawahara, 2012). We convolve with a box-function to take into account the change of Doppler shift during the 600 s exposures per frame. The signals were injected into the spectra, before performing telluric and/or stellar lines removal using SYSREM algorithm. The residuals for each SYSREM iteration were cross-correlated with the artificial signal itself.

The CCF of each order of each frame was aligned to the planet rest-frame according to the injected (e.g. the middle panel of Figure 13). Then, the aligned CCFs of all frames except the frames corresponding to the secondary eclipse were summed to create the total mean CCF. To examine the exoplanet signal for various combinations of SYSREM iteration numbers, the peak value of the total mean CCF in 2.5 km s was taken between the expected rest-frame radial velocity of the (0 km s) planet and divided by the its outside standard deviation (non-detection CCF value) to include the “noise” (henceforth ). The results for five injected signals are shown in Figure 11.

Figure 11: of all orders for various SYSREM iteration numbers. The vertical lines mark the SYSREM iteration number that made the signal strength optimum. The blue circle, yellow asterisk, green diamond, red square, and purple cross are the of the recovered injected signal with = [0.2, 0.4, 0.6, 0.8, 1.0] respectively. The number in the white box at the bottom left of each panel represents the order number and the color represents the CCD (blue or red).

The signal strength () behaves differently in every order, depending on the amount of telluric and/or stellar lines contamination and the strength of the injected signal (represented by value). There are several jumps in curves, which may be caused by spurious signal, thus by following the general trend of the recovered curves of 5 different injected signals the optimum SYSREM iteration number we chose the optimum SYSREM iteration numbers. For blue and red CCDs the optimum SYSREM iteration numbers are [2, 2, 8, 2, 1, 1, 5, 3, 1, 3, 5, 2, 7, 2, 8, 2, 2, 2] and [2, 3, 5, 2, 3, 4, 2, 8, 3, 9, 3, 5], respectively, for each of its orders. Figure 12 shows the level of telluric and/or stellar lines contamination of each order, which is represented by the standard deviation value of the standard reduced spectrum in that order, and its corresponding optimum SYSREM iteration number. One can see the tendency that the optimum SYSREM iteration number increases as exhibiting more telluric and/or stellar lines in the spectrum, except for one in the bottom right of the figure ( the order having many O lines).

Figure 12: Standard deviation of the spectrum of each order with their corresponding optimized SYSREM iteration number. The blue cross and the red circle indicate the spectrum in the blue and red CCDs, respectively. The value of the standard deviation represents the amount of telluric and/or stellar lines contamination for every order.

Using the optimum SYSREM iteration number, the mean CCF map was calculated for all of the spectrum model as shown in Figure 13 (left panel). Then the CCF was aligned to the rest-frame of the planet (middle panel). The TiO signal can be seen in the left panels as a positive (dark) signal with an arc shape for H- and FI- spec and as a negative (bright) one for NI-spec. In the middle panels of the figures, the planet signal was aligned such that it can be seen as a vertical dark/bright trail at the = 1.5 km s. The right panels show the mean CCF with the width of 6 pixels centered at = 1.5 km s (CCF). For FI- and H-models, the CCF curve exhibits a positive offset during the visible phase while there is no offset when the planet is behind the star. This feature supports the atmospheric origin of the CCF signal. The exposure time for a single frame is 600 s, which corresponds to 7 km s at the near-eclipse phase. This long exposure time should have smeared out the day-side spectrum by the change in the radial component of the orbital velocity of the planet, especially in the phases near the secondary eclipse. This smearing effect may account for the fact that the signal is only visible at the phase 0.35, while it should be stronger at the phase 0.35 because larger part of the dayside of the WASP-33b is visible.

Figure 13: Order-based optimized version of cross-correlation map for H-spec with log VMR= 8 (first row), FI-spec (second row), NI-spec (third row), and the injected signal (fourth row). The left panel shows the mean CCF map. The middle panel shows the aligned mean CCF map at = 237.5 km s and = 1.5 km s. The planet signal can be seen as a dark trail in the expected rest-frame of the radial velocity of the planet from the first observed phase until the appearance of the secondary eclipse phase (red dashed). The right panel shows the mean CCF for a 6 pixel column bin (CCF) centered on = 1.5 km s along the orbital phase. The solid black line shows the smoothed CCF with a 3 frame smoothing window. The blue dashed line shows the zero value of CCF.
Figure 14: Figure S/N map of H-spec with maximum peak at = 237.5 km s and = 1.5 km s, which gives 4.8- detection. The white dashed line is the expected Kp and Vsys from the previous studies. The top and right panels show the 1-dimensional cross section of the CCF peak along the and respectively. The black dash lines show the most significant signal, the white dash lines show the expected and , and the color bar grid interval is 1-, the white area also represents the 1- error of the detected signal.

Figure 14 shows a S/N map for H-spec with VMR= 8. The noise level in the S/N map after optimization is significantly suppressed. The strongest peak was found with 4.8- detection significance at = 237.5 km s and = 1.5 km s within an elliptic region of 1 sigma. The latter is consistent with the estimates from the stellar RV 2.7 and Collier Cameron et al., 2010). Because we already have the strong constraint on km s from the stellar spectrum (§2.7), we may assume this value as the prior of on the S/N map. Then we obtain a stronger constraint, km s. This is the first time that the orbital velocity of WASP-33b was dynamically measured. The orbital velocity of a planet depends on the orbital period and the mass of the host star. Since the measurement of orbital period is very precise for a transiting exoplanet, our observation provides the first model-independent measurement of the mass of the host star. Using the period of WASP-33b from the previous study (see Table 1) the measured mass of the host star is , which is larger than values previously reported.

3.3 Statistical Tests

Figure 15: Significance map after converting from the -value from Welch’s t-test using function for the most significant detected signal with the H-spec model spectrum of VMR= 8. The black dashes show the most significant signal, the white dashes show the expected and , and the color bar grid interval is 1-, the white area also represents the 1- error of the detected signal.
Figure 16: Histogram of the in-trail and out-trail mean CCF distribution, the in-trail distribution is slightly shifted from the out-of-trail distribution.
Figure 17: Q-Q plot of the out-of-trail distribution which shows that the distribution is Gaussian until 4-.

The in-trail CCF (henceforth in-trail signal) was compared with the out-of-trail CCF (henceforth out-of-trail signal) by performing Welch’s t-test to check if the mean is same assuming that two distributions were drawn from the same parent distribution, using the SciPy module in Python 2.7. We used equation (1) to calculate the expected RV for the same range to ones for the CCF maps, and took the 1 pixel mean CCF value at the closest RV to RV. The out-of-trail signal contains all mean CCF values except the in-trail signal.

The out-of-trail and in-trail signal histograms (width of the in-trail signal= 3 pixels) were plotted. The in-trail signal distribution shifted further from the out-of-trail signal distribution. The distribution of the in-trail signal is clearly shifted from the zero values. From the Q-Q plot (see Figure 17) it can be seen that the out-of-trail signal distribution is a Gaussian distribution until 4-, therefore we can safely convert the half -value to value of the detection significance using an error function29. The Welch’s t test shows that the in-trail distribution is deviated from the out-of-trail signal distribution by 5.0- (see Figure 15) in line with the S/N of this peak detection.

4 Discussion and Conclusions

We confirmed the previously claimed of the inaccuracy of the TiO line list for wavelengths shorter than 6300 (see Hoeijmakers et al., 2015). We also showed that the line list accuracy is enough high for longer wavelengths, which was considered in processing order-based optimization, therefore our analysis no longer suffers from the inaccuracy issue. By measuring the radial velocity of the host star we also confirmed the measurement by Collier Cameron et al. (2010), which differs by 4 km s from the SIMBAD database. This confirmation gives us a narrow search space to find the possible exoplanet signal.

We reported a TiO molecule signature detection in the day-side spectra of WASP-33b by 4.8- confidence level, which provided direct evidence of the existence of TiO in the atmosphere of the hot Jupiter. The detection levels for VMR= -8, -9, and -10 (H-spec) lie within 1- from the highest one (VMR= -8); moreover, in our analysis we did not use a self-consistent atmospheric model, thus the constraint on the VMR cannot be obtained directly from our result only. Our TiO molecular detection for H-spec and FI-spec confirmed the existence of stratosphere (thermal inversion layer) in the day-side of WASP-33 b, as previously has been claimed by several studies (e.g. von Essen et al., 2015; Haynes et al., 2015). The full inversion T/P profile has also been reported for another super hot Jupiter, WASP-121b by Evans et al. (2017), which resolved the emission spectral feature of HO at near-infrared wavelength, using HST. Our result is complementary to the TiO detection using low-dispersion spectroscopy in WASP-19b by Sedaghati et al. (2017), who were able to constrain the relative abundance of TiO (0.12 p.p.b) but could not provide information about T/P profile as they were only able to measure the transmission level of the molecules in the atmosphere. The constraint on the relative abundance level of TiO in WASP-33b can be obtained by analyzing it using a self-consistent atmospheric model and/or by combining low- and high-dispersion spectroscopy in order to introduce a more precise constraint of the relative abundance of each detected molecules and the T/P profile of the exoplanet atmosphere (Brogi et al., 2017).

By observing for about nine hours only using the 8.2 m Subaru telescope, we are able to detect significant signature of TiO in the atmosphere of WASP-33b. Our results demonstrate that high-dispersion spectroscopy is a powerful tool to characterize the atmosphere of an exoplanet, and show a promising potential of developing similar/more advanced techniques for the Infra Red Doppler instrument (Kotani et al., 2014) in the Subaru telescope, and for extra large telescopes facilities in the future, as suggested by several authors (e.g. Kawahara & Hirano, 2014; Snellen et al., 2015).

5 Acknowledgement

This work was based in part on data collected on HDS at the Subaru Telescope, which is operated by the National Astronomical Observatory of Japan. We thank our anonymous referee, whose insightful comments improved the manuscript. S.K.N. acknowledges support from Indonesia Endowment Fund for Education Scholarship. S.K.N also acknowledges Toru Yamada for fruitful discussions during the analysis of the data. H.K. is supported by Grant-in-Aid for Young Scientist (B) from Japan Society for Promotion of Science (JSPS), No. 17K14246 and the Astrobiology Center from NINS, No. AB291003. T.H. is supported by Grant-in-Aid for Young Scientist (B) from Japan Society for Promotion of Science (JSPS), No. 16K17660.

Appendix A Correction for Blaze Function Variation and Normalization of Spectra

As Winn et al. (2004) showed, observation using HDS suffered from variation of the blaze function during the observation. The continuum profile is not important in our analysis, but it is useful to correct this variation as a method of subtracting the continuum of each spectrum. We used the first frame as a reference spectrum and compared the rest of spectra by calculating the ratio of both spectra. By looking at the ratio of all spectra, we found wavelength variations with similar patterns for all orders, but these vary along the time of the observation (see Figure 18). The ratio spectra were then clipped to remove any outlier. In order to avoid any residual broad stellar absorption lines mismatch in the ratio spectra, we selected a 301-pixels smoothing window and applied smoothing, using PyAstronomy.pyasl.smooth30 with a flat window function. Each of the compared spectra was then divided by its smoothed ratio spectra, resulting in a spectra with shared blaze function with the reference spectrum, while conserving stellar lines and telluric line strength variations. A spline function was fitted to the reference spectrum using continuum manually and divided all spectra by it to get the normalized spectra.

Figure 18: Top panel: Spectra from two different epoch for order 93 (left panel), 92 (middle panel), 91 (right panel), which were taken at the beginning (blue box) and the end (orange triangle) of the observation. The continuum profile of the second spectrum after correction (red circle) matched with the first spectrum. Bottom panel: ratio of the spectra before (blue circle) and after (black star) correction. The orange line is the smoothed ratio profile used to correct the variation.

Appendix B SYSREM algorithm

As in Tamuz et al. (2005), the aim is to find the two sets of effective extinction coefficients () and air mass () that optimally describe the atmospheric absorption in each wavelength bin () of each frame (). By taking the observed air mass as the first input of , we search for the optimal that minimizes


where is the uncertainty of pixels calculated by taking the root sum square of the standard deviations of its frame, and the wavelength bin. Then, can be estimated by


Then, by using the estimated , the “optimized air mass” () can also be found that minimizes


by calculating


By using the “optimized air mass”, the “optimized coefficient” can be calculated and by performing this “criss-cross” iteration the stable value of , can be found. The residual can be calculated as


At this point, the first systematic effect has been removed; then, by performing a similar calculation to find , that minimize


the second and subsequent systematics can be calculated and removed. Note that does not actually represent the real extinction coefficient and air mass even for the first iteration, but a linear systematic effect that varies as a function of wavelength and time (or frame number).

Appendix C Cross section of molecules

The TiO line list from Plez (1998) was used, which include five different isotopes with nine electronics systems, and the partition function that was published by Barklem & Collet (2016). As Py4CATS only supports HITRAN- and GEISA-like databases, instead of using its extract module, we used our custom build Python script to extract the TiO line list and wrote them in HITRAN-like format, which then were used to calculate line-by-line cross-sections. The cross-section () of each line is a product of the line strength () and a normalized line profile function ()


where is the line broadening half width at half maximum (HFWHM), is the frequency, and is the line centroid position. We modified lbl2xs.py at the adjusting line parameter of the and section to enable it to calculate line strength using other partition function databases, and to include thermal (Doppler), natural, and van der Waals broadening for TiO in the line profile calculations. The line strength at temperature T () is calculated by using Equation (1) in Sharp & Burrows (2007). A Voigt function was used for the line profile, which is defined as


where are defined as dimensionless variables in terms of distance from the line centroid position, , and the ratio of Lorentz and Gaussian HWHM , :


The thermal broadening is expressed with a Gaussian profile ():


while natural and van der Waals broadening is expressed with a Lorentz profile ()


where is Boltzmann’s constant, is temperature, is the mass of the molecular absorber, and is the speed of light. is a sum of van der Waals () and natural () line broadening HWHM,


As there was no information available for the van der Waals line broadening width of both molecules, we calculated it by following Sharp & Burrows (2007) as




where is the lower rotational quantum number of each energy transition, is the FWHM of a transition at 1 atm () when , is a scale factor of the dependency of the broadening on , and as suggested in the paper we used cm, cm, and 0. The minimum criterion in Equation C8 means that the line broadening at is used for larger values. We used the formalism by Gray (1976) for the natural broadening width:



  1. Adopted from Kovcs et al. 2013
  2. Adopted from Smith et al. (2011)
  3. Adopted from Collier Cameron et al. (2010)
  4. Adopted from Collier Cameron et al. (2010)
  5. Adopted from Collier Cameron et al. (2010)
  6. Adopted from Kovcs et al. 2013
  7. Adopted from Collier Cameron et al. (2010)
  8. Adopted from Collier Cameron et al. (2010)
  9. Adopted from Collier Cameron et al. (2010)
  10. Adopted from Johnson et al. (2015)
  11. Adopted from Johnson et al. (2015)
  12. Adopted from Kovcs et al. 2013
  13. Adopted from Kovcs et al. 2013
  14. Adopted from Kovcs et al. 2013
  15. Adopted from Kovcs et al. 2013
  16. Adopted from Johnson et al. (2015)
  17. Adopted from Kovcs et al. 2013
  18. Adopted from Kovcs et al. 2013
  19. Adopted from Kovcs et al. 2013
  20. Adopted from Kovcs et al. 2013
  21. Adopted from Smith et al. (2011)
  22. Adopted from Collier Cameron et al. (2010)
  23. The Image Reduction and Analysis Facility (IRAF) is distributed by the US National Optical Astronomy Observatories, operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the National Science Foundation.
  24. http://www.subarutelescope.org/Observing/Instruments/HDS/index.html
  25. https://www.naoj.org/Observing/Instruments/HDS/hdsql/hdsql-cl-20170807.tar.gz
  26. https://bitbucket.org/lucashnegri
  27. see http://www.libradtran.org for details
  28. The nodal precession of WASP-33b caused the orbital inclination to evolve from 86.61 in 2008 to 88.70 in 2014 (Johnson et al., 2015), and based on Smith et al. (2011) analysis on its orbital eccentricity, we can safely assume a circular orbit with orbital inclination 90 to calculate RV
  29. The output from the scipy module is the two-tailed -value.
  30. http://www.hs.uni-hamburg.de/DE/Ins/Per/Czesla/PyA/PyA/index.html


  1. Barklem, P. S., & Collet, R. 2016, A&A, 588, A96
  2. Birkby, J. L., de Kok, R. J., Brogi, M., et al. 2013, MNRAS, 436, L35
  3. Birkby, J. L., de Kok, R. J., Brogi, M., Schwarz, H., & Snellen, I. A. G. 2017, AJ, 153, 138
  4. Brogi, M., de Kok, R. J., Albrecht, S., et al. 2016, ApJ, 817, 106
  5. Brogi, M., Line, M., Bean, J., Désert, J.-M., & Schwarz, H. 2017, ApJ, 839, L2
  6. Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2013, ApJ, 767, 27
  7. Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, Reviews of Modern Physics, 73, 719
  8. Burrows, A., & Sharp, C. M. 1999, ApJ, 512, 843
  9. Coelho, P. R. T. 2014, MNRAS, 440, 1027
  10. Collier Cameron, A., Pollacco, D., Street, R. A., et al. 2006, MNRAS, 373, 799
  11. Collier Cameron, A., Guenther, E., Smalley, B., et al. 2010, MNRAS, 407, 507
  12. Crossfield, I. J. M., Barman, T., & Hansen, B. M. S. 2011, ApJ, 736, 132
  13. de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82
  14. Diamond-Lowe, H., Stevenson, K. B., Bean, J. L., Line, M. R., & Fortney, J. J. 2014, ApJ, 796, 66
  15. Esteves, L. J., de Mooij, E. J. W., Jayawardhana, R., Watson, C., & de Kok, R. 2017, AJ, 153, 268
  16. Evans, T. M., Sing, D. K., Wakeford, H. R., et al. 2016, ApJ, 822, L4
  17. Evans, T. M., Sing, D. K., Kataria, T., et al. 2017, Nature, 548, 58
  18. Fortney, J. J., Lodders, K., Marley, M. S., & Freedman, R. S. 2008, ApJ, 678, 1419
  19. Gontcharov, G. A. 2006, Astronomical and Astrophysical Transactions, 25, 145
  20. Hansen, C. J., Schwartz, J. C., & Cowan, N. B. 2014, MNRAS, 444, 3632
  21. Haynes, K., Mandell, A. M., Madhusudhan, N., Deming, D., & Knutson, H. 2015, ApJ, 806, 146
  22. Hoeijmakers, H. J., de Kok, R. J., Snellen, I. A. G., et al. 2015, A&A, 575, A20
  23. Hubeny, I., Burrows, A., & Sudarsky, D. 2003, ApJ, 594, 1011
  24. Johnson, M. C., Cochran, W. D., Collier Cameron, A., & Bayliss, D. 2015, ApJ, 810, L23
  25. Kaeufl, H.-U., Ballester, P., Biereichel, P., et al. 2004, in Proc. SPIE, Vol. 5492, Ground-based Instrumentation for Astronomy, ed. A. F. M. Moorwood & M. Iye, 1218–1227
  26. Kawahara, H. 2012, ApJ, 760, L13
  27. Kawahara, H., & Hirano, T. 2014, ArXiv e-prints, arXiv:1409.5740
  28. Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, ApJ, 673, 526
  29. Knutson, H. A., Howard, A. W., & Isaacson, H. 2010, ApJ, 720, 1569
  30. Kotani, T., Tamura, M., Suto, H., et al. 2014, in Proc. SPIE, Vol. 9147, Ground-based and Airborne Instrumentation for Astronomy V, 914714
  31. Kovcs, G. and Kovcs, T. and Hartman, J. D. 2013, A&A, 553, A44
  32. Mazeh, T., Tamuz, O., & Zucker, S. 2007, in Astronomical Society of the Pacific Conference Series, Vol. 366, Transiting Extrapolar Planets Workshop, ed. C. Afonso, D. Weldrake, & T. Henning, 119
  33. Narita, N., Suto, Y., Winn, J. N., et al. 2005, PASJ, 57, 471
  34. Noguchi, K., Aoki, W., Kawanomoto, S., et al. 2002, PASJ, 54, 855
  35. Ofir, A., Alonso, R., Bonomo, A. S., et al. 2010, MNRAS, 404, L99
  36. Plez, B. 1998, A&A, 337, 495
  37. Schwarz, H., Brogi, M., de Kok, R., Birkby, J., & Snellen, I. 2015, A&A, 576, A111
  38. Sedaghati, E., Boffin, H. M. J., MacDonald, R. J., et al. 2017, Nature, 549, 238
  39. Sharp, C. M., & Burrows, A. 2007, ApJS, 168, 140
  40. Smith, A. M. S., Anderson, D. R., Skillen, I., Collier Cameron, A., & Smalley, B. 2011, MNRAS, 416, 2096
  41. Snellen, I., de Kok, R., Birkby, J. L., et al. 2015, A&A, 576, A59
  42. Snellen, I. A. G., Albrecht, S., de Mooij, E. J. W., & Le Poole, R. S. 2008, A&A, 487, 357
  43. Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Nature, 509, 63
  44. Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049
  45. Spiegel, D. S., Silverio, K., & Burrows, A. 2009, ApJ, 699, 1487
  46. Tajitsu, A., Aoki, W., Kawanomoto, S., & Narita, N. 2010, Publications of the National Astronomical Observatory of Japan, 13, 1
  47. Tajitsu, A., Aoki, W., & Yamamuro, T. 2012, PASJ, 64, 77
  48. Tamuz, O., Mazeh, T., & Zucker, S. 2005, MNRAS, 356, 1466
  49. von Essen, C., Mallonn, M., Albrecht, S., et al. 2015, A&A, 584, A75
  50. Winn, J. N., Suto, Y., Turner, E. L., et al. 2004, PASJ, 56, 655
  51. Zellem, R. T., Lewis, N. K., Knutson, H. A., et al. 2014, ApJ, 790, 53
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description