ACCESS: A Visual to Near-infrared Spectrum of the Hot Jupiter WASP-43b with Evidence of H₂O, but No Evidence of Na or K

Observations of exoplanetary atmospheres offer the possibility of understanding the atmospheric physical properties and chemical composition of those worlds, as well as providing clues to their formation and evolution histories (e.g., Moses et al. 2013; Öberg et al. 2013; Mordasini et al. 2016; Espinoza et al. 2017). The first comparison studies of exoplanetary atmospheres using transmission spectra (see, e.g., Sing et al. 2016; Crossfield & Kreidberg 2017) found evidence of a gradual transition between clear and cloudy atmospheres, but no clear correlation of that transition with other system parameters, such as planetary mass, gravity, effective temperature, or stellar irradiation levels, and the chemical composition of the star. Recently, Pinhas et al. (2019) reanalyzed the Sing et al. (2016) sample and concluded that the majority of hot Jupiters have atmospheres consistent with subsolar H₂O abundances, with the log of those values ranging from $-{5.04}_{-0.30}^{+0.46}$ to $-{3.16}_{-0.69}^{+0.66}$.

High-altitude clouds/hazes have been inferred in the atmosphere of a number of exoplanets from scattering slopes in the visible (e.g., Pont et al. 2008, 2013; Sing et al. 2009, 2011, 2013; Gibson et al. 2013) and from the damping of pressure-broadened alkali Na i and K i lines originating from deeper within the atmosphere (e.g., Charbonneau et al. 2002; Wakeford & Sing 2015). The first detections of exoplanets with potentially clear atmospheres have only recently been made, e.g., WASP-96b (Nikolov et al. 2018).

This field is currently at the point where a number of efforts are underway to identify what system parameters, if any, correlate with the observed atmospheric properties of exoplanets. For example, the relationship between chemical abundance in an exoplanet's atmosphere and planet mass is actively being explored (Fraine et al. 2014; Kreidberg et al. 2014; Helling et al. 2016; Sing et al. 2016). Ground-based (e.g., ACCESS,¹⁵ GPIES,¹⁶ VLT FORS2,¹⁷ LRG-BEASTS¹⁸ ) and space-based surveys (PanCET¹⁹ ) are working to provide homogeneous spectra of a statistically significant number of exoplanet atmospheres in the search for these correlations. In this paper, we present the ground-based visual to near-infrared (NIR) transmission spectrum of the hot Jupiter WASP-43b obtained as part of the ACCESS survey.

WASP-43b (M_p = 2.052 ± 0.053 M_J, R_p = 1.036 ± 0.012 R_J, ${T}_{\mathrm{eq}}={1440}_{-39}^{+40}$ K; Gillon et al. 2012) is a hot Jupiter discovered by Hellier et al. (2011) transiting a V = 12.4 K7V-type dwarf star (M_s = 0.717 ± 0.025 M_⊙, R_s = 0.667 ± 0.010 R_⊙, T_eff = 4520 ± 120 K; Gillon et al. 2012) every 0.81 days. WASP-43 is unusually active for its stellar type, as indicated by the presence of strong Ca H and K lines and perhaps due to star–planet interactions in this very short-period system (Staab et al. 2017). Spitzer secondary eclipse data in the 3.6 μm and 4.5 μm bands indicate brightness temperatures of 1670 ± 23 K and 1514 ± 25 K, respectively, which rule out a strong thermal inversion in the planet's dayside photosphere (Blecic et al. 2014). In addition, thermal emission observed in the K band (Chen et al. 2014) agrees with atmospheric models of WASP-43b, which predict poor day-to-night heat redistribution in an atmosphere with no thermal inversion present (Kataria et al. 2015).

The presence of water on the dayside of the planet was observed with Hubble Space Telescope (HST)/WFC3 emission measurements by Stevenson et al. (2014), with additional transmission observations by Kreidberg et al. (2014) finding water abundances comparable to solar values. In the visual regime, the Gran Telescopio Canarias's visual System for Imaging and low Resolution Integrated Spectroscopy (OSIRIS) instrument shows a tentative excess in R_p/R_s at Na i and a complete lack of one near the K i doublet (Murgas et al. 2014). That same study also notes an increasing trend of the planet-to-star radius ratio from 6200 to 7200 Å and a decreasing trend redward of 7200 Å. They attribute this pattern to the possible presence of VO and TiO.

In this work, we search further for the presence of Na i, K i, and Hα with new visual transit observations from the Magellan/Inamori-Magellan Areal Camera and Spectrograph (IMACS). In addition, we combine the HST/WFC3 transmission spectrum of WASP-43b (Kreidberg et al. 2014) with the new visual data to produce the full visual to NIR spectrum spanning 5317.90 Å–16420 Å that can further constrain the water absorption features present in the infrared spectrum and provide new information about water abundance. In our analysis, we find that the atmosphere of WASP-43b is best described by a clear atmosphere with water abundance consistent with solar. The planet's spectrum is also contaminated with stellar heterogeneity.

This paper is structured as follows: in Section 2, we present our Magellan/IMACS observations. In Section 3, we outline the data reduction process used in our observations and describe the selection of wavelength bins to optimize the signal-to-noise ratio (S/N), search for atomic features, and compare to other results. We present the detrended white and binned light curves for each data set. In Section 4, we give a qualitative analysis of the impact of observational stellar activity on the resulting combined transmission spectra from different visits. In Section 5, we present the final transmission spectrum and also compare to the results of Murgas et al. (2014). In Section 5.4, we present the results of a retrieval modeling analysis on the combined ACCESS and HST transmission spectrum to find the best-fit transmission model when the presence of a heterogeneous stellar photosphere is also taken into account. We summarize and conclude in Section 7.

2.1. General Setup

2.2. Data Collection

3.1. Reduction Pipeline

We reduced the raw data using the ACCESS pipeline described previously (Rackham et al. 2017; Bixel et al. 2019; Espinoza et al. 2019). The detailed functions of the pipeline, including standard bias and flat calibration, bad pixel and cosmic-ray correction, sky subtraction, spectrum extraction, and wavelength calibration are described in detail in Jordán et al. (2013) and Rackham et al. (2017). We briefly summarize the data reduction here.

We applied the wavelength solution found with the arc lamps to the first science image, and the remaining science image spectra were then cross-correlated with the first's. We calibrated the spectra from all stars to the same reference frame by identifying shifts between the Hα absorption line minimum of the median spectra and air wavelength of Hα, and interpolated the spectra onto a common wavelength grid using b-splines. We aligned all spectra to within 2 Å, which is less than the average resolution element of 4.7 Å, assuming an average resolution of R = 1200 and shortest wavelength coverage of 5600 Å on IMACS. We also subtract scattered light within IMACS along each slit for every exposure, as described by Espinoza et al. (2019).

The final results are sets of wavelength-calibrated and extracted spectra for the target and each comparison star that can be used to produce integrated (white light; Figure 1) or spectroscopically binned light curves. The series of white-light curves (WLCs) produced in this fashion informed which comparison stars to omit in the rest of the analysis on a per data set basis. Based on the deviations of each comparison star's flux from the general trend of WASP-43b's flux in Figure 1, we omitted comparison star 5 from the first 2015 data set (ut150224, Transit 1), comparison stars 5 and 6 from the second 2015 data set (ut150309, Transit 2), comparison star 5 from the 2017 data set (ut170410, Transit 3), and no comparison stars from the 2018 data set (ut180603, Transit 4).

Zoom In Zoom Out Reset image size

With this set of good comparison stars established, we made a normalized white-light curve (F_div) for each observation by dividing the target WASP-43b flux F_W43 by the sum of the N good comparison star fluxes $({F}_{{C}_{n}})$,

$$\begin{eqnarray}&&{F}_{\mathrm{div}}=\displaystyle \frac{{F}_{{\rm{W}}43}}{{\displaystyle \sum }_{n}^{N}{F}_{{C}_{n}}},\end{eqnarray} \tag{ 1 }$$

to produce the curves shown in Figure 2 for each data set. We omitted the remaining outliers (points with fluxes deviating at least 2σ from the median flux of 10 neighboring points) from the rest of our detrending procedures. Because the different wavelength bin schemes we explored (see Section 3.2.4) sometimes had different outlier points in their respective divided white-light curves, we omitted a superset of outlier points from all data sets to maintain uniformity and reduce systematics.

Zoom In Zoom Out Reset image size

We detrended our original data, with the selected data points and comparison stars for each transit epoch selected, using three different methods described in Section 3.2.

3.2. Light-curve Analysis

We applied three detrending and transit-fitting methods to our white light and wavelength-binned data to test our results and verify that they were not dependent on the method used. The three detrending methods we used were polynomial wavelet (poly), polynomial wavelet followed by a common-mode correction (poly+CMC), and a Gaussian process combined with principal component analysis (GP+PCA).

3.2.1. Polynomial Wavelet Detrending

We performed a simultaneous transit model and systematics detrending fit on each of the transit white-light curves shown in Figure 2 using our Markov Chain Monte Carlo (MCMC) code, described in detail in Rackham et al. (2017). To briefly summarize, the divided light-curve model $F(t,{\boldsymbol{\theta }})$ can be written as

$$\begin{eqnarray}&&F(t,{\boldsymbol{\theta }},{\boldsymbol{\alpha }})=f({\boldsymbol{\theta }})P(t,{\boldsymbol{\alpha }}),\end{eqnarray} \tag{ 2 }$$

where t is the time, $f({\boldsymbol{\theta }})$ is the analytic transit model described in Mandel & Agol (2002), ${\boldsymbol{\theta }}$ is the vector of the orbital and transit parameters (e, a/R_s, i, ω, T, t₀, b, u₁, u₂, R_p/R_s), defined in Table 2, and ${\boldsymbol{\alpha }}$ is a vector of polynomial coefficients (α₀, α₁, ..., α_m), where

$$\begin{eqnarray}&&P(t,{\boldsymbol{\alpha }})=\displaystyle \sum _{m=0}^{M}{\alpha }_{m}{t}^{m}\end{eqnarray} \tag{ 3 }$$

for an Mth-order polynomial that we fit to the out-of transit (OOT) flux of F_div. We assumed a quadratic stellar limb-darkening profile and sampled the limb-darkening coefficients u₁ and u₂ according to Kipping (2013) to allow uninformative (uniform) priors to be placed on the transformed coefficients while avoiding the risk of sampling nonphysical values. We fit the divided light-curve model to determine the most likely values for ${\boldsymbol{\alpha }}$. We obtained the final detrended model (F_det) by dividing through by the OOT best-fit flux model $P(t,{\boldsymbol{\alpha }})$, such that

$$\begin{eqnarray}&&{F}_{\det }=\displaystyle \frac{{F}_{\mathrm{div}}}{P(t,{\boldsymbol{\alpha }})}.\end{eqnarray} \tag{ 4 }$$

Table 2.  Literature System Parameters

Bin Scheme	1	2	3
Eccentricity (e)	0	0	0
Semimajor axis/stellar radius (a/Rs)	4.867 ± 0.023	4.752 ± 0.066	4.872
Inclination (radians) (i)	1.426 ± 0.0056	1.433 ± 0.00175	1.433
Planet/stellar radius uncertainty $({\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}})$	0.0018	0.00145	0.00043a
Longitude of periastron (ω)	π	π	π
Period (days) (T)	0.813473978 ± 3.5 × 10−8	0.81347459 ± 2.1 × 10−7	0.81347436
Reference	Hoyer et al. (2016)	Murgas et al. (2014)	Kreidberg et al. (2014)

Notes. Literature system parameters corresponding to each bin scheme.

^aWe used ${\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$ from Hoyer et al. (2016) combined transit data because they are based on ground-based values while values from Kreidberg et al. (2014) are space based.

Download table as:  ASCIITypeset image

We performed a simultaneous MCMC fitting of transit parameters with PyMC (Salvatier et al. 2016). The likelihood was determined through the wavelet method described in Carter & Winn (2009). The fitted parameters include the midtransit time (t₀), planet-to-star radius ratio (R_p/R_s), three coefficients (a_i, i = 0, 1, 2) for the second-order polynomial used to fit the baseline OOT trend, two parameters for the transformed quadratic limb-darkening coefficients (q₁, q₂), one noise parameter (σ_w) for uncorrelated "white" noise, and one parameter (σ_r) for correlated "red" noise. We found that a second-order polynomial fit the OOT flux better than a lower order function.

We sampled t₀, σ_w, σ_r, q₁, and q₂ with uniform priors; α_i with a Gaussian prior with a width set by bootstrapping the uncertainty on α_i following Rackham et al. (2017); and R_p/R_s with a Gaussian prior with spread $5\times {\sigma }_{{R}_{{\rm{p}}}/{R}_{{\rm{s}}}}$, where R_p/R_s and ${{\boldsymbol{\sigma }}}_{{{\boldsymbol{R}}}_{{\rm{p}}}/{{\boldsymbol{R}}}_{{\rm{s}}}}$ are taken from the literature. We used five chains, each composed of 100,000 steps, and started them at the estimated location of maximum a posteriori probability, using an additional 30,000 steps for burn in. We thinned the chains by sampling them at 10× their autocorrelation function half-life before combining them to produce the final posterior distributions. Each bin scheme has a corresponding source for the fixed system parameters, and we list them in Table 2. We fit the white-light curve while keeping the midtransit time t₀ and average transit depth R_p/R_s free.

We produced the binned light curves following the same procedure for producing the WLCs, with the only difference being that we kept the midtransit time t₀ found in the WLC analysis fixed when performing the fits for R_p/R_s in each wavelength bin for bin schemes 1, 2, and 3, described in Section 3.2.4.

3.2.2. Common-mode Correction

Following Sedaghati et al. (2016), we divided the detrended WLC obtained in the polynomial detrending method by its best-fit model to produce a common-mode correction (CMC) residual. We then divided that residual through the binned light curves to remove wavelength-independent variations. Next, we applied the CMC correction to each bin scheme.

3.2.3. Gaussian Process and Principal Component Analysis

Gaussian Process (GP) regression is a powerful tool for modeling data in the machine-learning community (Rasmussen & Williams 2005) that has started to gain more popularity in the exoplanets field (see, e.g., Aigrain et al. 2012; Gibson et al. 2012). Gibson et al. (2012) provide a good overview to this methodology applied to exoplanet transit light curves. Applying this methodology for a collection of N measurements $({\boldsymbol{f}})$, such as the flux of a star measured over a time series, the log marginal likelihood of the data can be written as

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{log}{ \mathcal L }({\boldsymbol{r}}| {\mathsf{X}},{\boldsymbol{\theta }},{\boldsymbol{\phi }}) & = & -\displaystyle \frac{1}{2}{{\boldsymbol{r}}}^{\top }{{\rm{\Sigma }}}^{-1}{\boldsymbol{r}}-\displaystyle \frac{1}{2}\mathrm{log}\left|{\rm{\Sigma }}\right|\\ & & -\,\displaystyle \frac{N}{2}\mathrm{log}(2\pi ),\end{array}\end{eqnarray} \tag{ 5 }$$

where ${\boldsymbol{r}}\equiv {\boldsymbol{f}}-T({\boldsymbol{t}},{\boldsymbol{\phi }})$ is the vector of residuals between the data and analytic transit function T; ${\mathsf{X}}$ is the N × K matrix for K additional parameters, where each row is the vector of measurements ${{\boldsymbol{x}}}_{n}=({x}_{n,1},\,\cdots \,{x}_{n,K})$ at a given time n; ${\boldsymbol{\theta }}$ are the hyperparameters of the GP; ${\boldsymbol{\phi }}$ are the transit model parameters; and Σ is the covariance of the joint probability distribution of the set of observations ${\boldsymbol{f}}$. In our analysis, we used six systematics parameters: time, FWHM of the spectra on the CCD, air mass, position of the pixel trace through each spectra on the chip of the CCD, sky flux, and shift in wavelength space of the trace. We used the Python package batman (Kreidberg 2015) to generate our analytic transit model. From here, the log posterior distribution $\mathrm{log}{ \mathcal P }({\boldsymbol{\theta }},{\boldsymbol{\phi }}| {\boldsymbol{f}},{\mathsf{X}})$ can be determined by placing explicit priors on the maximum covariance hyperparameters and the scale-length hyperparameters. From ${ \mathcal P }$, the transit parameters can then be inferred by optimizing with respect to ${\boldsymbol{\theta }}$ and ${\boldsymbol{\phi }}$.

We accomplished the above optimization problem with the Bayesian inference tool PyMultiNest (Buchner et al. 2014) and computed the log likelihoods from the GP with the george (Ambikasaran et al. 2014) package. We implemented this detrending scheme by simultaneously fitting the data with an exponential squared kernel for the GP under the assumption that points closer to each other are more correlated than points farther apart. The PCA methodology follows from Jordán et al. (2013) and Espinoza et al. (2019), where M signals, S_i(t), can be extracted from M comparison stars and linearly reconstructed according to the eigenvalues, λ_i, of each signal. This allows for the optimal extraction of information from each comparison star to inform how the total flux of WASP-43 varies over the course of the night. We Bayesian model averaged the principal components together, which were determined by fitting with one, then two, up to M principal components, to create the final detrended WLC and model parameters of interest. We present the associated WLCs in Figure 3, the best-fit parameters in Table 3, and the associated corner plots in the Appendix. Based on the quality of the fits, discussed in Section 5.4, and to streamline our work, we show only the results of this method for the transmission spectra that informed our retrieval analysis.

Zoom In Zoom Out Reset image size

Table 3.  Fitted WLC Values from GP+PCA Detrending Method Shown in Figure 3

Parameter	Definition	Transit 1	Transit 2	Transit 3	Transit 4
Rp/Rs	Planet radius/star radius	${0.15854}_{-0.00074}^{+0.00079}$	${0.15800}_{-0.00278}^{+0.00218}$	${0.15436}_{-0.00164}^{+0.00167}$	${0.16030}_{-0.00118}^{+0.00127}$
d	Transit depth (ppm)	$25134\pm 250$	$24963\pm 880$	$23828\pm 516$	$25695\pm 406$
t0 − 2450000	Midtransit (JD)	${7077.72325}_{-0.00004}^{+0.00004}$	${7090.73888}_{-0.00008}^{+0.00008}$	${7854.59100}_{-0.00008}^{+0.00008}$	${8273.53019}_{-0.00005}^{+0.00005}$
P	Period (days)	${0.81347}_{-0.00000}^{+0.00000}$	${0.81347}_{-0.00000}^{+0.00000}$	${0.81347}_{-0.00000}^{+0.00000}$	${0.81347}_{-0.00000}^{+0.00000}$
a/Rs	Semimajor axis/star radius	${4.92738}_{-0.02850}^{+0.02856}$	${4.97737}_{-0.06297}^{+0.05644}$	${4.90935}_{-0.05815}^{+0.05881}$	${4.85836}_{-0.03367}^{+0.03458}$
b	Impact parameter	${0.65644}_{-0.00720}^{+0.00698}$	${0.65270}_{-0.01372}^{+0.01433}$	${0.67153}_{-0.01576}^{+0.01335}$	${0.66218}_{-0.00818}^{+0.00822}$
i	Inclination	${82.34509}_{-0.11874}^{+0.11731}$	${82.46517}_{-0.25848}^{+0.23358}$	${82.13717}_{-0.24363}^{+0.26874}$	${82.16642}_{-0.14579}^{+0.14519}$
q1	LD coeff 1	${0.70515}_{-0.23016}^{+0.19778}$	${0.60216}_{-0.17155}^{+0.23895}$	${0.50975}_{-0.19881}^{+0.30047}$	${0.72406}_{-0.19668}^{+0.18106}$
q2	LD coeff 2	${0.38994}_{-0.19375}^{+0.10395}$	${0.30680}_{-0.18818}^{+0.17608}$	${0.41212}_{-0.24787}^{+0.21415}$	${0.31602}_{-0.16639}^{+0.10556}$

Notes. We share the associated corner plots in the Appendix. Note: we computed transit depths directly from R_p/R_s.

Download table as:  ASCIITypeset image

We applied the same methodology on a wavelength bin by wavelength bin basis to produce the simultaneously fitted light curves for bin scheme 3 and create our final transmission spectrum, shown in Figure 8. We used the open source package ld-exosim²⁰ to determine that a square-root limb-darkening law was the most appropriate for WASP-43 and incorporated this into our GP analysis. We present the final transmission spectrum using this method in Figure 8.

3.2.4. White-light Curve and Binning Schemes

We applied the detrending methods described in Sections 3.2.1–3.2.3 to the following wavelength-binning schemes in our analysis:

• 1.  
  Bin scheme 1: a set of uniform bins centered around the air wavelength values of key spectral features (Na i-D, Hα, K i, Na i-8200) to produce a transmission spectrum focused on these features.
• 2.  
  Bin scheme 2: binning and system parameters identical to those in Murgas et al. (2014) to directly compare our transmission spectra with the ones presented in that study.
• 3.  
  Bin scheme 3: similar 230 Å binning and system parameters to Kreidberg et al. (2014) to combine our visual measurement with their NIR measurements made with HST. We used this binning scheme to perform atmospheric retrievals described in Section 5.4. We centered 230 Å bins on the vacuum wavelength locations of Na i-D, Hα, K i, and Na i-8200 using smaller bins when necessary to have at least two wavelength bins between each feature.

We applied the poly+CMC detrending method described in Section 3.2.1 to all three bin schemes and applied GP+PCA described in Section 3.2.3 to bin scheme 2 (25 nm) and bin scheme 3 because of their similar wavelength binning and wavelength coverage in the visual (25 nm versus 23 nm). This allowed us to directly compare detrending methods between our study and Murgas et al. (2014).

We discuss each binning scheme in more detail in Sections 5.2, 5.3, and 5.4. We also include the literature values for system parameters used in each bin scheme in Table 2. For bin scheme 1, where no associated literature values are being used for comparison, we adopt the most up-to-date values from Hoyer et al. (2016).

Before combining the transmission spectra from each night, we first considered the impact of stellar photospheric heterogeneity, which can have an observable effect on transmission spectra (Pont et al. 2008, 2013; Sing et al. 2011; Oshagh et al. 2014; Zhang et al. 2018), even if magnetically active regions are not occulted by the transiting exoplanet (McCullough et al. 2014; Rackham et al. 2018, 2019; Apai et al. 2018). Qualitatively, global variations in stellar activity could manifest themselves as an overall dimming or brightening of the star, which could lead to significant variations in transit depths. Changes in photometric activity roughly correlate with the covering fraction of starspots, which in turn can modulate the luminosity of the star and impact observed transit depths (Berta et al. 2011). Furthermore, those variations can be wavelength dependent, leading to slopes with spurious spectral features in the transmission spectrum.

The white-light transit depths we observed (Table 3) varied by as much as 1869 ppm between transits. To account for this offset between the data sets, we investigated the contribution due to stellar activity. Changes in photometric activity roughly correlate with the covering fraction of starspots, which in turn can modulate the luminosity of the star and impact observed transit depths (Berta et al. 2011). To assess the brightness variation of WASP-43 over the time frame of our observations, we used 2012 February 15–2018 May 21 activity data from 990 out-of-transit V-band images of WASP-43b taken by the Ohio State University's All-Sky Automated Survey for Supernovae²¹ (ASAS-SN) program (Shappee et al. 2014; Kochanek et al. 2017).

The ASAS-SN photometric activity was sampled much more coarsely than our transit observations, so we used a regression routine (Alam et al. 2018) to fit the data and estimate the amplitude of the photometric variation induced by stellar activity during each of the four transit epochs of WASP-43b. Following Alam et al. (2018), we used a negative log-likelihood kernel $(\mathrm{ln}{\boldsymbol{L}})$ for the objective function given by

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{ln}({\boldsymbol{L}}) & = & -\displaystyle \frac{n}{2}\mathrm{ln}(2\pi )-\displaystyle \frac{1}{2}\mathrm{ln}(\det {\boldsymbol{K}})\\ & & -\,\displaystyle \frac{1}{2}{\left({\boldsymbol{y}}-\mu \right)}^{\top }{{\boldsymbol{K}}}^{-1}({\boldsymbol{y}}-\mu ),\end{array}\end{eqnarray} \tag{ 6 }$$

where ${\boldsymbol{y}}$ is the data, μ is the model, n is the number of observations, and ${\boldsymbol{K}}$ is the covariance matrix. ${\boldsymbol{K}}$ describes the correlation weight between all possible pairs of photometric measurements and is populated with the GP kernel to quantify the correlation of pairs of observations. We used a gradient-based optimization routine to find the best-fit hyperparameters and used the 15.6 day stellar rotation period from Hellier et al. (2011). Figure 4 shows the GP regression model for the relevant ASAS-SN data.

Zoom In Zoom Out Reset image size

Figure 4 shows the complete ASAS-SN light curve, and Figure 5 details the ASAS-SN photometry near each of our transit epochs. Overall, the relative flux from the photometric monitoring varies by as much as 3% from the median value obtained from the GP. Observations with ASAS-SN are too coarse to effectively sample the photometric activity during times of transits. However, these data still give us a rough idea of differences in stellar flux between epochs.

Zoom In Zoom Out Reset image size

We conclude that the photometric activity data alone are not enough to constrain the contribution of unocculted heterogeneities on the surface of WASP-43b to the resulting transmission spectrum. Nonetheless, we argue that changes in disk coverage by unocculted heterogeneities likely drive the white-light light-curve depth variations that we observe between transit epochs. For this reason, we calculate and apply transit depth offset corrections as described in Section 5.1 before building the final transmission spectrum. In Section 5.4, we model the possible contribution of an unocculted heterogeneous photosphere to the resulting transmission spectrum without relying on photometric monitoring data.

5.1. Combining Nights

5.2. Bin Scheme 1: Species-dependent Binning

In this binning scheme, we set the wavelength bin sizes based on the absorption bandwidths of features of interest, in particular Na i-D, Hα, K i, and Na i-8200. We set the minimum bin size for a given feature to be equal to the full width of its observed stellar absorption line, including the contributions from the wings of the line. This gives a set of four bin widths equal to 60, 10, 60, and 40 Å for the respective species listed above. For each species, we mapped a region covering five times its bin width above, centered on the air wavelength to resolve any potential peaks, and we used larger bins to cover the rest of the spectrum. We produce the combined spectra in Figure 6 following this procedure for the (poly+CMC) detrending scheme.

Zoom In Zoom Out Reset image size

We observed an apparent peak near the Na i-8200 line, but the fact that the only point far from the baseline is also far from the air wavelength for this species indicates that this peak is most likely due to residuals from water tellurics. We also observed a potential K i peak 1σ above the median in bin two, 10 Å immediately redward of its air wavelength location, but believe that this is due to residual telluric absorption as well. Furthermore, we do not detect this peak at all in the GP+PCA detrended data. We also do not detect an absorption peak near Na i-D or Hα. To set this bin scheme apart from the other two bin schemes, we used more up-to-date system parameters from Hoyer et al. (2016).

5.3. Bin Scheme 2: Comparison to Ground-based Study

We applied the same detrending and combining methods using the system parameters reported by Murgas et al. (2014) to compare our transmission spectra with those from their similar study of WASP-43b. We adopted the following four wavelength bin schemes identified in their work: (i) 100 Å bins ranging from 5445 to 8845 Å, (ii) 250 Å bins ranging from 5300 to 9050 Å, (iii) 750 Å bins ranging from 5300 to 9050 Å, and (iv) 180 Å bins centered near the K i 7665 and 7699 Å doublet to produce Figure 7.

Zoom In Zoom Out Reset image size

Qualitatively, the shapes of the combined spectra in each binning scheme tend to follow the same slight upward curving slope near 7000 Å seen in Murgas et al. (2014). Unlike their finding, we do not observe an excess near 5892.9 Å that would indicate the presence of Na i in the atmosphere of WASP-43b.

5.4. Bin Scheme 3: Combining with NIR Study

We used the same bin width and system parameters (Table 2) from Kreidberg et al. (2014) to combine their NIR transmission spectrum with our visual spectrum. We found that although the polynomial+CMC detrending method tended to produce smaller error bars on average in the transmission spectrum (180 ppm versus 210 ppm), GP+PCA does a better job overall at fitting systematics in the light curves, resulting in an average standard deviation across the transit epochs of 396 ppm versus 274 ppm in the WLC model residuals. The average GP+PCA WLC residuals were also larger than the average white noise of 146 ppm estimated in the poly+CMC detrended WLCs, so the GP+PCA procedure was not overfitting the data. For these reasons, we adopted the GP+PCA detrended transmission spectrum, shown as black circles in Figure 8, as our final spectrum (reported in Table 7), which we combined with the NIR transmission spectrum from Kreidberg et al. (2014) for the atmospheric retrieval analysis.

Zoom In Zoom Out Reset image size

With the final combined visual+NIR transmission spectrum from Section 5.4, we searched for signals in the atmosphere of WASP-43b. We used a Bayesian atmospheric retrieval code based on the same nested sampling Bayesian inference software we used to perform the GP+PCA detrending, PyMultiNest. The details of the retrieval code are given in our previous study of WASP-19b (Espinoza et al. 2019), and we only briefly summarize the methodology here.

Following the semianalytical formalism from Bétrémieux & Swain (2017) and Heng & Kitzmann (2017), we assume an isothermal and isobaric atmosphere, with an optically thick base region with radius (R_p/R_s)₀ and reference pressure P₀, which we interpret as the cloud-top pressure. Above this region is an optically thin planetary atmosphere with average temperature T that can have either (i) a set of atomic and molecular species and/or (ii) a scattering haze defined by ${\sigma }_{\mathrm{haze}}{(\lambda )=a{\sigma }_{0}(\lambda /{\lambda }_{0})}^{{\gamma }_{\mathrm{haze}}}$ (MacDonald & Madhusudhan 2017), where σ₀ = 5.31 × 10⁻²⁷ cm² is the Rayleigh scattering cross section of H₂ at the reference wavelength λ₀ = 350 nm, and a and γ_haze are free parameters. We constrain γ_haze to be between 0 (uniform opacity) and −4 (Rayleigh scattering) to allow for a better constraint on a. Transmission spectra from separate studies can be combined by retrieving for an offset between the different data sets. A detailed overview of the retrieval framework is given in Appendix D of Espinoza et al. (2019).

Additionally, we explored the impact of a heterogeneous stellar photosphere on the observed transmission (Pinhas et al. 2018) by following the formalism described by Rackham et al. (2018, 2019). To summarize directly from the schematic in Figure 1 of Rackham et al. (2018): during a transit, exoplanet atmospheres are illuminated by the portion of a stellar photosphere immediately behind the exoplanet (the transit chord). Changes in transit depth must be measured relative to the spectrum of this light source. However, the light source is generally assumed to be the disk-integrated spectrum of the star. Any differences between the assumed and actual light sources will lead to apparent variations in transit depth.

In this framework, stellar contamination of the transmission spectrum from unocculted starspots and faculae is considered by placing constraints on the allowed spot- and faculae-covering fractions using a set of rotating photosphere models and then translating the covering fractions of potential stellar contamination in the transmission spectrum. We incorporate this transit light-source effect into our retrieval framework with a three-parameter model for the stellar photosphere to fit this simultaneously with the planet's atmosphere. The three parameters are T_chord, the effective temperature of the transit chord; T_het, the mean effective temperature of the heterogeneous features not occulted by the transit chord; and f_het, the fraction of the projected stellar disk covered by these heterogeneous features. The impact of these heterogeneities on the transmission spectrum is expressed by the wavelength-dependent corrective factor _λ on the transit depth, where

$$\begin{eqnarray}&&{\epsilon }_{\lambda }\equiv {\left[1-{f}_{\mathrm{het}}\left(1-\displaystyle \frac{{S}_{\lambda ,\mathrm{het}}}{{S}_{\lambda ,\mathrm{chord}}}\right)\right]}^{-1}={\left(\displaystyle \frac{{R}_{{\rm{p}},\lambda }}{{R}_{{\rm{p}},\lambda ,0}}\right)}^{2}.\end{eqnarray} \tag{ 7 }$$

Here, R_p,λ and R_p,λ,0 are the apparent and actual planetary radius measured at wavelength λ, respectively; S_λ,het is the spectrum of the unocculted photosphere determined by T_het; and S_λ,chord is the spectrum of the portion of the photosphere inside of the transit chord. Following previous studies (e.g., McCullough et al. 2014; Rackham et al. 2017), we use PHOENIX stellar spectra (Husser et al. 2013) to model the emergent spectra of the photospheric components.

As Rackham et al. (2018) note, this formalism assumes that the transit chord can be described by a single emergent spectrum. Although this is not guaranteed for any one of our transits, they note that this formalism also holds for transits in which an occulted spot- or faculae-crossing event is present in the transit signal above the observational uncertainty and taken into account in the transit modeling. We also explored more complex models including multiple spot- and faculae-covering fractions, but we found that the data did not warrant the additional complexity of such a model.

Our combined retrieval approach uses the posterior Bayesian evidence $Z\equiv {\mathbb{P}}(D| H)$ computed by PyMultiNest, which is the probability of the data D given the hypothesis H to perform model comparisons. This property of the nested sampling algorithm allows us to study how complex our models have to be to explain the observed distortions to the light curve (such as number of spots) via the posterior odds, ${\mathbb{P}}({H}_{n}| D)/{\mathbb{P}}({H}_{k}| D)$, where the joint probability ${\mathbb{P}}({H}_{n}| D)=P(D| {H}_{n}){\mathbb{P}}({H}_{n})$, with ${\mathbb{P}}({H}_{n})$ being the prior probability on the hypothesis H_n. If we approximate model n and model k as having the same prior distribution on their respective hypothesis H_n, H_k, then the posterior odds simplify to just the ratio of the evidences,

$$\begin{eqnarray}&&\displaystyle \frac{{\mathbb{P}}({H}_{n}| D)}{{\mathbb{P}}({H}_{k}| D)}=\displaystyle \frac{{\mathbb{P}}(D| {H}_{n})}{{\mathbb{P}}(D| {H}_{k})}\equiv \displaystyle \frac{{Z}_{n}}{{Z}_{k}}.\end{eqnarray} \tag{ 8 }$$

In log space, this is the difference of the logs of each evidence and is denoted as ${\rm{\Delta }}\mathrm{ln}Z\equiv \mathrm{ln}{Z}_{n}-\mathrm{ln}{Z}_{k}$. In the context of this study, Δln Z is a measure of how statistically different a given model (Z_n) is from a flat atmospheric model (Z_k ≡ Z_flat). Trotta (2008) and Benneke & Seager (2013) review how these log odds translate to frequentist significance hypothesis testing. We note from that work that absolute log odds below 1 are usually considered inconclusive, near 2.5 can be interpreted as moderate evidence, and higher than 5 can be interpreted as highly significant. It is important to caution though that frequentist hypothesis testing has only one null hypothesis, whereas proper Bayesian model comparison considers a range of possible hypotheses, which limit the comparison with frequentist methods.

6.1. Retrieval Results

We fit the combined GP+PCA ACCESS detrended spectrum and the HST/WFC3 NIR spectrum from Kreidberg et al. (2014) using the retrieval code described in the previous section. We fit for an offset between the two data sets, and a range of models. Those include combinations of clear and cloudy/hazy atmospheres with Na, K, and H₂O and contamination of the planet's transmission spectrum by stellar surface heterogeneity. Table 4 shows the prior distributions of the parameters used in the models. We used the prescription provided by Benneke & Seager (2013, Table 2) to interpret our relative log evidences (Table 5). Based on their values, we find that they all register as a strong detection, with no one model being statistically more likely than another, although the model with the largest log evidence relative to a flat atmosphere (Δln Z = 8.26) is the one including stellar heterogeneity, combined with a clear atmosphere with H₂O (but no Na or K). We adopt that model as the one that best fits the data and show it in Figure 9. From this model, we estimate an average spot contrast of 132 K ± 132 K and covering fraction of ${0.27}_{-0.16}^{+0.42}$. The corner plot for that model solution is shown in the Appendix. The parameters of that model and their uncertainties are summarized in Table 6.

Zoom In Zoom Out Reset image size

Table 4.  Priors Used in Retrieval Models

Model Component	Parameter	Units	Description	Prior Distribution
Offset	offset	ppm	Offset between Magellan/IMACS and HST/WFC3 data	Normal (Mean Depth, 1000 ppm)
Base	(Rp/Rs)0	⋯	Radius corresponding to the top of the cloud layer or τ ≫ 1	Uniform (0.8, 1.2)
	P0	bar	Reference pressure at (Rp/Rs)0	Log uniform (10−4, 1)
	T	K	Average temperature planet atmosphere	Uniform (0, 1500)
Atomic features	X	⋯	Mixing ratio of species X	Log uniform (10−14, 1)
Haze	a	⋯	Amplitude of the haze cross-section power law	Log uniform (10−10, 1020)
	γ	⋯	Index of the haze cross-section power law	Uniform (−10, 0)
Stellar photosphere	Tocc	K	Average temperature of the transit chord	Uniform (4000, 5000)
	Thet	K	Average temperature of the heterogeneous surface features	Uniform (4000, 5000)
	Fhet	⋯	Fraction of the unocculted photosphere covered by spots	Uniform (0, 1)

Download table as:  ASCIITypeset image

Table 5.  Δln Z for Various Models Relative to a Flat Spectrum for the ACCESS GP+PCA Detrended Data, Combined with the Kreidberg et al. (2014) Data

Model	H2O	Na+H2O	Na+K+H2O
clear	4.98	6.42	6.08
haze	5.43	5.76	5.60
spot	8.26	7.90	7.74
spot+haze	6.92	6.70	6.59

Note. The largest value of 8.26 corresponds to the best retrieved transmission spectrum model shown in Figure 9.

Download table as:  ASCIITypeset image

Table 6.  Retrieved Parameters for the Best-fit Retrieval Model Shown in Figure 9

Parameter	Description	Value
Fhet	Spot-covering fraction	${0.27}_{-0.16}^{+0.42}$
T	Planet temperature (K)	${352.91}_{-125.08}^{+206.14}$
Thet	Spot temperature (K)	${4169.40}_{-92.33}^{+86.10}$
Tocc	Occulted temperature (K)	${4300.91}_{-77.16}^{+94.04}$
logH2O	Water volume mixing ratio	$-{2.78}_{-1.47}^{+1.38}$
logP0	Reference pressure (bar)	$-{1.67}_{-1.16}^{+1.08}$
offsetMagellan/IMACS	Offset between data sets (ppm)	${25749.40}_{-105.70}^{+101.54}$
f	Planet radius normalization	${0.99}_{-0.00}^{+0.00}$

Download table as:  ASCIITypeset image

We have collected, extracted, and combined transmission spectra of WASP-43b from Magellan/IMACS over four transit epochs spanning the years 2015 to 2018. We combined this with IR data from HST/WFC3 to create a transmission spectrum with a total wavelength coverage of 5318–16420 Å. We analyzed the combined spectrum in a dynamic nested sampling framework with NIR data from Kreidberg et al. (2014), extending up to 16420 Å to search for the presence of different species. Assuming a water volume mixing ratio of 6.1 × 10⁻⁴ for a planetary atmosphere with solar abundances (Kreidberg et al. 2014), our retrieval yields a log H₂O volume mixing ratio of $-{2.78}_{-1.47}^{+1.38}.$ (${2.72}_{-0.09}^{+65.26}$ solar). Our retrieved water abundance is consistent with the 1σ range found by the joint transmission and emission spectrum analysis in Kreidberg et al. (2014, 0.4–3.5 solar) and phase-curve analysis in Stevenson et al. (2014). Our retrieved planetary temperature is also consistent with the ranges predicted from the more recent 2.5D theoretical phase-curve retrievals of the terminator region in Irwin et al. (2019, Figure 9). In the visual spectrum, we do not observe a statistically significant excess of K i given the data, as reported in Murgas et al. (2014). We also do not observe the presence of Na i or Hα in our combined spectra.

Our analysis also investigates the contribution of stellar heterogeneities to observed transmission spectra, given that WASP-43 is an active star. The best-fitting model to our transmission spectrum calls for the presence of spots in the surface of the star being a more favored model than a model with atmospheric hazes. However, we do note that the impact on the spectrum from both hazes and the contamination from surface stellar heterogeneities can be degenerate, and given the quality of the current data, it is not possible to fully distinguish between the two.

As the sample of available transmission spectra of different planets increases, it will become useful to compare their measured spectra and other system parameters (e.g., stellar irradiation levels, spectral type, metallicity, planet density) to look for possible correlations. Uniform data sets for a wide range of planets, such as the ones ACCESS is building, will be crucial in this analysis.

The results reported herein benefited from support, collaborations, and information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. This paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile. We thank the staff at the Magellan Telescopes and Las Campanas Observatory for their ongoing input and support to make the ACCESS observations presented in this work possible. A.J. acknowledges support from FONDECYT project 1171208 and by the Ministry for the Economy, Development, and Tourism's Programa Iniciativa Científica Milenio through grant IC 120009, awarded to the Millennium Institute of Astrophysics helped improve the quality of this work.

Facilities: Magellan: Baade - , Smithsonian Institution High Performance Cluster (SI/HPC). -

Software: Astropy (Astropy Collaboration et al. 2013), corner (Foreman-Mackey 2016), Matplotlib (Hunter 2007), NumPy (Oliphant 2006; van der Walt et al. 2011), PyMC (Salvatier et al. 2016), Multinest (Feroz et al. 2009), PyMultiNest (Buchner et al. 2014), SciPy (Virtanen et al. 2019), SPOTROD (Béky et al. 2014), batman (Kreidberg 2015), george (Ambikasaran et al. 2014).

Below are the additional figures referenced in the text.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size