EDP Sciences
Free Access
Issue
A&A
Volume 598, February 2017
Article Number A131
Number of page(s) 12
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201629473
Published online 13 February 2017

© ESO, 2017

1. Introduction

The first transiting hot Jupiter, HD 209458b, was discovered more than a decade ago (Charbonneau et al. 2000). Since then a few thousand transiting exoplanets have been detected with various transit surveys from ground-based facilities (e.g., Hartman et al. 2004; Pollacco et al. 2006), and from space (e.g., Borucki et al. 2010; Koch et al. 2010). This has provided a fertile ground for studying the exoplanet population from a global point of view (Sing et al. 2016), and has opened a window to a deeper characterization of planetary systems.

Transiting systems allow the study of their atmospheres via reflection and transmission spectroscopy (see the review paper by Burrows 2014, and references therein). During primary transit, the exoplanet blocks part of the stellar light. If the planet has an atmosphere, an additional fraction of the stellar light can be absorbed by the materials present in its atmosphere. Since this absorption takes place at discrete wavelengths, the planet size and the transit light-curve depth will appear larger or smaller, depending on the wavelength of the observation. As a result, the variations in the transit depth as a function of wavelength can reveal the presence of different atomic and molecular species (e.g., Charbonneau et al. 2002; Tinetti et al. 2007; Vidal-Madjar et al. 2011; Crossfield et al. 2011; Désert et al. 2009; Hoeijmakers et al. 2015; Nikolov et al. 2015), clouds (Gibson et al. 2013; Kreidberg et al. 2014), and hazes (Pont et al. 2013).

Most of the spectro-photometric measurements of exoplanet atmospheres have been performed through space-based observations (e.g., Deming et al. 2013; Sing et al. 2016). Although these instruments have superior stability and sensitivity, their low spectral resolution can only constrain some atmospheric models. So far, high-resolution spectroscopy has only been performed using ground-based facilities, and despite the added complications posed by the Earth’s atmosphere, ground-based high-resolution studies have added valuable contributions to our understanding of exoplanetary atmospheres. For example, at resolutions of R ~ 105, the absorption lines of individual chemical species can be spectroscopically resolved (Birkby et al. 2013; Rodler et al. 2012; de Kok et al. 2013; Lockwood et al. 2014; Brogi et al. 2012), as can the orbital motion, diurnal rotation of the planet, and exo-atmospheric wind speed (Snellen et al. 2010, 2015; Louden & Wheatley 2015; Brogi et al. 2012). Even large-scale dynamics in the upper atmosphere have already been detected (see, e.g., Kulow et al. 2014; Ehrenreich et al. 2015).

Studies of the optical transmission spectra (mainly performed by the Hubble Space Telescope, HST) have revealed that many hot Jupiters have featureless spectra with strong Rayleigh scattering slopes toward short wavelengths (Pont et al. 2008; Sing et al. 2011; Huitson et al. 2012; Nikolov et al. 2015; Sing et al. 2016). This has been attributed to the presence of high-altitude clouds and haze layers that effectively obscure absorption features by making the upper atmosphere opaque. In these cases, only species that are present in the uppermost layers of the atmosphere can be observed. Indeed, models indicate that only a handful of such species are expected to be present in the optical region (see, e.g., Fortney et al. 2010; Seager & Sasselov 2000; Brown 2001), in consequence making neutral sodium one of the most intensively studied alkali metals to date.

Exoplanetary sodium absorption was first detected by Charbonneau et al. (2002) in the atmosphere of HD 209458b by analyzing spectro-photometric data obtained by the HST during transit. Then, Redfield et al. (2008) claimed the first ground-based detection of sodium in HD 189733b and in the same year, Snellen et al. (2008) reported the ground-based detection of sodium in HD 209458b. Altogether, sodium has been detected in a handful of hot Jupiter atmospheres: HD 209458b (Charbonneau et al. 2002; Snellen et al. 2008), HD 189733b (Redfield et al. 2008; Jensen et al. 2011; Huitson et al. 2012), WASP-17b (Wood et al. 2011; Zhou & Bayliss 2012), and XO-2b (Sing et al. 2012) are classical examples. Recently, Wyttenbach et al. (2015) reported the ground-based detection of sodium in HD 189733b using observations collected during three transits with the HARPS spectrograph. Cauley et al. (2016) also analyzed the pre-transit and in-transit phases at different atomic absorption lines, consisting of sodium.

In this paper, similar to Redfield et al. (2008) and Wyttenbach et al. (2015), we present the results of our efforts in the detection of atmospheric sodium in the transmission spectrum of the hot Jupiter HD 189733b. We analyze a single transit observed with the high-resolution Ultraviolet and Visual Echelle Spectrograph (UVES) at ESO’s Very Large Telescope and we apply a detailed modeling of the systematics. We also estimate the equivalent width of the exoplanetary sodium line using the deformations introduced in the excess light curves by the changing radial velocity of the planet, and with it we estimate the abundance of sodium in the atmosphere. In Sect. 2 we describe our observations and the data reduction steps. We continue in Sect. 3 with a detailed description of the data analysis, the transmission spectroscopy method, and the modeling procedure. We show our results and the discussions in Sect. 4, and conclude in Sect. 5.

Table 1

Adopted values for the orbital and physical parameters of HD 189733 during the fitting procedures in this work.

2. Observations and data reduction

2.1. Our target

The hot Jupiter HD 189733b orbits a bright (V ~ 7.7) and active K-type star every ~2.2 days. HD 189733b has an atmosphere with a scale height of about 200 km (Désert et al. 2011). Thanks to its large scale height and bright host star, this exoplanet has turned into one of the favorite targets for atmospheric characterization. Especially in the case of ground-based studies, HD 189733b is an excellent “testbed” to explore new exoplanet atmospheric observation and characterization techniques. The orbital and physical parameters of this system adopted in this work are summarized in Table 1. To be consistent with current data, we adopt the planet-to-star radius ratio obtained by the Sing et al. (2011) HST observations around the sodium wavelength. However, we note that within the precision of our data other reasonable values of RP/RS would have produced similar results to the ones presented in this work.

2.2. Observing log and instrumental setup

thumbnail Fig. 1

Sky conditions during our observations. The airmass is indicated with a green dashed line and the red cross points indicate the seeing, in arcsec. The measured strength of several strong telluric water lines are shown with solid blue lines. These values are artificially shifted to fit inside the plot. Gray areas show the time between transit ingress and egress. Dark gray is the transit at full depth.

Open with DEXTER

During the night of July 1, 2012, a transit of HD 189733b was observed using UVES (Dekker et al. 2000) mounted on Kueyen, UT2 (Program ID: 089.D-0701 A). During the observations 244 spectra were acquired: 67 exposures before ingress, 88 exposures during transit, and 89 exposures after egress. The first 29 spectra have an exposure time of 30 s, while all of the rest were exposed for 45 s. The data were taken using the dichroic beam splitter with central wavelengths around 760 nm in the red arm and 437 nm in the blue arm. Our target lines are located in the red arm, where the slit width was 0.7′′ and the spectral resolution is about 60 000. We estimate a signal-to-noise ratio around the sodium lines of approximately 100. The sky conditions (airmass, seeing, and telluric water) are shown in Fig. 1. As the figure shows, the depths of water lines follow a trend similar to the seeing. However, the water lines do not clearly correlate with airmass. We calculate the strength of the water lines by measuring and averaging the equivalent widths of the six strongest water lines and the values for airmass and seeing are extracted from the image headers.

2.3. Data reduction

The initial data reduction and extraction of 1D spectra were performed using the UVES pipeline, version 5.2.0 (Czesla et al. 2015). The subsequent analysis related to this work was conducted using Interactive Data Language (IDL) and Python 2.7. The sequence of the reduction processes are summarized below.

2.3.1. Spectrum normalization

The first step was to normalize the continuum of the spectra in the region around the sodium lines. During the initial reduction we removed the instrumental blaze function from each individual spectral order. However, the spectra were still not properly normalized since we observed variations in the continuum. To correctly normalize the spectra we first obtained the average spectrum of all the available reduced spectra. Then we divided this average spectrum into 60 wavelength regions within a range of 400 Å around the sodium lines. By linearly interpolating the maximum values within each bin we obtained the shape of the continuum in the averaged spectrum. Finally, each individual spectrum was divided by this interpolation.

2.3.2. Spectrum alignment

Owing to the Earth’s rotation and instrumental instabilities, the wavelength solution of our spectra slightly drifts (up to about ±1 km s-1). To correct for these shifts we identified 80 stellar absorption lines in a range of 400 Å around the sodium lines. We then fitted a Gaussian profile to each one of the line cores to obtain the position of each line center. For each line we then calculated the time-averaged line position and thus obtained the deviation of each line in each order of each exposure. The strongest component is a wavelength-independent offset caused by the entire spectrum being shifted on the detector. There may also be wavelength-dependent deviations caused by a magnification, i.e., a scaling of the image of the spectra on the detector, so we performed a linear fit to the values of the misalignment of all lines in the exposure to obtain a 2D map of the misalignment (shift map) as a function of wavelength. By re-interpolating all exposures to this map, all absorption spectra were aligned, and thus placed in the same reference frame. The amplitude of the fitted linear shift is plotted in Fig. 2, which shows that the spectra drift in time according to a saw-tooth pattern, also identified and described by Czesla et al. (2015). The value of the shift in pixels rises to about ±0.5 px. To study the wavelength stability of this data set, Czesla et al. (2015) also obtained the radial velocity variations of the telluric water lines in time. Since they show the same saw-tooth behavior, we consider this to be an instrumental effect that is likely produced by the movement of the stellar seeing disk inside the slit. After removing the telluric water lines (see Sect. 2.3.3), we aligned all spectra to a common rest-frame using the shift map.

thumbnail Fig. 2

Misalignment of each spectrum with respect to a reference spectrum taken in the middle of the observing run.

Open with DEXTER

2.3.3. Telluric correction

The Earth’s telluric water and oxygen lines are removed using the telluric absorption model described in Husser & Ulbrich (2014) and Husser et al. (2016). This model determines the parameters for the stellar and Earth atmospheres simultaneously, and fits the widths and depths of the lines. Each spectrum is then divided by its corresponding telluric model to correct for telluric effects. It should be noted here that visual inspection of the telluric sodium lines in the spectra at the expected location revealed no telluric sodium feature. In addition, the object distance is relatively close (~19 pc), which agrees with the observed lack of interstellar sodium.

3. Data analysis and modeling

3.1. Excess light curve

Since the absorption cross-section of species in the atmosphere is wavelength-dependent, the planetary radius as measured by the depth of the transit light curve is wavelength-dependent as well. We therefore can infer the presence of an absorbing species by obtaining the transit light curve at that wavelength at which a given species is expected to absorb (integration band), and compare it with the transit light curve at a wavelength where the exoplanet atmosphere is transparent (reference band). This comparison is made by computing the ratio of the two light curves, referred to as the “excess light curve”, where the additional absorption at the wavelength of interest as the transit progresses can be seen (Snellen et al. 2008). In the case of sodium, we are interested in finding the excess absorption around the optical sodium D1 and D2 lines. We specifically introduce a central band (1)

where Fcenter describes the integrated flux inside the target line that is expected to change in time (t) in the case of exoplanetary atmospheric absorption, and λ0 is the central wavelength of the target absorption line. We then define left and right reference bands (2)and (3)where M, N, O, and P are suitably chosen coefficients that determine the location and width of the reference bands in the continuum close to the sodium lines, and Δλ is half the integration band width. From Fleft, Fright, and Fcenter we construct the relative line flux, Fline, defined as (4)The parameter Fline is the flux ratio of the integration and reference bands, which determines the excess light curve. In this fashion, all flux variations that affect the central and the left and right reference bands in the same way cancel out, while variations that affect only the central band such as absorption by the exoplanetary atmosphere remain.

thumbnail Fig. 3

Integration bands and derived excess light curve. Left: gray shaded areas show the integration bands with passband of 1 Å, 1.5 Å, and 3 Å, centered at the core of each sodium line. The red and blue intervals show the passbands of the reference bands for D2 and D1, respectively. Right: the raw excess light curve shown with the red circles for the 1.5 Å integration band inside the sodium D2 line. Dashed vertical lines indicate times of ingress and egress.

Open with DEXTER

In our analysis we first locate the centroids of each sodium D-line in every reduced spectrum by fitting a Voigt profile. Then we chose flux integration passbands of 1 Å, 1.5 Å, and 3 Å centered on each line to characterize the extra Na I absorption from the planetary atmosphere. The smallest passband width is set by the signal-to-noise ratio in the deep cores of the Na I lines; the aperture needs to be wide enough to include enough photons to be sensitive to small variations in the depth of the Na I lines. The upper limit in the width of the passband is given by the value at which the planetary Na I signal blends with the noise. Specifically, our 3 Å integration aperture upper limit is determined empirically by comparing the standard deviations in the light curves for each passband to that of the reference continuum. The standard deviation of the reference continuum light curve, where one does not expect to have any planetary atmospheric signal, is 0.83 parts per thousand (ppt), which we attribute to random noise. In the case of the Na I line passbands, the standard deviations of the light curves are 1.9 ppt, 1.4 ppt, and 0.85 ppt for the 1 Å, 1.5 Å, and 3 Å passbands, respectively. These values clearly indicate the presence of an additional signal in the 1 Å and 1.5 Å passbands, while the noise level of the 3 Å passband is similar to that of the reference continuum, placing the signal and the noise at the same level. The three integration bands are illustrated in Fig. 3, left, together with the sections of the spectrum we used as reference continuum. For a comparative analysis of our results with those of other authors, our integration passbands are close to the integration passbands of Snellen et al. (2008), Albrecht (2008), and Wyttenbach et al. (2015). Our reference bands are fixed to the same values for all the integration bands, and are placed closer to the upper wings of the sodium lines to minimize the difference between the limb-darkening (LD) coefficients in the integration and reference bands (see Sect. 3.2.2). We also place them in a location free of strong absorption features. Taking all of this into account, our reference passbands have a width of 1 Å and are placed on each side of each sodium line, as indicated in Fig. 3. The reference bands between the lines are taken to be the same for both sodium D-lines, otherwise for the D2 line a fraction of the reference band would fall inside the deep iron line in between the sodium lines.

Finally, we use Eq. (4) to obtain the relative strength of the sodium absorption features in each exposure. The excess light curves we obtain by using this approach on the sodium D2 line for the integration bands of 1.5 Å are shown in Fig. 3, right, and we call it the “raw excess light curve”. As the figure shows, the raw excess light curve does not look transit-like and is probably highly affected by external physical, environmental, or instrumental effects. In the next sections we investigate the source of these effects and use them to produce a detailed model fit to the data.

3.2. Model components and parameters

Before drawing any conclusions on the HD 189733b atmosphere, the external effects present in the raw excess light curves have to be taken into consideration. In total, we identify three main effects: the occurrence of a stellar flare (Czesla et al. 2015 and Klocová et al., in prep.) starting close to mid-transit time (A); a wavelength-dependent limb-darkening (B); and the observed changes in the line profile produced when the planetary sodium line moves inside the stellar sodium line, induced by the planetary orbital motion (C). In the next sections we explain each effect separately and specify their fitting parameters. The magnitudes of the individual errors for the excess light curves are determined in two ways, by computing the standard deviation of data points between the 30th and 100th exposures, and by computing the standard deviation of the residual light curves produced by subtracting a least-squares fit of the four model components arranged as in Eq. (9) to the data. For a more conservative approach, the values of the errors are assigned by choosing the largest value between the two approaches.

3.2.1. Stellar flare (component A)

HD 189733 has been characterized as a highly active K-type star (e.g., Poppenhaeger et al. 2013). As a result, our observations are prone to being contaminated by stellar flares and/or spots. Indeed, we detect a rising pattern of flux that starts right after the mid-transit point, identified as a stellar eruption. Usually, a clear evidence of a stellar flare can be derived from the measurements of the equivalent widths of emission lines in the cores of Ca II H and K and Hα lines originating in the chromospheric layers (e.g., Klocová et al., in prep.). By analyzing these spectral lines, we confirmed the presence of a flare, which – as predicted – effectively took place close to mid-transit time (see top panel of Fig. 11 in Czesla et al. 2015). Although this feature is identified using the Ca II H and K lines, stellar flares might affect other spectral lines. For instance, neutral sodium in the lower stellar atmosphere can be contaminated by this kind of eruption (Cessateur et al. 2010). Therefore, in our specific case of study it is important to consider this time-dependent change in the core of the line as part of our model budget. To mitigate the flare, we made use of the variations in the equivalent width of the Ca II K line (Klocová et al., in prep.), and we call this model component the flare profile, flrCa(t). To quantify to what extent the flare correlates with our raw sodium excess light curve, we used a Pearson’s correlation analysis. The six computed values of Pearson’s coefficient, computed from the two sodium lines times three integration bands, range between 0.25 and 0.45, and their corresponding p-values are in five of the six cases smaller than 10-3. For the 3 Å integration band and the D2 line, the p-value is 0.02. In these cases, the computed p-values indicate the probability of an uncorrelated system producing data sets that have a Pearson correlation at least as extreme as that computed from these data sets. In other words, they indicate that the null hypothesis, which states that the raw excess light curve and the flare are uncorrelated, can be rejected with high significance. Thus, there is enough statistical evidence that the flare has to be taken into consideration. To mitigate the stellar activity present in our excess light curves, we fitted a scaled version of the Ca II K flare profile to the raw excess light curves with the expression (5)where FLscale is the fitting parameter that scales the flare, and flrNa is the model component that accounts for the stellar flare. Our best-fit model for the Na D2 line at 1.5 Å integration width, showing its isolated flare component, can be seen in Fig. 4A as a solid black line, plotted along with the raw excess light curve in red points.

3.2.2. Wavelength dependent limb-darkening effect (component B)

thumbnail Fig. 4

Individual model components plotted over the excess light curve for Na D2 at 1.5 Å. Panel A): best-fit flare model component. In this panel the green dashed lines indicate the beginning and the end of the transit; panel B): best-fit differential limb-darkening model component; panel C): best-fit changing RV model component, which includes the exoplanetary atmospheric excess absorption in addition to a bump at the center. The final best-fit model, which is the combination of all the model components, is shown in Fig. 6.

Open with DEXTER

Table 2

Limb-darkening coefficients (u1, u2) for the sodium D2 and D1 integration bands (core) and fixed reference (ref.) bands with the specified wavelength ranges (λ-range).

Stellar limb darkening is wavelength-dependent. Thus, the limb-darkening coefficients inside the broad and deep lines are different from the neighboring continuum. Since the construction of the excess light curve relates the core of the Na D-lines to their near continuum, the difference in limb-darkening values will have an impact on the shape of the excess light curve (see, e.g., Charbonneau et al. 2002; Sing et al. 2008; Czesla et al. 2015). We can formulate this effect using the relation (6)where LD is the wavelength dependent limb-darkening model, LCcore; LCleft and LCright are the transit light-curve models from Mandel & Agol (2002) using the orbital parameters tabulated in Table 1, which have their own limb-darkening coefficients at the core of the sodium line, and at the left and right reference passbands, respectively (see the limb-darkening model in Fig. 4B and Fig. 2 of Czesla et al. 2015). To take this effect into consideration we first computed the limb-darkening coefficients in all relevant (integration and reference) passbands using the PHOENIX angle-resolved synthetic spectra (Hauschildt & Baron 1999; Husser et al. 2013) implementing the method detailed in von Essen et al. (2013). PHOENIX requires effective temperature, surface gravity and metallicity of the star as input, for which we adopted the values Teff = 4900 K, [Fe/H] = 0, and log g = 4.5, closely matching the stellar parameters of HD 189733 listed in this work in Table 1. The derived angle-dependent intensities were then fitted with a quadratic limb-darkening law, (7)where I(μ)/I(0) corresponds to the normalized stellar intensity, u1 and u2 correspond respectively to the linear and quadratic limb darkening coefficients, and μ is cosθ, where θ is the angle between the normal to the stellar surface and the line of sight to the observer. In this fashion, μ = 0 refers to the limb, while μ = 1 corresponds to the center of stellar disk. To carry out the fit between the integrated stellar intensities and Eq. (7), we consider the μ values between 0.1 and 1 to avoid the steep intensity gradient appearing close to the stellar limb that, if accounted for, would have to be fitted with an added exponential growth to the quadratic function. Since the orientation of the HD 189733b orbit does not produce a nearly grazing transit, the time that the planet spends in the 0-0.1 μ-region is very short and our simplification is justified. The derived values for the sodium D2 and D1 lines in their corresponding integration and reference bands are listed in Table 2. Here it should be noted that since the integration band widths are, in this work, between 1 and 3 Å, we re-computed the PHOENIX spectra, increasing the spectral sampling to ~0.1 Å.

As previously mentioned, using the computed limb-darkening values we produce three light-curve models (Mandel & Agol 2002) and calculate the wavelength dependent limb-darkening model introduced in Eq. (6). The sensitivity of our data is not high enough to consider the limb darkening coefficients as fitting parameters. Therefore, this model is taken into account but with all their values fixed. The limb-darkening model is shown in Fig. 4B as a black line, overplotted onto the raw excess light curve.

3.2.3. Planetary radial velocity (component C)

thumbnail Fig. 5

Top: schematic shift of the planetary sodium line inside the stellar sodium line during the transit. The size of the planetary line has been exaggerated. The gray area corresponds to 1 Å integration band. Bottom: effect of the planetary line shift on the excess light curve with the integration passband of 1 Å, causing a bump that reaches its maximum level around the mid-transit point. We note that the x-axis is also proportional to the transit time, with zero set on the mid-transit point.

Open with DEXTER

Because the planet is moving during its transit on a curved orbit, its radial velocity changes. Hot Jupiters often have orbital speeds on the order of 100 km s-1, which can result in an in-and-out transit radial velocity variation larger than ±10 km s-1. This change in radial velocity causes the entire transmission of the exoplanet to be Doppler shifted, and has already been used to detect molecules in the atmospheres of hot Jupiters (e.g., Snellen et al. 2010; Brogi et al. 2012; Rodler et al. 2012, using very high-resolution spectra, R ~ 100 000). In our case as illustrated in Fig. 5-top, the sodium line originating in the planetary atmosphere shifts between the lower wing and the core of the stellar sodium line. During this shift, the planet atmosphere always absorbs a certain fraction of starlight. In other words, the relative absorption of stellar flux by the planet atmosphere is always the same. This causes the excess light curve to increase at transit center; when the stellar flux is at minimum level, the planetary sodium line is close to the center of the core of the stellar line and the radial velocity of the planet is close to 0. Near ingress and egress, the planetary sodium line is not located in the direct vicinity of the stellar line core. The planetary sodium line at these times is embedded in a part of the stellar line wing where the stellar flux is higher than in the core. Here it absorbs the same fraction of starlight, but this absorption is larger in absolute terms. This causes a “bump” in the bottom of the transit, which is illustrated in Fig. 5-bottom. This effect was first mentioned by Albrecht (2008) and can also be tentatively seen in the results of the sodium core analysis by Snellen et al. (2008), Zhou & Bayliss (2012), Wyttenbach et al. (2015), and Cauley et al. (2016). In the particular case of HD 189733b, during transit the radial velocity of the exoplanet changes between –16 km s-1 and +16 km s-1. This results in a Doppler shift of approximately ±0.31 Å around the core of the stellar sodium line. We note that with this amount of shift, the exoplanetary sodium absorption line is always embedded in the stellar sodium line of the K-type. To model the radial velocity effect, we time-averaged all out-of-transit spectra to obtain the shape of the stellar sodium lines at high signal-to-noise, and used this as the stellar sodium line template. We represented the planetary line as a Gaussian profile with depth, ANa, and width (Gaussian sigma), σNa, both considered in this work as fitting parameters. This Gaussian profile is injected into the time-averaged sodium lines calculating the expected radial velocity of the planet (VR) during each exposure (see Eq. (8)). In other words, we create 244 spectra by multiplying the moving planetary sodium line Gaussian profiles by the time-averaged normalized stellar spectrum. We call this new series of spectra the “convolved spectra”. The radial velocity of the planet at each exposure is calculated as (8)but is converted into wavelength shift to position the center of the Gaussian profile over the stellar sodium template. Here, P is the orbital period, a is the semi-major axis, t is the time of observation for each frame, and tc is the mid-transit time. The convolved spectra were then integrated following Eq. (4). The result is a numerical model that we call the “initial RV model”. In the initial RV model, even the exoplanetary Gaussian profiles taking place outside transit are also taken into account. However, in transmission spectroscopy the planetary atmosphere is invisible before ingress and after egress. To extract the time span in which the transit takes place and also to correct for the shape of the RV model during ingress and egress, the initial RV model is multiplied by a transit model that has an offset of zero and a normalized depth of one. The result gives us the final RV model. Figure 5 is a mock model that illustrates the changing radial velocity effect on the excess light curve for an integration passband of 1 Å. In the figure, the planetary and stellar sodium lines are both represented by Gaussian profiles. Finally, our RV model, as a component of the best-fit model, is plotted over the data in Fig. 4C. We note that this model is not completely symmetric since the two stellar sodium lines are not fully symmetric.

thumbnail Fig. 6

First row: raw excess light curves on sodium D2 for integration passbands of 1 Å, 1.5 Å, and 3 Å inside the core of each sodium line. The best-fit model is plotted on top in black. Second row: residuals for each passband D1.

Open with DEXTER

3.3. Final model and best-fit parameters

The overall model we use to fit the data must be a combination of the model components just introduced (A, B, C) plus a normalization constant (offset or D) for the flare component. None of the model components can reproduce the data when treated individually, thus a combination is required. We choose a model consisting of a multiplication of all of the components plus the offset. The final model is shown in Eq. (9). Since we are dealing with very small effects, the multiplication or summation of the components would both give the same results (see also Czesla et al. 2015): (9)The final model (Eq. (9)) has four fitting parameters in total: the flare scaling parameter (FLscale), the amplitude (ANa) and width (σNa) of the planetary sodium Gaussian profile, and the normalization constant (offset). Throughout this work, to explore the best-fit parameters and their associated uncertainties we perform a Markov chain Monte Carlo (MCMC) analysis, using the affine invariant ensemble sampler emcee (Foreman-Mackey et al. 2013). We employ 100 walkers, with 360 chains each, where the initial positions are synthesized from a Gaussian distribution about our best estimates. All the free parameters have uniform prior imposed. We allow a burn-in phase of ~50% of the total chain length, beyond which the MCMC is converged. The posterior probability distribution is then calculated from the latter 50% of the chain.

4. Results and discussion

4.1. Best-fitting model and errors

Table 3

Best-fit values for the model parameters obtained from the Na D1 and Na D2 excess light curves in the three integration bands.

Our best-fitting values of the MCMC procedure are summarized in Table 3. To illustrate the quality of our fit, Fig. A.1 shows the posterior distributions and the correlation plots for all the fitting parameters, computed here as an example for the 1.5 Å excess light curve around the sodium D2 line. The five remaining excess light curves produce similar plots. As the posterior distributions in Fig. A.1 show, the σNa and ANa parameters are degenerate and the rest of parameters are not correlated at all or only very weakly. However, in this work we use them to calculate the equivalent width of the sodium line (equivalently an area, nothing more than a direct multiplication of both parameters). Therefore, the degeneracy between parameters does not affect our results.

The best-fit model for sodium D2 passbands and the corresponding residuals are shown in Fig. 6. As the figure shows, the transit depth, the flare profile amplitude, and the bump strength decrease when the integration passbands increase. In smaller passbands, the planetary signal is more pronounced because there is a higher contrast between the exoplanetary absorption and stellar absorption. Moreover, in smaller integration passbands fewer points are added together and thus – compared to wider passbands where the random noise cancels out better – the light curve is more scattered.

Table 4

Statistical tests for some model combinations for the 1.5 Å passband.

To ensure that Eq. (9) is the best representation of the data, we performed some statistical tests. We made use of the Bayesian Information Criterion (Schwarz 1978), BIC = χ2 + klnN, where k is the number of model parameters and N is the number of data points. The BIC assess model fits penalized for the number of estimated parameters. We also compute chi-square (χ2) and reduced chi-square (χ2/ν), being ν the number of degrees of freedom equal to the total number of data points minus the number of fitting parameters for each model combination. Producing this statistical analysis for all the possible combinations of model components would not be very efficient since some of the possible combinations are trivially not the best-fit answer. Therefore, we exclude some of them and pre-selected five different model combinations (modules) using visual inspection as a criterion, and calculated the previously mentioned parameters for each one. Our results for the passband of 1.5 Å is shown in Table 4. Tests on other passbands also showed similar results. Based on the values of the table, the models without the flare (component A) and exoplanetary atmosphere absorbing signal (component C) are very poor. Without considering component B, the χ2, χ2/ν and BIC values did not change much and even slightly decreased. It is important to keep in mind that the differential limb-darkening is a physical effect that inevitably exists. Without this effect the exoplanetary signal must be overestimated. Therefore, although the goodness of the fit did not change much, component B was also included in defining the best model combination. Thus, briefly, from the minimization of the three statistical tools, we conclude that a good representation of the data is given by a combination of the model components A, B, and C plus the normalization constant D. The final models are introduced in Eq. (9) and the best-fit model for all of the passbands are shown in Fig. 6. As can be seen in this figure and in Fig. 4A, near the end of the observation the flare model does not perfectly fit the data. One probable reason is the shape of the flare at Ca II K lines is slightly different from the shape in sodium lines.

Furthermore, we tested whether assuming a constant model for the limb darkening effect has an impact on our results. A change in the LD model would be the product of miscalculated limb darkening values. Since there is no empirical way to test the accuracy of the LD values, we carried out the following test: any change in the limb darkening values would mostly change the amplitude of the two features centered around the mid-transit point, either by increasing or decreasing it. Therefore, we ran our whole fitting algorithm two more times, arbitrarily increasing this amplitude by 50% and decreasing it by the same amount with respect to its original shape. After computing the excess depths in the exact same fashion as explained in this work but with the LD model component changed, we observed no significant difference between different results when 1σ errors are being considered. Thus, the precision of our data does not allow for a limb darkening fitting.

thumbnail Fig. 7

Excess light curves in three integration passbands for both sodium D1 and D2 after correcting for the flare and differential limb-darkening components.

Open with DEXTER

4.2. Exoplanetary sodium line

Similar to Albrecht (2008) and Zhou & Bayliss (2012), our data clearly show a reduction of the exo-atmospheric absorption in the middle of the transit which is caused by the radial velocity shifts of the planetary sodium signal inside the broad stellar sodium line. We use this effect to model the sodium D-lines in the atmosphere of HD 189733b. According to our model parameters, we estimated the planetary sodium line depth, ANa, and a measurement of the line width, σNa. Our best-fit values, along with 1σ errors, are summarized in the last two columns of Table 3 for both planetary sodium D1 and D2 lines in the three integration bands. We find the values of the line depth and line width for the two narrower passbands to be consistent with each other within the given errors. Obviously this should be the case, since the shape of the planetary signal moving inside the stellar sodium line should be independent of the choice of passband. However, amplitudes and depths at 3 Å are slightly below the error bars of the other two passbands. As we mentioned before (in Sect. 3.1), with this data set at 3 Å the quality of the data is insufficient most probably due to the flare and short cadence of the data. In here, that must be the reason for very large error bars and underestimated planetary signal. Therefore, when calculating the average of the σNa and ANa, we exclude the results of the 3 Å data. The average line depth and width are Å and . We also estimated the equivalent width to be the product of the averaged line depth and the averaged line width, which is ~(0.0023 ± 0.001)Å.

In principle, the width of the line potentially contains the information on the line broadening sources (e.g., pressure broadening, rotational broadening, winds). However, here σNa and ANa are degenerate and therefore we cannot robustly interpret the physical conditions of the atmosphere directly from them. It is only the product of these two values that indicates a physical meaning from equivalent width of the exoplanetary sodium that blocks the star. We note that in principle it is possible to break this degeneracy by comparing the sodium line profile to the profile obtained from the complete transmission spectrum. The complete spectrum can be obtained from the division of in-transit spectra by the out-of-transit spectra, similar to the methods of Wyttenbach et al. (2015) and Redfield et al. (2008). However, we did not use that method because of very high intrinsic variability and activity of this star, especially during the time of this observation, which prevented us from achieving the required alignment between the spectra. In our upcoming paper, we will investigate this further by analyzing a planet orbiting a less active star.

4.3. Excess light-curve depths

For a better visualization of the outcome, in Fig. 7 we show all of the excess light curves after removing the flare and differential limb-darkening components. The only model plotted over the data in this figure is the planetary RV model component. In our models, we did not use any transit light-curve model directly and so the amount of the additional depth that has been added to the light curve here can be achieved by measuring the equivalent width of the exoplanetary sodium line. The additional excess to the light curve is then the equivalent width divided by the integrated flux of each passband divided by the reference band level. Table 5 shows the values of the excess light curves for the three integration bands (1 Å, 1.5 Å, and, 3 Å) as the average of sodium D2 and D1 lines. As expected, these light curves show that the depth of the excess light curve increases with the decrease of the width of the integration passband. This is due to a higher contribution of the star in the wider passbands compared to that of the planet.

Table 5

Computed values for the relative absorption depth in [%] of the light curve.

In all cases the derived excess depth obtained from the Na D2 line for a given integration band is deeper than the one obtained from the Na D1 line. This has also been found in Huitson et al. (2012). One possible scenario is explained by looking at the strength of the absorption signal, which is proportional to the contrast between the stellar sodium line and the planetary sodium line. In the case of Na D2 the stellar line is deeper and therefore the contrast is smaller. Thus, the signal in Na D2 appears larger than the transmission signal in Na D1.

4.4. Potassium and neutral calcium lines

Potassium and calcium are also among the alkali metals. After sodium, for the computed effective temperature of HD 189733b, models of exoplanet atmospheres predict potassium as a prominent absorption feature (Seager & Sasselov 2000; Fortney et al. 2010). To investigate the presence of these alkali metals in the HD 189733b atmosphere, we apply the same method to study the potassium line at 7699 Å and a pronounced neutral calcium line at 6122.22 Å. However, we were not able to detect any signature of exoplanetary potassium and calcium in their excess light curves within their noise. We attribute this to the quality of the data rather than to the chemical composition of the exoplanet atmosphere: the computed standard deviation of our potassium excess light curve is 1.2 ppt, while the amplitude of potassium derived from models is predicted to be 0.24 ppt (Pont et al. 2013). Therefore, our data are not precise enough to produce reliable results around these wavelengths. For the calcium excess light curve, the standard deviation of the excess light curve is 1.0 ppt. Although the data is slightly more accurate, in agreement with Pont et al. (2013) we did not detect any calcium absorption. This could probably mean that the quality of the data in both observations are not sufficient to detect the neutral calcium. An alternative explanation could be that the calcium abundance is not in the limit of detection.

thumbnail Fig. 8

Our averaged excess depths for each integration band centered at the sodium lines (red filled stars) compared to previous works.

Open with DEXTER

4.5. Comparison with previous work

Previous measurements of the excess absorption of the sodium in HD 189733b for passbands between 1 and 80 Å, along with our results, are shown in Fig. 8. Bibliographic values are taken from Wyttenbach et al. (2015), Huitson et al. (2012), Redfield et al. (2008), Jensen et al. (2011), and Cauley et al. (2016). Huitson et al. (2012) used low-resolution spectra of STIS/HST and obtained the excess depths for wide ranges of integration bands from 3 to 80 Å around the sodium feature (black filled diamonds). Wyttenbach et al. (2015) used high-resolution data obtained with HARPS and performed ground-based transmission spectroscopy using integration bands of 0.75, 1.5, 3, 6, and 12 Å. They used two approaches: one by fitting a transit model on the excess light curve (green filled circles) and the other by investigating the residuals of the in-transit divided by the out-of-transit spectra (in-out division; blue filled squares). Redfield et al. (2008) and Jensen et al. (2011) both used the same set of data taken by the high-resolution spectrograph mounted at the Hobby Eberly Telescope (HET) and tried to detect the exo-atmosphere by applying the in-out division approach (yellow unfilled triangle and purple cross). Finally, using the same approach, Cauley et al. (2016) analyzed a single transit observation of HiRES on Keck (unfilled brown circle). In our work we introduced a new approach for an accurate measurement of the exoplanetary transmission signal. Compared to previous measurements this approach is unbiased by stellar differential limb-darkening effect, stellar flaring activity, and the changing RV bump. Our estimated values of the additional absorption by the exoplanetary atmospheric sodium are indicated by the red filled stars in Fig. 8. In comparison to other studies, our data points show stronger absorption with larger error bars. The main reason for this difference must be the correction of the bulgy shape in the middle of the transit excess light curve caused by differential limb-darkening and changing radial velocity effects. In other words, not considering an increase in flux near the mid-transit time would underestimate the absorption depth, thus producing smaller values. The error bars in our work are larger by about a factor of 3. This is reasonable since our model has more components and therefore consists of fewer degrees of freedom.

The differential limb-darkening model calculations performed in this work can be compared with those computed by Czesla et al. (2015). The purpose of this paper is to consider the line cores to actually search for the planet signatures. Therefore, the cores of the lines are investigated, where the effects of the flare and the planets are present. Czesla et al. (2015) only considered the wings of the sodium lines where the effects of stellar activity and the planet are not present (or only marginally). Furthermore, we use ratios rather than differences, but as shown in Czesla et al., this is a small effect for small-amplitude variability (Eq. (7) in Czesla et al. 2015). Moreover, Czesla et al. (2015) use Kurucz models, while we use PHOENIX angle resolved synthetic spectra. A difference between PHOENIX and Kurucz is that the former is calculated in spherical symmetry, the latter in plane parallel, which is justified because the atmosphere in main sequence stars is thin and the curvature does not play a role. However, getting the LD correct for our needed order of precision, the plane parallel geometry is probably insufficient. In this respect, the phoenix model is considered to be better.

4.6. Rough estimation of Na number density in the HD 189733b atmosphere

Stellar abundances can be obtained from equivalent widths (EWs) of spectral lines, using a relation called the “curve of growth”. For small equivalent widths (i.e., weak lines) the equivalent width of a line is proportional to its number density. In this work we computed the equivalent width of HD 189733b as the product between the line depth and the line width, ln (EW) = ln(ANa × FWHM) = ln() = –5.02, corresponding as expected to the weak line regime. Assuming a scale height of 100 km above the clouds, the number density of sodium in the atmosphere of HD 189733b is 104 atoms/cm3, close to the value that Heng et al. (2015) estimated.

4.7. Further considerations in the model budget

Despite carrying out a meticulous analysis of the model components in this work, there are other aspects that might have a second-order impact on our data. For example, stellar spots and stellar rotation are two physical effects that deform the stellar spectral line shape and, therefore, might influence the derived excess light curves. In this work these effects were ignored for the two reasons. First, with respect to stellar rotation, HD 189733 is a slow rotating star (Winn et al. 2007, ~12 days) and since the exposures are equally distributed along the transit this effect, if any, should be averaged out in our all-transit integrations (Dravins et al. 2015; Wyttenbach et al. 2015). Second, to measure the influence of occulted and unocculted stellar spots, in addition to the rotational period and amplitude of the modulation of the total stellar flux (which would give us an idea of the spot coverage at the moment of the observations) we would have needed simultaneous broad-band observations (Pont et al. 2013; McCullough et al. 2014). Indeed, the use of simultaneous observations with low-resolution instruments is a good way to detect spot crossing events. A question that we will be able to answer with new high-resolution observations simultaneous to broad-band data is whether the change in excess depth observed in this work could have been created by unaccounted spot signal. With only the present data set we cannot assess this.

5. Conclusions

We used high-resolution spectra of UVES to measure the sodium absorption in the transmission spectrum of the hot Jupiter HD 189733b. In this work we applied a new approach based on the changing radial velocity of the exoplanet. To this end, we analyzed six excess light curves in sodium D1 and D2 that were integrated in three different bands. We modeled a combination of three main effects on each data set and extracted some information about the exo-atmospheric sodium. These effects are a stellar flare starting close to mid-transit time, the differential limb-darkening effects, and the planetary sodium line profile that is moving inside the stellar sodium line because of the change in the radial velocity of its orbit.

We confirm the ground-based detection of sodium in the upper atmosphere of HD 189733b. Close to the mid-transit time a bump appears in the excess light curves, which is caused by the radial velocity changes of the planet while moving in its orbit with respect to its parent star. By modeling this effect we estimate the equivalent width of the exoplanetary sodium line by an estimation of its width and depth. In the future, with the improvement in the quality of observations, this might be used to better constrain the physical conditions in the upper atmosphere of exoplanets as part of the atmospheric models. The average equivalent width of the two sodium lines is 0.0023 ± 0.0010 Å. The average depth of the excess light curve is 0.72 ± 0.25% at 1 Å, 0.34 ± 0.11% at 1.5 Å, and 0.20 ± 0.06% at 3 Å, in good agreement with previous research.

Finally, in this work we learn the following. First, in high-resolution transmission spectroscopy at any wavelength, Ca H & K or Hα should be always observed simultaneously when monitoring flaring activity consistency. Second, depending on some parameters, such as the temperature of the star, the width of the stellar sodium line, and the orbital velocity of the planet, the bump that appears in the mid-point of the excess light curves can originate from changing radial velocity effects and/or from the differential limb-darkening effect. A careful treatment of these two factors must be carried out. Finally, in an observation with short cadence and a large number of data points, the radial velocity bump can be better studied.

Acknowledgments

S. Khalafinejad acknowledges funding by the DFG in the framework of RTG 1351 and thanks I. Snellen, S. Albrecht, and the anonymous referee for their suggestions and advice. Accomplishment of this work was obtained through assistance, advice, and scientific discussions with many other people, especially at Hamburg Observatory and Center for Astrophysics (CfA). S. Khalafinejad would like to specifically thank M. Holman, S. Czesla, M. Güdel, B. Fuhrmeister, A. K. Dupree, M. Payne, H. Müller, P. Ioannidis, and M. Salz for their valuable support and assistance. C. von Essen acknowledges funding for the Stellar Astrophysics Centre provided by The Danish National Research Foundation (grant No. DNRF106). T. Klocová acknowledges support from RTG 1351 and DFG project CZ 222/1-1.

References

Appendix A: Additional figure

thumbnail Fig. A.1

Posterior distributions of the model parameters fitted in this work in the shape of histograms, along with their correlation plots. Values are computed from the 1.5 Å excess light curve around the sodium D2 line. Here FL is the flare scaling parameter, Sig is the width of the exoplanetary Gaussian profile, Amp is the amplitude of the exoplanetary Gaussian profile, and off is the normalization constant for the flare model.

Open with DEXTER

All Tables

Table 1

Adopted values for the orbital and physical parameters of HD 189733 during the fitting procedures in this work.

Table 2

Limb-darkening coefficients (u1, u2) for the sodium D2 and D1 integration bands (core) and fixed reference (ref.) bands with the specified wavelength ranges (λ-range).

Table 3

Best-fit values for the model parameters obtained from the Na D1 and Na D2 excess light curves in the three integration bands.

Table 4

Statistical tests for some model combinations for the 1.5 Å passband.

Table 5

Computed values for the relative absorption depth in [%] of the light curve.

All Figures

thumbnail Fig. 1

Sky conditions during our observations. The airmass is indicated with a green dashed line and the red cross points indicate the seeing, in arcsec. The measured strength of several strong telluric water lines are shown with solid blue lines. These values are artificially shifted to fit inside the plot. Gray areas show the time between transit ingress and egress. Dark gray is the transit at full depth.

Open with DEXTER
In the text
thumbnail Fig. 2

Misalignment of each spectrum with respect to a reference spectrum taken in the middle of the observing run.

Open with DEXTER
In the text
thumbnail Fig. 3

Integration bands and derived excess light curve. Left: gray shaded areas show the integration bands with passband of 1 Å, 1.5 Å, and 3 Å, centered at the core of each sodium line. The red and blue intervals show the passbands of the reference bands for D2 and D1, respectively. Right: the raw excess light curve shown with the red circles for the 1.5 Å integration band inside the sodium D2 line. Dashed vertical lines indicate times of ingress and egress.

Open with DEXTER
In the text
thumbnail Fig. 4

Individual model components plotted over the excess light curve for Na D2 at 1.5 Å. Panel A): best-fit flare model component. In this panel the green dashed lines indicate the beginning and the end of the transit; panel B): best-fit differential limb-darkening model component; panel C): best-fit changing RV model component, which includes the exoplanetary atmospheric excess absorption in addition to a bump at the center. The final best-fit model, which is the combination of all the model components, is shown in Fig. 6.

Open with DEXTER
In the text
thumbnail Fig. 5

Top: schematic shift of the planetary sodium line inside the stellar sodium line during the transit. The size of the planetary line has been exaggerated. The gray area corresponds to 1 Å integration band. Bottom: effect of the planetary line shift on the excess light curve with the integration passband of 1 Å, causing a bump that reaches its maximum level around the mid-transit point. We note that the x-axis is also proportional to the transit time, with zero set on the mid-transit point.

Open with DEXTER
In the text
thumbnail Fig. 6

First row: raw excess light curves on sodium D2 for integration passbands of 1 Å, 1.5 Å, and 3 Å inside the core of each sodium line. The best-fit model is plotted on top in black. Second row: residuals for each passband D1.

Open with DEXTER
In the text
thumbnail Fig. 7

Excess light curves in three integration passbands for both sodium D1 and D2 after correcting for the flare and differential limb-darkening components.

Open with DEXTER
In the text
thumbnail Fig. 8

Our averaged excess depths for each integration band centered at the sodium lines (red filled stars) compared to previous works.

Open with DEXTER
In the text
thumbnail Fig. A.1

Posterior distributions of the model parameters fitted in this work in the shape of histograms, along with their correlation plots. Values are computed from the 1.5 Å excess light curve around the sodium D2 line. Here FL is the flare scaling parameter, Sig is the width of the exoplanetary Gaussian profile, Amp is the amplitude of the exoplanetary Gaussian profile, and off is the normalization constant for the flare model.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.