Confirmation of Asymmetric Iron Absorption in WASP-76b with HARPS

In recent years, high-dispersion spectroscopy has transformed the field of exoplanet atmospheric characterization by resolving signatures of atmospheric dynamics in hot Jupiters. Snellen et al. (2010) first observed a possible net blueshift of CO absorption on HD 209458b and attributed it to a global wind moving from the hot dayside to the cooler nightside of the planet. These high-velocity winds were predicted by 3D global circulation models and are thought to homogenize the temperatures on the daysides and nightsides on these strongly irradiated gas giants (Showman et al. 2008). Since the first detection, global winds have been detected on many hot Jupiters with velocities up to 10 km s⁻¹ (e.g., Alonso-Floriano et al. 2019; Casasayas-Barris et al. 2019; Nugroho et al. 2020).

Ehrenreich et al. (2020) further advanced the use of high-dispersion spectroscopy to study atmospheric dynamics by resolving an asymmetric absorption signature in two individual transits of WASP-76b (West et al. 2016). This was possible for the first time thanks to the radial velocity precision of ESPRESSO on ESO's 8 m Very Large Telescope (see also Tabernero et al. 2020 for a subsequent analysis of the data set, cataloging many atoms). Ehrenreich et al. (2020) find that while at the start of the transit an iron absorption signal is visible near the planet's rest-frame velocity, at the end of the transit it is significantly blueshifted to −11 km s⁻¹. The cause of this signature was explained by a combination of planetary rotation and a global wind coming from the dayside, necessitating that the iron condenses out on the planet nightside. During a transit, the angle under which the planet is observed changes by almost 30°, meaning that the irradiated hemisphere on the morning side is visible early on during the transit, while the evening side is only seen during the second half of the transit. Assuming that there is no or very little iron absorption from the nightside of the planet owing to rainout, the components of iron absorption from the morning and evening side rotate into and out of view, showing different levels of radial velocity shift due to the rotation and global wind counteracting and reinforcing each other, respectively (Ehrenreich et al. 2020). This is the first time that a chemical gradient has been directly observed across the surface of an exoplanet, and it demonstrates the future potential of high-dispersion spectroscopy to observationally constrain condensation and 3D circulation modeling.

Evidence for a similar asymmetry in iron has been seen in WASP-121b and MASCARA-2b as well, which further solidifies the presence of condensation signatures in many ultrahot Jupiters. Bourrier et al. (2020) find a tentative offset in WASP-121b's systemic velocity derived from the entire cross-correlation mask (which is dominated by iron lines) from the beginning to the end of the transit. Using ESPRESSO, Borsa et al. (2021) validated the presence of the blueshift variations in the iron signal during the transit of WASP-121b, finding a blueshift of −2.80 ± 0.28 km s⁻¹ at the beginning of the transit and −7.66 ± 0.16 km s⁻¹ at the end. Finally, Hoeijmakers et al. (2020a) find evidence of a strengthening absorption signature from iron in MASCARA-2b from the beginning to the end of the transit, which could also indicate changing chemistry across the surface of the planet.

In this paper, we aim to confirm the detection of the iron signature in WASP-76b using HARPS archival data. In Section 2 we present the data, and in Section 3 we explain how we created the model iron spectrum. In Section 4 we describe the cross-correlation analysis between the data and the model. Next, we present the results of the iron condensation feature in Section 5. Finally, we summarize our conclusions in Section 6.

WASP-76 is a bright F7-type main-sequence star. The relevant parameters of the system are given in Table 1. The data utilized in this work consist of time-series spectra taken during four different transits of WASP-76b. All of the spectra were obtained with the High Accuracy Radial velocity Planet Searcher (HARPS; Mayor et al. 2003) at ESO's La Silla 3.6 m telescope. HARPS has a spectral resolving power of ∼115,000 and covers a wavelength range from 378 to 691 nm. All of the data are publicly available and were downloaded from the ESO archive.¹ Table 2 lists the observing details of the four separate nights. The transit on night 1 was published in Brown et al. (2017) and was used to derive an estimate of the Rossiter–McLaughlin effect. The transits on nights 2 and 3 were part of the Hot Exoplanet Atmospheres Resolved with Transit Spectroscopy (HEARTS) survey and were previously published in Seidel et al. (2019) as a focused study on the sodium doublet. Night 4 has not been previously published.

Table 1.  WASP-76 System Parameters

Parameter	Value	Reference
Star
Spectral type	F7	W16a
V-mag	9.52 ± 0.03	CDS SIMBAD
M* (MSun)	1.458 ± 0.021	E20b
R* (RSun)	1.756 ± 0.071	E20
Teff (star; K)	6250 ± 100	W16
vsys (km s−1)	−1.16	E20
$v\sin i$ (km s−1)	1.48 ± 0.28	E20
K* (m s−1)	${116.02}_{-1.35}^{+1.29}$	E20
Planet
Mp (MJ)	${0.894}_{-0.013}^{+0.014}$	E20
Rp (RJ)	${1.863}_{-0.083}^{+0.070}$	E20
log g	2.85	calculated
Teq (K)	2160 ± 40	W16
Porb (d)	1.809886	W16
T0 (d)	2456107.85507	W16
a (au)	0.0330 ± 0.0005	W16
Kp (km s−1)	196.52 ± 0.94	E20
i (deg)	${88.0}_{-1.6}^{+1.3}$	W16
vrot (km s−1)	5.1	calculated

Notes.

^aWest et al. (2016). ^bEhrenreich et al. (2020).

Download table as:  ASCIITypeset image

The data were processed with the HARPS Data Reduction Pipeline (DRS v3.5 for night 1 and v3.8 for nights 2, 3, and 4), which performs flat-fielding, 2D extraction, and wavelength calibration. The pipeline produces 1D blaze-corrected, stitched spectra. The wavelength solutions have already been corrected for the barycentric velocity and are given in air. The spectra were subsequently shifted to the star's rest frame using the known system velocity (v_sys) and corrected for the reflex motion of the star due to the planet (K_*) at each observation. The spectra were then interpolated onto the same wavelength grid (of length 308,100 pixels), masking any missing values at the wavelength extremes due to the different velocity shifts to create a 308,100 × N wavelength grid in time for each night.

Even though the blaze functions have been corrected by the pipeline, small differences in the slope of the stellar spectra due to imperfect blaze removal can cause problems when the stellar spectrum is removed in a later stage of the analysis. As we do not want to remove the shape of the spectrum, since it contains information about the response function of the spectrograph, transmission through Earth's atmosphere, and the signal from the host star, we followed a similar procedure to Hoeijmakers et al. (2019) and divided the 1D spectra into 50 parts, each of size 6162 × N, and normalized each wavelength section to the mean flux in time.

To remove any contamination from cosmic rays, we performed a 3σ clipping procedure. We flagged any values that were more than 3σ different from the other time-series spectra at each wavelength pixel and interpolated over them in the wavelength direction. We also performed a masking of entire wavelength columns that had high standard deviations owing to telluric lines, bad pixels, or deep stellar lines. We tested masking different percentages of pixels and chose the value that most decreased the noise in the final cross-correlation function without significantly altering the peak. This value ranged from about 5% (nights 2, 3, and 4) to 15% (night 1) of all pixels in the four different nights.

We created a model transmission spectrum of iron using petitRADTRANS (Mollière et al. 2019), in a method similar to Kesseli et al. (2020). petitRADTRANS is a radiative transfer code that creates emission or transmission spectra at low or high resolution. The high-resolution mode includes accurate line opacities for a large number of atoms, ions, and molecules. We also tested using the same F9 radial velocity stellar mask as was used in Ehrenreich et al. (2020), but we found a stronger signal from the planet when we used a model atmosphere specifically tuned to the parameters of the planet (signal-to-noise ratio (S/N) of final detection increased from ∼5 to 9).

petitRADTRANS requires as input a temperature–pressure profile, the planet's surface gravity and radius, the abundances of the requested species, and the mean molecular weight (MMW) of the atmosphere. We used the planet radius and surface gravity that are listed in Table 1, as well as an isothermal temperature–pressure profile at the equilibrium temperature of the planet. We experimented with different temperatures and found that changes of 150 K could change the height of the retrieved signal by small amounts (∼50 ppm) but did not change the velocity position. We implemented an equilibrium chemistry code from Mollière et al. (2017) to determine the MMW of the planet at each altitude and assumed a constant volume mixing ratio (VMR) of 10⁻⁵ for Fe at all altitudes. While a constant VMR at each altitude is not realistic, neither this assumption nor an order-of-magnitude change to the VMR significantly alter our recovered planet signal. A VMR of 10⁻⁵ is motivated by previous chemical modeling of hot Jupiters (Lodders 2002; Visscher et al. 2010; Hoeijmakers et al. 2020b). Finally, we added continuum opacity from H₂ − H₂ collisions, H₂ − He collisions, and H⁻.

We reduced the resolution of the model to the resolution of the HARPS spectrograph by convolving it with a Gaussian kernel with an FWHM of 2.7 km s⁻¹. We did not perform any steps to take into account the planet's rotation, as the asymmetry in the iron signal is thought to be due to iron being present mostly on one limb of the planet.

We followed the steps of Hoeijmakers et al. (2020b) for the cross-correlation analysis and used the following equation:

$$\begin{eqnarray}&&c(\nu ,t)=\displaystyle \sum _{i=0}^{M}{x}_{i}(t){T}_{i}(\nu ),\end{eqnarray} \tag{ 1 }$$

where c(ν, t) is the resulting cross-correlation grid for each radial velocity value (ν) at each time (t). We define x_i(t) as the observed spectrum at each time and T_i(ν) as the model at each velocity. The model was shifted from −250 to +250 km s⁻¹ in 1 km s⁻¹ steps and then interpolated onto the same wavelength grid as the data. We normalized the model such that ${\sum }_{i=0}^{M}{T}_{i}(\nu )=1$. In this way, the cross-correlation function sums all of the pixels (i) in the spectrum, weighting them by their absorption.

At this point the contribution from the host star still needs to be removed, which dominates the cross-correlation function. In order to remove the host star, we divided each cross-correlation function by the average out-of-transit cross-correlation function. After this step, we were left with the residual amplitude, which corresponds to the transit depth due to atmospheric Fe.

The cross-correlation matrices for each night are shown in Figure 1. Each night has an excess along the exoplanet's expected velocity, but in all of the individual nights there are significant residuals around 0 km s⁻¹ owing to a combination of the Doppler shadow and reduced S/N in the iron line cores. By combining all of the transits together, we reduced much of this noise and were able to more easily determine the shape of the trace. To combine the nights together, we interpolated the cross-correlation grids onto a new cross-correlation grid that was uniformly sampled along the phase axis in steps of 0.004. We chose this step size to be slightly larger than the average phase step size of each night, which ranged from about 0.002 to 0.0038, so that none of the values would be extrapolated. We then averaged the four nights together, weighting each by its average S/N and the number of spectra during the transit (Table 2). The left panel of Figure 2 shows the averaged cross-correlation matrix from all four nights. The exoplanet's trace can now clearly be seen.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. The Doppler Shadow

To determine whether the shape of the iron absorption is indeed asymmetric during the transit, we shifted the exoplanet into its expected rest frame using the value of K_p from Table 1 (panel (d) of Figure 2). Between phases of 0.0 and +0.04 (the second half of the transit) all of the cross-correlation functions peak around −8 km s⁻¹, except for two noisy pixels at phases 0.028 and 0.022. Before a phase of −0.01 the peaks shift towards 0 km s⁻¹, and at the very beginning of the transit (phase = −0.04) the peak of the cross-correlation function is at +1.5 km s⁻¹.

We combined the 2D cross-correlation function grid into a single cross-correlation function by averaging the phases during the transit together (Figure 3). The peak residual amplitude of the feature is at about 500 ppm. If we average the beginning of the transit together and the end of the transit together separately, we find that the beginning of the transit has a peak radial velocity of +1.5 km s⁻¹, while the end of the transit has a peak at −7.5 km s⁻¹. The velocities of this asymmetry are similar to the ones found in Ehrenreich et al. (2020).

Zoom In Zoom Out Reset image size

We calculated uncertainties on these radial velocity shifts by fitting Gaussians to the 1D cross-correlation function of the beginning of the transit and the end of the transit separately (red and blue curves shown in the right panel of Figure 3). The Gaussian for the beginning of the transit has a center at +0.04 ± 1.47 km s⁻¹, while the Gaussian for the end of the transit has a center at −6.9 ± 0.75 km s⁻¹. The stated uncertainties in the measurement are calculated by dividing the standard deviation of the Gaussian by the S/N of the line. As a secondary test, we also checked these values by taking the average and standard deviation of the individual values in each each night at the beginning and end of the transits. For each night, we recorded the value of any peak present between −25 and +10 km s⁻¹ for each cross-correlation function between phases +0.042 and +0.002 for the second half of the transit and −0.026 and −0.042 for the first half of the transit. The average radial velocity in the first half of the transit was +0.3 ± 1.07, while the average radial velocity in the second half of the transit was −7.1 ± 0.74 km s⁻¹. The uncertainties in the peak are given by the standard deviation divided by the square root of the number of measurements. Both methods return very similar results and show that the peaks are not consistent within two standard deviations, and so we consider this a confirmation of the asymmetric iron absorption feature.

We find that the cross-correlation peak during the second half of the transit is stronger than the peak during the beginning of the transit by about 100 ppm. This is consistent with the results from Ehrenreich et al. (2020) and aids in the hypothesis that monatomic iron condenses out of the atmosphere on the cooler morning side of the planet and on the even colder nightside. The iron signal from the four nights of HARPS data presented here and the iron signal presented in Ehrenreich et al. (2020) both have noise levels of about 50 ppm pixel⁻¹. By taking into account this noise level, the wavelength steps and resolutions of the two spectrographs, and the signal amplitude and assuming unresolved iron lines, we find that 7.3 HARPS transits are required to achieve similar results to one ESPRESSO transit.

We analyzed four publicly available transits of WASP-76b from the HARPS instrument on ESO's La Silla 3.6 m telescope. To search for atomic absorption from Fe, we created a model of the transmission spectrum using petitRADTRANS. We then cross-correlated the model with the data using a method similar to Hoeijmakers et al. (2020b). We resolved a small Doppler shadow signature and removed it by modeling the stellar spectra at each phase.

We found clear evidence of absorption from Fe at an S/N of ∼9 and an average residual amplitude of 500 ppm across the transit. Our analysis shows an asymmetric velocity shift from the morning to the evening side of the planet. We divided the transit into phase bins and determined the radial velocity of the peak and uncertainty of this peak for the beginning and end of the transit separately. At the beginning of the transit, the signal is best fit by a Gaussian with a peak of +0.04 ± 1.47 km s⁻¹, while during the second half of the transit the peak is at −6.9 ± 0.75 km s⁻¹. The morning and evening radial velocities that we find are not within two standard deviations of each other, further confirming the results from Ehrenreich et al. (2020).

Additional chemical modeling and observations are required to determine the process that is removing iron from the gas phase on the nightside of WASP-76b (i.e., condensation into metallic iron or other iron molecules). FeH is predicted to be the second most abundant iron carrier at this temperature after monatomic Fe (Visscher et al. 2010), and so some Fe gas could condense into FeH gas, but Kesseli et al. (2020) did not find any evidence of FeH in WASP-76b. More observations of Fe condensation and other Fe molecules on a range of exoplanets will aid in the physical interpretation of this signature.

A.K. and I.S. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program under grant agreement No. 694513. This paper is based on data obtained from the ESO Science Archive Facility under request No. 569981 by Aurora Kesseli. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

Facility: ESO:3.6 m (HARPS). -

Software: astropy (Astropy Collaboration et al. 2018), matplotlib (Hunter 2007), numpy (Van Der Walt et al. 2011).