The GTC exoplanet transit spectroscopy survey X. Stellar spots versus Rayleigh scattering: the case of HAT-P-11b
F. Murgas G. Chen E. Pallé L. Nortmann G. Nowak




Stellar spots versus Rayleigh scattering: the case of HAT-P-11b



Instituto de Astrofísica de Canarias (IAC), E-38205 La Laguna, Tenerife, Spain

              fmurgas@iac.es 
              Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38206 La Laguna, Tenerife, Spain
              Key Laboratory of Planetary Sciences, Purple Mountain Observatory, Chinese Academy of Sciences, Nanjing 210008, China






 
  
  
  Rayleigh scattering in a hydrogen-dominated exoplanet atmosphere can be detected from ground or space based telescopes, however, stellar activity in the form of spots can mimic Rayleigh scattering in the observed transmission spectrum. Quantifying this phenomena is key to our correct interpretation of exoplanet atmospheric properties.
  
   We use the 10 m telescope Gran Telescopio Canarias (GTC) to carry out a ground-based transmission spectra survey of extrasolar planets to characterize their atmospheres. In this paper we investigate the exoplanet HAT-P-11b, a Neptune-size planet orbiting an active K-type star.
  
   We obtained long-slit optical spectroscopy of two transits of HAT-P-11b with the Optical System for Imaging and low-Intermediate-Resolution Integrated Spectroscopy (OSIRIS) on August 30 2016 and September 25 2017. We integrated the spectrum of HAT-P-11 and one reference star in several spectroscopic channels across the λ∼ 400-785 nm region, creating numerous light curves of the transits. We fit analytic transit curves to the data taking into account the systematic effects and red noise present in the time series in an effort to measure the change of the planet-to-star radius ratio (R_p/R_s) across wavelength.
   
   By fitting both transits together, we find a slope in the transmission spectrum showing an increase of the planetary radius towards blue wavelengths. A closer inspection to the transmission spectrum of the individual data sets reveals that the first transit presents this slope while the transmission spectrum of the second data set is flat. Additionally we detect hints of Na absorption in the first night, but not in the second. We conclude that the transmission spectrum slope and Na absorption excess found in the first transit observation are caused by unocculted stellar spots. Modeling the contribution of unocculted spots to reproduce the results of the first night we find a spot filling factor of δ=0.62^+0.20_-0.17 and a spot-to-photosphere temperature difference of Δ T = 429^+184_-299 K.
  
   

   

   The GTC exoplanet transit spectroscopy survey X.
    F. Murgas12
          G. Chen 123
          E. Pallé12
          L. Nortmann12
          G. Nowak12

    Received August 10, 2018; accepted December 22, 2018
========================================================================================================






§ INTRODUCTION


A popular technique to study the atmospheres of transiting exoplanets is transmission spectroscopy. This method consists in measuring the change in the planetary radius across wavelengths during a transit event. In the past decade this procedure has been used to detect several planetary atmospheric features such as atomic and molecular absorption lines (e.g., , , , , , , , , ) and scattering features (e.g., , , , , , , ). Thanks to their relatively large atmospheric scale heights most of these detections have been accomplished for Jupiter analogs, however there has been several efforts in recent years to study the atmospheres of exoplanets with masses around the Saturn/Neptune regime and lower (e.g., , , , , , , , , , , ). 

Despite the proven success of transmission spectroscopy, there is still a risk of attributing to the planet certain features whose origin may be on star. One of these features is Rayleigh scattering. Stellar activity in the form of spots can introduce signals that mimic the effect of an increase of the measured planetary radius toward the blue (e.g., , , ).

The exoplanet HAT-P-11b is a Neptune-size planet discovered by <cit.> orbiting a relatively bright (V = 9.47 mag) active K dwarf star ( and references therein) with a period of 4.88 days. With an initial measured radius of R_p = 0.422 ± 0.014  R_J (R_p = 4.73 ± 0.16  R_⊕) and a mass of M_p = 0.081 ± 0.009  M_J (M_p = 25.8 ± 2.9  M_⊕), HAT-P-11b had the smallest radius discovered by a ground based survey at the time of discovery. Using Kepler data, <cit.> improved the planetary radius measurement to R_p = 4.31 ± 0.06  R_⊕, a smaller radius than the reported in the discovery. <cit.> took radial velocity (RV) measurements during a transit and found that the orbit of HAT-P-11b is highly misaligned with respect to the rotational axis of its host star with an angle of 103^+26_-10 deg. Using Subaru RV data <cit.> found an angle 103^+22_-18 deg, in agreement with <cit.>. <cit.> studied the radio emission of the system at 150 MHz during an eclipse and detected radio emission which they attributed to the planet, although it has not been confirmed by follow up observations to date. <cit.> detected the secondary eclipse of this planet using data from Kepler, obtaining an eccentricity of e=0.26459^+0.00069_-0.00048 and a geometrical albedo similar to Neptune of 0.39± 0.07. Using a decade of RV measurements, <cit.> found a second planet orbiting HAT-P-11; with a mass similar to Jupiter (M_psin i=1.6± 0.1; M_J), a period of 9.3 years, and an eccentricity of e=0.6±0.03; HAT-P-11c could explain the orbital tilt of HAT-P-11b.

<cit.> obtained a transmission spectrum of HAT-P-11b using Hubble Space Telescope (HST) instrument Wide Field Camera 3 (WFC3) in the wavelength range of 1.1-1.7 μm coupled with Spitzer observations at 3.6 and 4.5 μm. The data showed evidence of water absorption at 1.4 μm. <cit.> re-analyzed HST/WFC3 transit observations of several gaseous exoplanets, including the HAT-P-11b data set taken by <cit.> confirming the absorption feature of H_2O.

Here we present the optical transmission spectrum of HAT-P-11b obtained with the 10.4 m Gran Telescopio Canarias (GTC). This paper is organized as follows. In <ref> we describe the observations and data reduction, in <ref> we describe the light curve fitting process, in <ref> we show the transmission spectrum of HAT-P-11b and discuss our results. Finally in <ref> we present the conclusions of this work.




§ OBSERVATIONS AND DATA REDUCTION


The data were taken with the 10.4 m telescope GTC at the Observatorio Roque de los Muchachos, in Spain. Using its Optical System for Imaging and low-Resolution-Integrated Spectroscopy (OSIRIS, ) instrument in its long-slit spectroscopic mode.

Two transits of HAT-P-11b were observed with OSIRIS in the nights of August 30 2016 (hereafter N1) and September 25 2017 (hereafter N2). For both observations we were able to put HAT-P-11 and one reference star inside the slit (see Table <ref> for the star properties and Figure <ref> for the acquisition image and an example of a raw science frame for each night). The data were collected using the R1000B grism (λ∼ 350 - 787 nm of coverage) and a custom built slit of 40 arcsec of width. The instrumental configuration for both nights is described in Table <ref>. For the first night the exposure time was set to 4.25 seconds and the total time of observation was 4.27 hours (from 20:43 UT to 01:00 UT) which translated into 693 science images. At 21:57 UT we lost 9 minutes of data due to problems with the telescope software, after the problem was solved the stars were reacquired at the same initial position inside the slit. For the night of September 25 2017 the exposure time was set to 3.50 seconds and the total time of observation was 5.07 hours (from 20:14 UT to 01:16 UT) which translated into 803 science images. Figure <ref> shows the count level for the pixel with the maximum flux inside the extraction aperture for each star and night[saturation levels taken from <http://www.gtc.iac.es/instruments/osiris/media/OSIRIS-USER-MANUAL_v3_1.pdf>].





The data were reduced using the same approach described in <cit.>. The science images were calibrated using a standard procedure (bias and flat correction), the wavelength correction was done using the 2D arc lamp images to create a pixel-to-wavelength transformation map. The spectra were extracted using an optimal extraction algorithm (). A fixed extraction aperture size of 42 and 37 pixels were adopted for N1 and N2 respectively. The aperture size selected for the extraction was applied to both target and reference star and is the one that delivered the lowest scatter in the points out of transit of the white light curve, i.e. the curve produced by integrating the flux of both target and reference star in the λ 399-785 nm region (excluding the region near the O_2 telluric line, λ 757-768 nm). The mid exposure time was computed using the mid point of the UT values of the opening and closening of the shutter, these mid exposure times were transformed to Barycentric Dynamical Time with the help of the <cit.> code.





Figure <ref> shows an example of the extracted spectrum of HAT-P-11 and the reference star. Note that the reference is 2.78 mag fainter than the target in V band, so the flux of the reference star is arbitrarily multiplied in the figure by a factor of 14.



   


§ DATA ANALYSIS


The data for both nights were fitted following the procedure of <cit.>. The analytic transit curve was created using the models of <cit.> through the  implementation of <cit.>. For the limb darkening coefficients (LDC) we adopted a quadratic limb darkening law; the coefficient values u_1 and u_2 were computed interpolating ATLAS (<cit.> synthetic stellar spectrum with the code of <cit.>. The coefficients were estimated using the stellar parameters of HAT-P-11 presented in <cit.> (T_eff=4750 K, [Fe/H] =0.3 dex, and log g_⋆  (cgs)= 4.5). To model the white light curve transit we set as free parameters: the planet-to-star radius ratio R_p/R_s, the quadratic limb darkening coefficients u_1 and u_2, the central time of transit T_mid (for N1 and N2), the orbital semi-major axis over stellar radius a/R_s, and the orbital inclination i. For the spectroscopic light curves the orbital parameters i and a/R_s were fixed to the values found for the white light curve. The orbital eccentricity e, period, and argument of the periastron ω were fixed to the values presented in <cit.> for all the curves analyzed here.

We created a common-mode noise model by dividing the white light curve by the best fitting transit model, and then each spectroscopic light curve was divided by this common-mode noise. To model the systematic effects present in the data, we used the Gaussian Processes (GP) python implementation provided by  (). The use of GP to model the red noise has become a standard procedure in light curve analysis (e.g., ). For all the red noise sources, we selected a squared exponential kernel described by


    k(x_i,x_j) = A^2 exp[ -∑^N_α=1( x_α,i-x_α,j/L_α)^2 ]

where x_i,j are the input vectors associated with the origin of the red noise, in this case a time component and seeing variation measured by the full width at half maximum (FWHM) of the spectral profile in the spatial direction of the target.

Near the position of the target there is a faint star whose flux is inside the spectral extraction aperture for both nights. GTC was able to observe a transit of HAT-P-11b in August 12 2017, unfortunately the target was saturated and in non-linearity range in most of the in-transit images so we decided to not use this data set in the transmission spectroscopy analysis. The seeing on that night was stable enough to allow us to separate the flux of HAT-P-11 and the contaminant star (see Figure <ref>). The flux of the contaminating star was extracted from the non-saturated pre-transit frames of the August 12 2017 data set (see Appendix). We included the flux contamination level by this contaminant in the MCMC procedure to compute an accurate planet-to-star radius ratio (R_p/R_s) across the probed wavelength range. 



A likelihood function was evaluated iteratively using the python Markov chain Monte Carlo (MCMC) suite  (). For the white light curve we used normal priors for the following parameters: the LDC with values centered around the predicted coefficients and with a normal distribution width of 0.1 plus the LDC constrains by <cit.>, the orbital inclination and semi-major axis over stellar radius (mean centered around  values), and the flux contamination fraction of the contaminant star. For the spectroscopic light curves we used normal priors for the LDC and the flux contamination fraction of the contaminant star. For the rest of the parameters we used uniform priors for both white light and spectroscopic curves.

The MCMC for the white light curve consisted of 90 independent chains (or walkers) and run for 3000 iterations. For the spectroscopic curves we used 32 independent chains and 3000 iterations. After discarding the first 500 iterations as burn-in period, we computed the median and 1σ values of the posterior distribution for each parameter.




§ RESULTS AND DISCUSSION





 §.§ White light curve

The white light curves for N1 and N2 are shown in Figures <ref> and <ref>. Since the reference star is ∼ 2.7 magnitudes fainter than HAT-P-11 in V band, the reference star was more affected by seeing variations during the observations and it is the major source of point-to-point variation in the light curve. The seeing variations affect the number of pixels where the point spread function (PSF) is contained and this effect is correlated with the FWHM of the stellar flux. Our light curve modeling is able to reproduce this noise source and we can recover the transit for each night. Once we subtract the transit model from the light curves, we reach a root mean square (RMS) of the residuals of 519 ppm and 600 ppm for N1 and N2 respectively. The transit parameters and 1σ uncertainties fitted using a joint analysis of N1 and N2 are listed in Table <ref>.










 §.§ Transmission spectrum

The GTC/OSIRIS transmission spectrum for the joint fit of N1 and N2 is shown in Figure <ref>. Our joint fit results are consistent with a pure Rayleigh scattering feature when compared to a flat model (straight line) and  () models with clear (i.e. no clouds) and cloudy atmospheres. The cross-section of the particles that produce the Rayleigh scattering follow a power law σ = σ_0(λ / λ_0)^α, where λ_0 is the cross-section size at the reference wavelength λ_0 and α=-4 (). The change of measured planet-to-star radius ratio across wavelength is related to the planetary atmospheric scale height and α by 


    d R_p/d lnλ = α H = ακ_B T_p/μ_m g_p

where H is the atmospheric scale height, κ_B the Boltzmann constant, T_p the temperature of the atmosphere of the planet, μ_m the mean molecular weight of the planetary atmosphere, and g_p is the planet's surface gravity.

Assuming a hydrogen-dominated atmosphere (μ_m=2.37), a planetary temperature of T_p = 878 ± 15 K (), and adopting g_p=14.08 ± 1.17  m/s^2 (computed using the planetary mass of  and radius of ), we find α_joint = -9.86 ± 6.19 for the measured d R_p / d lnλ using the joint fit of N1 and N2. It is worth mentioning that the relatively large uncertainty in α_joint comes from the stellar radius of HAT-P-11 (eq. <ref> needs the planetary radius, not the planet-to-star radius ratio). Considering the uncertainties, the α value for the joint fit is consistent with Rayleigh scattering from a hydrogen dominated atmosphere, although the measured value departs from the expected α=-4. This led us to examine in closer detail each observing night.

We re-analyzed the transits for N1 and N2, this time leaving the planet-to-star radius ratio as free parameter for each night independently. Figure <ref> presents a comparison of the transmission spectrum of HAT-P-11b for each night fitted independently along with the results of the joint fit. The data from N1 presents a greater R_p/R_s versus lnλ slope; if we use the same μ_m, T_p, and g_p values as before we get α_N1 = -21.85 ± 9.25. In contrast the transmission spectrum for N2 is relatively flat and featureless with a slope slightly increasing towards redder wavelengths. One possible explanation for the discrepancy between data sets is stellar activity.










 §.§ The impact of stellar spots

The host star of HAT-P-11b is known to be active, in the discovery paper <cit.> measured a mean Ca ii H & K emission activity index of ⟨ S ⟩= 0.61. <cit.> studied the evolution of the S-index of HAT-P-11 from 2008 to 2017 using Keck/HIRES data, they found evidence of an activity cycle close to 10 years of duration and measured a mean S-index of ⟨ S ⟩= 0.58 ± 0.04. The mean S-index measured by <cit.> is equivalent to an activity index of log R'_HK = -4.35 (see  for the definition of the R'_HK index), a value representative for stars exhibiting significant stellar activity.

One way in which stellar activity affects the star is the appearance of spots in its surface. There are several transits of HAT-P-11b observed by Kepler where the planet crossed regions with spots (for example see  Fig. 23). <cit.> used Kepler data to constrain the orbital obliquity and identified two regions with spots which were interpreted as active regions of the star. Also using Kepler data <cit.> found evidence of a long-lived spot region in HAT-P-11 that was eclipsed by the planet 6 times and discovered that the orbital period of HAT-P-11b and the rotational period of the host star are in a period ratio of 6:1. <cit.> studied the spot area coverage and spot distribution across the stellar surface, deducing that 3^+6_-1 % of the surface is covered by spots and that they are distributed around 16^∘± 1^∘ latitude.

Occulted and unocculted stellar spots can affect the flux measurement during a transit, leading to a wrong estimate of the transit depth (). According to <cit.> the median amplitude of spot crossing events in Kepler band is close to 500 ppm, by visual inspection we did not detect any spot crossing event in both of our data sets which both have rms of the residuals of the order of 500 ppm in the white light curve. The major problem preventing us from detecting occulted spots is the faintness of the reference star, which introduced strong systematics correlated with seeing. Even if there were spot crossing events in our data sets, their effect on the measured transmission spectrum would be negligible due to the modeling of common-mode noise together with Gaussian Processes (CMN removes any abrupt change common in wavelength, GP treat the spot as time-correlated noise). However, there is still a chance that unocculted spots can be responsible of producing a signal similar to Rayleigh scattering in the transmitted spectrum of the planet for N1. 

Since spots are regions with different temperature than the average surface temperature, they can introduce a slightly different stellar flux level in some wavelengths creating an effect similar to Rayleigh scattering. To model the effect of unocculted spots during the HAT-P-11 b transit of N1 we followed the approach of prior works (e.g., , , , ) in which the observed planet-to-star radius ratio (without including the effect of limb darkening) can be modeled by 


    ( R̅_̅p̅/R̅_̅s̅)^2 = ( R_p/R_s)^2 1/1 - δ (1 - F_ν(spot)/F_ν(phot))

where R̅_̅p̅/R̅_̅s̅ is the observed planet-to-star radius ratio, R_p/R_s is the true planet-to-star radius ratio (i.e. not affected by spots), δ is the variable that takes into account the stellar surface area covered by spots or filling factor, F_ν(spot) is the flux from the stellar spots, and F_ν(phot) is the flux from the stellar photosphere.

To model the flux from the star (and spotted regions) we used PHOENIX stellar models () for HAT-P-11 (the model with stellar parameters of T_eff=4780 K, log g_⋆  (cgs)= 4.5, [Fe/H] =0.5 dex) and compared it with a grid of spectra with temperatures ranging from 2700 K to 7000 K.

We performed a MCMC fitting procedure using equation <ref> applied to the transmission spectrum of N1 and N2, after removing the last data point (bin centered at 777.5 nm) which presents a different R_p/R_s than the rest of the bins probably due to noise added by the nearby telluric O_2 line (λ=760.5 nm, see Fig. <ref>). We want to emphasize that the main noise source in this bin is the contaminant star whose wavelength is shifted due to misalignment in the slit (see Fig. <ref>). If there were no such a contaminant star, the majority of the telluric O_2 band would have been excluded from our flux integration (see Fig. <ref>), and the remaining effects should have been corrected by the reference star.

The MCMC was launched with  using 80 independent chains and two separate runs: a burn-in run with 10000 iterations and a main run with 25000 iterations. Using the posterior distributions to get the median values and 1σ uncertainty range, we get for N1: δ_N1 = 0.62^+0.20_-0.17 and T_spot  N1=4351^+184_-299 K. <cit.> reported a filling factor of 3^+6_-1 % for HAT-P-11, a value that the authors suggest as a lower limit due to the assumption made that the stellar spot coverage probed by the transit is representative of the spot coverage of the whole star and the sensitivity of their measurements. In the same article, they claim that HAT-P-11's filling factor is close to the values of the active K star OU Gem presented in <cit.> with numbers ranging from δ≤ 0.04 up to 0.35. Our median value for the filling factor is higher than the values for OU GEM, however the lower limit of the filling factor presented here (δ_N1  low = 0.45) is close to the numbers presented in <cit.> and typical values seen in active mid-K dwarfs (e.g., , ) and M dwarfs (e.g., , ). The difference in temperature between the star and the spot is Δ T = T_phot-T_spot  N1= 429^+184_-299 K (considering a stellar model with a temperature of T_phot = 4780 K), which is almost half the temperature presented in <cit.> who measured a Δ T ≈ 900 K based on the largest star spot crossing event observed by their Kepler and Spitzer runs. The surface-to-spot temperature difference of the largest spot crossing event observed by <cit.> is consistent with values found in active K dwarfs, which usually are in the range of 1000 - 1500 K (e.g., , ).

The relative large value and uncertainties found in δ are probably caused by two factors: 1) the uncertainties of the measured planet-to-star radius ratio allowed several combinations of filling factors and spot temperatures to reproduce the observed transmission spectrum, and 2) the filling factor and the temperature of the stellar spot (which gives F_ν(spot)) are correlated, thus, producing a very wide posterior probability distribution for the δ parameter.

The second data set could also be affected by bright and dark unocculted spots whose effects cancel out producing the measured flat transmission spectrum (see Fig. <ref>). It is possible to extend Eq. <ref> to include the effects of bright and dark unocculted spots (see ), however due to our uncertainties and the strong correlation between the filling factor and temperature of the spot we decided to use the same MCMC procedure applied for N1 for simplicity. For the flat transmission spectrum of N2 we get δ_N2 = 0.19^+0.30_-0.14 and T_spot  N2=4724^+422_-1488 K, a lower filling factor compared to N1 and a spot temperature closer to the adopted temperature of the star.



As a final test we used equation <ref> to compute the filling factor and spot temperature necessary to go from the transmission spectrum with a slope of N1 to the featureless transmission spectrum of N2. Using the same MCMC fitting procedure as before and without any constrains in the temperature of the spot we get δ = 0.042^+0.278_-0.036 and T_spot=4802^+1073_-891 K. If we put the constraint that the spot temperature must be higher than the average stellar surface temperature (i.e., bright spots only), we get δ = 0.045^+0.317_-0.041 and T_spot=5078^+2016_-263 K. Both results are consistent with each other considering the uncertainties, although the second test preferred higher spot temperatures (Δ T = 298^+2016_-263 K) compared to the unconstrained temperature test (Δ T = 22^+1073_-891 K). The low filling factor values obtained in these tests indicate that a relatively small amount of bright spots can suppress the slope found in N1 producing the flat transmission spectrum of N2.

Another way to establish if unocculted spots affected the transmission spectrum is to check the measured transit depth around some atomic lines. For example spots can create an excess in the transit depth at the Na i 589.0 and 589.6 nm doublet while the transit depth decreases in H_α (656.3 nm) (see ). We fitted (individually for N1 and N2) several light curves with a bin size of 5 nm of width around the Na i doublet, but this time using the mean FWHM of the absorption lines as input for the GP instead of the FWHM of the spectral profile used in the previous presented results. The reason for this change is that using the FWHM of the absorption lines delivers better results when modeling the red noise in curves produced with narrow wavelength bins.

Figure <ref> shows the planet-to-star radius ratio measured around the Na i doublet for N1 (top panel) and N2 (bottom panel), compared to the same  models presented in Fig. <ref> and the expected R_p/R_s using the best fitting model of Eq. <ref> for N1 and N2. The first dataset (N1) presents an excess in the transit depth for the bin centered at the Na i doublet in agreement with the atmospheric models for HAT-P-11b and the unocculted spot model. In the second data set (N2) the excess in R_p/R_s at the bin centered at the Na i doublet disappears, and the overall transmission spectrum is more consistent with a flat model, as predicted by the unocculted spot model for N2.

In the case of H_α our measurements are not precise enough to detect the decrease in R_p/R_s around this line and differentiate it from a flat model (see Fig. <ref>).



Based on the transmission spectrum slope, unocculted spot modeling, and the transit depth excess detected in Na i in N1 and not in N2, we are confident that the measured transmission spectrum of N1 was affected by stellar spots mimicking a feature similar to Rayleigh scattering.




§ CONCLUSIONS


We present here an analysis of two primary transits of the Neptune-size exoplanet HAT-P-11b taken with GTC/OSIRIS instrument. We were able to obtain long-slit spectra of HAT-P-11 and one reference star for the nights of August 30 2016 and September 25 2017. Integrating the flux of both stars in different passbands, we created several light curves to measure the change in transit depth across the wavelength range of λ 400-785 nm. The light curves for both nights were fitted jointly and individually using a Bayesian MCMC procedure that takes into account the systematic effects present in the data.

The GTC/OSIRIS transmission spectrum of HAT-P-11b obtained with a joint MCMC fit presents a slope towards blue wavelengths. The joint fit results are consistent with a pure Rayleigh scattering feature when compared to a flat model (straight line) and synthetic atmosphere models with clear and cloudy atmospheres. Assuming that this slope was caused by small particles present in the atmosphere of the planet and adopting a planetary equilibrium temperature of T_p = 878 ± 15 K and a molecular weight of μ_m=2.37, we measured α_joint = -9.86 ± 6.19. Considering uncertainties this α_joint value is in the expected Rayleigh scattering range caused by hydrogen-dominated atmospheres (α = -4), although strong assumptions were made when computing this parameter (temperature of the planet and molecular weight). A closer inspection to each individual night revealed that the transmission spectrum for HAT-P-11b taken on August 30 2016 has a significant slope when compared to the transmission spectrum on the night of September 25 2017 which is mostly featureless and flat. 

The difference in the observed transmission spectrum between the two data sets could be attributed to stellar activity, more specifically to a greater number of unocculted stellar spots during N1. When modeling the amount of stellar surface area covered by dark spots and their temperature to reproduce the results of the August 2016 data set, we find δ_N1 = 0.62^+0.20_-0.17, T_spot  N1=4351^+184_-299 K and Δ T = T_star - T_spot  N1 = 429^+184_-299 K. Although our uncertainties in the filling factor are relatively large, the lower limit is in agreement with values seen in active mid-K dwarfs and M stars. The origin of the large uncertainties in the filling factor comes from the fact that δ and T_spot are strongly correlated.

Additionally we modeled the filling factor and spot temperature needed to suppress the slope of the first data set to obtain the flat transmission spectrum of the second night for two cases: 1) setting the spot temperature completely free, and 2) forcing the spot temperature to be greater than the average stellar surface temperature (i.e., bright spots). For the first case we get δ = 0.042^+0.278_-0.036 and T_spot=4802^+1073_-891 K and for the second δ = 0.045^+0.317_-0.041 and T_spot=5078^+2016_-263 K. Based on these results a relative small amount of bright spots present in the star in the second night could produce the observed flat transmission spectrum, assuming that the spectrum of N1 is a closer representation to the true spectrum of the planet.

Since stellar spots can also affect the intensity of some spectral lines we also compared the transit depth around the Na i doublet (λ 589.0 nm and 589.6 nm) for both nights, finding that the transit in the first data set is deeper around the Na i doublet than the nearby continuum while for the second night the transmission spectrum is flat; this is consistent with a major presence of unocculted dark spots during the transit of N1 compared to N2. The same analysis was performed for H_α (λ 656.3 nm), however our transit depth measurements around this line are not precise enough to detect a change between N1 and N2.

Based on the fact that the slope in the transmission spectrum is found in N1 but not N2 and the comparison between the transit depth around Na i lines for N1 and N2, we conclude that the transmission spectrum for N1 was affected by stellar activity in the form of unocculted dark spots. The second data set could also be affected by bright and dark unocculted spots whose effects cancel out producing the measured flat transmission spectrum for N2, or simply due to a low stellar activity epoch.

In light of our results for HAT-P-11b, we advocate here for repeated observations of exoplanet transmission spectroscopy studies, to account for possible stellar variability. Simultaneous/contemporary activity monitoring is crucial to establish the true nature of Rayleigh scattering features detected in exoplanets orbiting around active stars.


  Based on observations made with the Gran Telescopio Canarias (GTC), installed in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias, in the island of La Palma. G.C. also acknowledges the support by the National Natural Science Foundation of China (Grant No. 11503088) and the Natural Science Foundation of Jiangsu Province (Grant No. BK20151051).

  Software:  (),  (),  (),  (),  (<https://github.com/hpparvi/PyDE>). This research made use of Astropy, a community-developed core Python package for Astronomy (, ). Correlation plot for N1 and N2 spot modeling done with  ().














aa




  

  

§ ADDITIONAL FIGURES


  

 §.§ Saturated images for August 12 2017

The night of August 12 2017 GTC observed a transit of HAT-P-11b, however the target was saturated in a large fraction of images and we decided to not use this data set to obtain the transmission spectrum of this planet. During this night the seeing conditions were good enough to separate the flux of HAT-P-11 and a nearby star, hence we used this data set to establish the flux contamination of including this nearby star inside the extraction aperture. Figure <ref> shows the flux level of the pixel with the maximum flux inside the extraction aperture during the night (top panel) and the saturated pixels along wavelength of the time series (bottom panel). The flux of the contaminant star was extracted using the pre-transit non-saturated images.

  
  
  

 §.§ Light curve fitting tests

  We compared the transmission spectrum of HAT-P-11b obtained by fitting the spectroscopic light curves with Gaussian Processes and without GPs for each individual nights. We tested several models with different baseline functions and parameters. For the final comparison we selected the models with the lowest Bayesian Information Criterion (BIC) value. The two non-GP baseline models selected were
  
    p_1 = c_0 + c_1S_y + c_2S_y^2 + c_3t + c_4t^2 + c_5t^3

  and  
  
    p_2 = c_0 + c_1S_y + c_2S_y^2 + c_3P_y + c_4P_y^2 + c_5t + c_6t^2 + c_7t^3

  where c_i are fitted coefficients, t is the time, S_y is the FWHM in spatial direction, and P_y is the pixel position of the center of the Point Spread Function (PSF) of the spectral profile. To account for the effect of red noise in the data we used the time-averaging method ().

  Figure <ref> presents the transmission spectrum obtained with and without GPs for the baseline models p_1 and p_2. Both methods show consistent results and the slope in the transmission spectrum does not change significantly.
  
   
  
  

 §.§ Spectroscopic light curves

    In this section of the appendix we present the spectroscopic light curves of HAT-P-11b for the nights of August 30 2016 and September 25 2017.

    
    
    
    

   

 §.§ Correlation plots

   Here we present the correlation plots of the MCMC fitting procedure of Eq. <ref> for N1 and N2.
   
    
    
    
    


