ACCESS: Ground-based Optical Transmission Spectroscopy of the Hot Jupiter WASP-4b

On the first three nights, we used a setup consisting of a 400–1000 nm spectroscopic filter, a 300 lines/mm grism with a blaze angle of 175, and 10'' wide by 20'' long spectral slits for the target and reference stars. Twelve reference stars were observed, although only one was used in the final data reduction for reasons discussed in Section 3.3. Most of the stellar spectra were dispersed across two chips. While the lunar sky background was minimal on the first and third nights, it was substantial during the second, and we take this into account in our data-reduction pipeline. Finally, the observations on 2013 September 24 (Transit 1) did not commence until shortly after ingress, so we were only able to observe a partial transit.

Our instrument setup for the fourth observation was modified from the previous three. We used a 570–980 nm order-blocking filter for the purpose of eliminating higher-order interference toward red wavelengths, 10'' wide by 10'' long slits, and a 150 lines/mm grism with a blaze angle of 188. This setup allowed for more tightly dispersed spectra that fall on a single chip, thereby reducing detector-to-detector variations in the spectra and avoiding chip gaps. The moon was in full phase and separated from the target by 40°, contributing significantly to the ambient sky background, which we again account for in our data-reduction pipeline.

In some of our IMACS data sets for this and other ACCESS targets, we have noticed that spatially uniform sources of light inherit a nonuniform profile in the spatial direction once the mask and dispersive element are in place. This effect also widens the spectra of point sources. We have attributed this to internal scattering within the instrument and have found that it is more common in newer images: of the four data sets for WASP-4, only Transits 3 and 4 are affected. We measure the extent of the effect using our flat-field images, finding that the profile of spatially uniform light peaks at the slit center and appears ∼10% fainter at the slit edges. Unless we account for this, we will underestimate the sky background when we extract the stellar spectra.

For the two affected data sets, we use quartz lamp exposures taken at the beginning of the night to model the scattered light profile. The peaked profiles are not consistent with common symmetrical functions (e.g., Gaussian), nor are they well described by a classical high-order polynomial, because high-order polynomials commonly fail at the edges of their fitted intervals (an effect known as Runge's phenomenon). Instead, we use Chebyshev polynomials, which are not as susceptible to this effect, and we find that sixth-order polynomials are sufficient to match the data.

The polynomials are fitted independently for every row along the dispersion direction of every slit in the flat-field images. Then, in each science image we fit the amplitudes of the polynomials using the background level in the outermost 14 px in each row along the dispersion direction (8 px for Transit 4 due to shorter slits). Finally, we subtract the fitted polynomial from each row.

Applying the flat-field correction for Transits 1, 3, or 4 does not significantly change our binned light curves or the transit spectra that we derive in Section 4. However, applying the correction for Transit 2 introduces ∼2%-level, nonlinear time-dependent trends into the out-of-transit baseline flux of the binned light curves blueward of 530 nm. These trends remain even when normalizing by the reference star, as discussed in the following. Therefore, for consistency, we choose not to apply the flat-field correction to any of our data sets.

The shape of the target's light curve is complicated by instrumental and atmospheric effects, such as changes in air mass or transparency. To calibrate out these effects, we simultaneously observed 12 reference stars of comparable optical apparent magnitudes and color ratios using multislit masks. Of these, two spectra reached the saturation limit of the detector and were not usable.

We use our highest quality data set (Transit 3) to determine which of the 10 remaining reference stars to use in our light curve analysis. Our primary consideration is the shape of the out-of-transit baseline for the target light curve when normalized by each star. Eight of the stars leave residual long-term trends at the ∼1% level. Two stars leave trends at the ∼0.1% level, and they happen to be the only two stars occupying the same pair of detector chips as the target. The point-to-point scatters of the baseline points for the two resulting light curves are 0.4 and 0.6 mmag, which are both smaller than the scatters for the eight light curves with larger long-term trends; we therefore further restrict our set of potential reference stars to these two. We perform a similar analysis of the two remaining stars using our three other data sets and discover that one of the two stars reliably produces a flat baseline, while the other star introduces trends at the ∼0.5%–1% level. Therefore, we discard this star as well.

We are left with a single reference star, 2MASS J23341836-4204509, which reliably produces a flat out-of-transit baseline with a low point-to-point scatter in all of our data sets. We note that this reference star was also used in a similar analysis by Huitson et al. (2017). The spectral type of the reference star has not been previously published, but using parallax measurements from Gaia DR2 (Gaia Collaboration 2018), we calculate that it is intrinsically brighter than the target by ∼80% from 400 to 1000 nm. While this wavelength range does not capture all of the emitted light, it does cover the emission peak for late F and all G stars. Using a simple ${L}_{* }\propto {M}_{* }^{3.5}$ scaling relation and assuming that our wavelength range covers most of the emitted light suggest that the reference star is ∼20% more massive than the G7V 0.85 M_⊙ target (Gillon et al. 2009), consistent with an early G or late F dwarf.

Magellan-Baade is equipped with an atmospheric dispersion compensator (ADC), which corrects for the effects of differential atmospheric refraction as the air mass of our target changes over the course of the night. If left uncorrected, differential refraction would affect our stellar wavelength solutions in a nonlinear manner, "stretching" or "shrinking" the spectra over the course of the night. We test the accuracy of the ADC correction by measuring the distance in calibrated wavelength space between Na and Hα in the target and reference spectra as a function of time, which we calculate to vary by no more than 1 Å over the course of the night. Furthermore, the difference between the measured distance for the target and reference stellar spectra changes by no more than 0.5 Å, and less than 0.1 Å on nights with higher quality data. Since the long-term change in this measurement is much lower than the frame-to-frame scatter (up to 5 Å), we reason that any effect on the measured transit depth will be negligible.

During Transits 1 and 2, we opted for longer, low-noise readout modes and short exposure times to avoid detector nonlinearity. As a result, we spent only ∼20% of our time gathering light. By comparison, we achieved an efficiency of ∼70% and ∼50% during Transits 3 and 4, respectively, as we chose the fastest readout mode and took longer exposures. The difference in read noise between the fastest (5.6 e-) and slowest (2.8 e-) modes is negligible compared to the Poisson noise (>100 e-), and as we show in Section 3.8, the detector is sufficiently linear for our purposes provided the pixel values remain near or below half-well.

Based on this experience, we advise that observers prioritize longer exposure times with fast readouts for transit spectroscopy with IMACS.

Transits 2 and 4 suffer from high sky background values due to the target's proximity to the full-phase moon. The effects can be measured by comparing the scatter of the residuals in Figure 2. Scaled for exposure time, the magnitudes of the scatter for Transits 1 and 3 agree to within 15%, while the scatter of Transit 2 is ∼50% larger. The scatter of Transit 4 is similarly ∼50% larger, although Transit 4 was observed in a redder filter with less sky contamination.

We consider whether slit losses may have been significant during our observations. For Transits 1–3, the slits were 10'' wide and 20'' long, while for Transit 4 the slits were 10'' wide and 10'' long. We perform a least-squares fit of a Moffat profile to the point-spread function of a bright but unsaturated star in each of the field acquisition images, and then we integrate to determine how much light would lie outside the slit if placed over the star. For every night, the fraction of light lost is more than an order of magnitude lower than the Poisson noise.

We observed Transits 3 and 4 with longer integration times to reduce the relative readout overhead and improve our signal-to-noise ratios (S/N). However, this increased the maximum pixel values of the target spectrum from ∼15,000 to ∼30,000 ADU, nearly half of the full well value (∼65,000 ADU). This raises the question of whether the detector's response could have become nonlinear (i.e., the gain was not uniform over the range of pixel values), which could bias the mean-subtracted transit spectra that we extract later in Section 4.

In principle, nonlinearity should not be an issue for measuring the relatively small changes in transit depth from bin to bin, as the gain should not be expected to change over the small range of changing values. However, the pixel values at the peak of the target spectrum varied by up to 15% from frame to frame and 40% over the course of each night, due mostly to seeing variations that changed the point-spread function. Therefore, it is worthwhile to determine the threshold beyond which nonlinearity could significantly affect our results.

To do so, we modify the raw data for Transit 3 to reflect a 0.1% or 1% linear increase in the gain over the range 0–35,000 ADU. We then reduce the modified data and extract the mean-subtracted transit spectra for both gain prescriptions, which we plot alongside the original spectrum in Figure 1. We find an average effect of 0.07σ per bin for a 0.1% increase in gain, and 0.23σ per bin for a 1% increase.

Zoom In Zoom Out Reset image size

So long as the detector is linear to 0.1% from zero to half-well, nonlinearity effects should have a negligible impact on our results. Even in the case of 1% nonlinearity, the effect is ≪1σ except in a few bins. We are therefore confident that our results are robust to realistic levels of nonlinearity, and we advise that future observers target a similar range of peak pixel values (∼35,000 ADU) so as to optimally balance observational efficiency and detector linearity.

The average values of R_p/R_s in the spectra of Transits 1–3 agree to within 1.5σ, but Transit 4 disagrees at the ∼5σ level. This mirrors the deviation in the fitted orbital parameters for Transit 4 discussed in Section 4.3, and the causes are likely the same.

Regardless of the reason for this disagreement, we are primarily interested in the change in R_p/R_s with wavelength. We subtract the weighted mean $\overline{{R}_{p}/{R}_{s}}$ from the individual spectra for Transits 2–4, and then calculate the average values of ΔR_p/R_s in each bin, weighted by the uncertainties for each night. The combined transmission spectrum is presented in Figure 4.

Zoom In Zoom Out Reset image size

Because of the large uncertainties and incomplete baseline coverage for Transit 1, we do not include its spectrum in the combined result; see Section 4.6. We also correct the spectrum of Transit 3 for known stellar contamination before incorporating it into the combined spectrum, as detailed in Section 4.7.

Two of the fitted bins have been excluded from our primary figures and results: the 450–470 nm bin of Transit 2, and the 875–900 nm bin of Transit 4. The values for both bins deviate by >3σ from the average value of their respective spectra, and both bins are at the low S/N ends of the stellar spectrum. For reference, these bins are included in the figures and table in Appendix A and B.

The 450–470 nm bin of Transit 2 lies ∼3σ below the average value of the spectrum; we can think of no astrophysical explanation for this effect. Furthermore, the likelihood of a statistical 3σ outlier across all of our 82 bins is ∼0.3%. Most likely, our polynomial model is inadequate to describe the low S/N systematic trends in this bin.

The 875–900 nm bin of Transit 4 lies ∼4σ above the weighted mean. While this wavelength range does correspond with water absorption, such an offset is not observed in the other spectra. Furthermore, the value of this bin depends strongly on whether a linear or quadratic polynomial is used to model the statistics, as discussed in Section 4.4. Again, it seems most likely that our systematics model is inadequate for this low-S/N bin.

We also use our MCMC code to simultaneously fit for the transit depth, the inclination (i), and the semimajor axis (a/R_s), while keeping the eccentricity and period fixed. The results of these fits are listed in Table 3. The parameters from Transits 1–3 are generally consistent with those measured by Hoyer et al. (2013), although the values from the partial Transit 1 are poorly constrained.

Table 3.  Fitted Orbital Parameters When Modeling the Planet Radius, Inclination, and Semimajor Axis Jointly

	Rp/Rs	a/Rs	i (deg)
Hoyer+13	${0.15445}_{-0.00025}^{+0.00025}$	${5.463}_{-0.020}^{+0.025}$	${88.52}_{-0.26}^{+0.39}$
Transit 1	${0.15386}_{-0.00174}^{+0.00135}$	${5.386}_{-0.157}^{+0.088}$	${89.33}_{-0.73}^{+0.47}$
Transit 2	${0.15617}_{-0.00167}^{+0.00164}$	${5.439}_{-0.092}^{+0.070}$	${87.90}_{-0.85}^{+1.24}$
Transit 3	${0.15388}_{-0.00047}^{+0.00045}$	${5.473}_{-0.038}^{+0.021}$	${88.94}_{-0.52}^{+0.53}$
Transit 4	${0.15654}_{-0.00065}^{+0.00067}$	${5.520}_{-0.030}^{+0.026}$	${89.46}_{-0.60}^{+0.37}$

Note. For comparison we include the values derived by Hoyer et al. (2013).

Download table as:  ASCIITypeset image

The parameters from Transit 4 deviate 2–3σ from the previously measured values. It is not clear why this is the case. Stellar variability may be a partial contributor; as discussed in Section 6.2, WASP-4 is variable at the 6 mmag level, and this can lead to changes in the apparent transit depth. The deviation may also be related to losses that are due to the shorter slits employed during Transit 4, as discussed in Section 3.7.

Because of this large deviation, we also try fitting the spectrum of Transit 4 by fixing the inclination and semimajor axis to the fitted values, rather than those in the literature. The result is only a ≪1σ change in the binned values for ΔR_p/R_s, so we believe the spectrum for this night is robust despite the disagreement in the orbital parameters. For consistency with the other data sets, we keep the parameters fixed to the literature values in the remaining sections.

We find that the trends that remain in the target light curves after dividing by the reference starlight curves correlate with air mass, which tends to reach its minimum near the center of the transit. Nevertheless, no physically motivated (e.g., exponential) function of air mass is able to describe this trend better than a simple second-order polynomial in time, so we opt to use time-dependent polynomials to characterize our long-term systematics.

The long-term trends in many of our binned light curves cannot be adequately described by a linear polynomial in time, so we use a quadratic polynomial instead. In principle, any degeneracy between the transit depth and the polynomial value should be incorporated into our MCMC uncertainties. To measure the impact of the model choice on our results, we apply both linear and quadratic systematics models and compare the values of ΔR_p/R_s for each.

For Transits 1–3, the average absolute change in ΔR_p/R_s between the linear and quadratic models is approximately 0.5σ per bin and no larger than 1.6σ in any bin. For Transit 4, the average change is approximately 1.5σ per bin, but the reddest bin (which has the lowest S/N) changes at the 3σ level. Because of this large inconsistency between systematics models, we exclude the reddest bin from our results.

Some authors use the residuals of their white light curve models to calibrate out common (i.e., wavelength-independent) systematics from each bin (e.g., Sing et al. 2013; Nikolov et al. 2015; Huitson et al. 2017). In our case, we find that such a correction is unwarranted. By eye, it is not apparent that our binned light curves share any common systematics (see Appendix A). Furthermore, when we apply the common mode correction, we find that the impact on the resulting transit spectra is negligible compared to the uncertainties.

As discussed in Section 4, we fit a parameter σ_r to characterize the level of correlated ("red") noise in our data. While the ratios in Table 1 suggest a large amount of systematic noise, this does not necessarily imply correlated noise; the systematic noise may still be captured by the white noise parameter.

To assess the level of red noise in our data, we construct Lomb–Scargle periodograms of the residuals to each binned light curve model, and then we use bootstrap resampling to calculate the false alarm probability (FAP; e.g., VanderPlas 2018) for signals with periods between the exposure times and the duration of our observations. Figure 5 shows an example periodogram for one of the bins of Transit 3. We find that only six out of 82 binned light curves possess signals with FAP < 10%, and only two have FAP < 1%. This suggests that there are few strong periodic signals in our light curves. Nevertheless, we conservatively choose to apply the red noise parameterization, which results in a ∼15% increase in our uncertainties.

Zoom In Zoom Out Reset image size

The light curves of Transits 1–3 all suffer from a partial lack of phase coverage, due either to poor weather, observing window constraints, or spot-crossing events that must be masked in the model fit. In principle, our MCMC approach to fitting the model should adequately incorporate this lack of information into the uncertainties on R_p/R_s. However, it is worthwhile to estimate the extent to which this missing information inflates our uncertainties.

The full light curve of Transit 4 allows us to investigate the effects of poor phase coverage in the previous three transits. Specifically, we remove data points from the binned light curves of Transit 4 to mimic the phase coverage of Transits 1–3, and we then fit for the transit depth and calculate the resulting increase in the uncertainties on the binned values of R_p/R_s. Our results are summarized as follows:

• 1.  
  Mock Transit 1: We remove the preingress baseline and most of the first half of the transit. This results in a 150% increase in uncertainties, as well as a 4.5σ decrease in the wavelength-averaged radius. The ΔR_p/R_s values change at the >1σ level, and the change appears to be wavelength dependent, introducing a slope on the order of ∼1σ from blue to red.
• 2.  
  Mock Transit 2: We remove the points near midtransit that correspond to a spot-crossing event. This results in a 50% increase in uncertainties and a 0.5σ decrease in the average radius. The ΔR_p/R_s values change at the ≲0.2σ level.
• 3.  
  Mock Transit 3: First, we remove most of the postegress baseline, which results in a 50% increase in uncertainties with no considerable change in the average radius and ≲0.5σ changes in ΔR_p/R_s. Next, we exclude only those points that correspond to the spot-crossing event in Transit 3, finding effects of similar magnitude. Removing both portions of the light curve, however, does not increase the uncertainties any further.

The results for Transit 1 lead to concerns over whether the lack of a preingress baseline could introduce a bias into the shape of our transit spectrum. Furthermore, the large uncertainties mean that this transit contributes little to the combined spectrum. We therefore opt to exclude the spectrum of Transit 1 from the combined result.

Meanwhile, the results for Transits 2 and 3 suggest that little bias has been introduced into the transit spectrum, due to the lack of phase coverage in each, while the uncertainties on ΔR_p/R_s should be larger by ∼50% than they would be given full phase coverage.

The light curves of Transits 2 and 3 feature prominent spot occultation features that we exclude from the light curve model fit. While doing so allows us to fit the apparent transit depth without simultaneously modeling the spot, it does not eliminate the impact of the occulted spot on the resulting transit spectrum. Since the spot is cooler than the rest of the photosphere, it breaks the transit model's assumption that the photosphere is azimuthally symmetric and amplifies the fitted transit depth in a wavelength-dependent manner (e.g., Deming et al. 2013). In previous works we have referred to this as the transit light source effect (TLSE); here we derive an approximate correction (ignoring limb-darkening) for the TLSE due to a single spot.

Assuming the photosphere to have specific flux F_λ,phot and radius R_s, the spot to have specific flux F_λ,spot and radius R_spot, and the planet to have apparent radius R_p and true radius R_p,0, we find the observed relative transit depth will be

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{R}_{p}}{{R}_{s}}\right)}^{2}=\displaystyle \frac{{F}_{\lambda ,\mathrm{phot}}{R}_{p,0}^{2}}{{F}_{\lambda ,\mathrm{phot}}({R}_{s}^{2}-{R}_{\mathrm{spot}}^{2})+{F}_{\lambda ,\mathrm{spot}}{R}_{\mathrm{spot}}^{2}}.\end{eqnarray} \tag{ 1 }$$

By rearranging this equation, we derive a corrective factor for the transit depth:¹⁰

$$\begin{eqnarray}&&{\epsilon }_{\lambda }\equiv {\left(\displaystyle \frac{{R}_{p}}{{R}_{p,0}}\right)}^{2}={\left[1-\left(1-\displaystyle \frac{{F}_{\lambda ,\mathrm{spot}}}{{F}_{\lambda ,\mathrm{phot}}}\right){\left(\displaystyle \frac{{R}_{\mathrm{spot}}}{{R}_{s}}\right)}^{2}\right]}^{-1}.\end{eqnarray} \tag{ 2 }$$

Later, in Section 5, we find the feature in Transit 3 to be best described by a model of a single spot with radius ${R}_{\mathrm{spot}}/{R}_{s}\,={0.27}_{-0.02}^{+0.03}$ and a temperature contrast of T_s − T_spot = 100 ± 5 K. While this feature could in principle be due to a mixture of occulted spots and faculae, our data are insufficient to constrain more complex models, so we use the parameters of the single-spot model for the purposes of correcting for the TLSE in this spectrum. We interpolate a PHOENIX (Husser et al. 2013) stellar photosphere model grid onto the values in Table 2 to compute F_λ,phot and F_λ,spot at their respective temperatures, and then we calculate _λ from Equation (2).

In Figure 6 we plot ^1/2, which is the corrective factor for the planet radius, as well as for the original and corrected spectra for Transit 3. The magnitude of the correction from 450 to 900 nm is approximately equal to the binned uncertainty. The corrective factor is listed in Appendix B, binned to match our final spectrum. This corrective factor is not applied to the spectrum of Transit 3 in Figure 3, but is applied before creating the combined spectrum in Figure 4.

Zoom In Zoom Out Reset image size

We do not derive a similar correction for Transit 2 because it is unclear (see Section 5) whether this feature is best modeled by a small or large spot, and the slope introduced by the large-spot model is six times steeper across the wavelength range than that of the small spot. Furthermore, the large-spot model would only introduce a ∼0.3σ increase in the measured radius from the reddest to bluest bins.

Finally, it is worth noting that this effect is just as prominent for spots outside the transit chord. Equation (2) can be further generalized by replacing ${({R}_{\mathrm{spot}}/{R}_{s})}^{2}$ with F_het, the areal spot-covering fraction of the photosphere, and using an average spot temperature to compute F_λ,spot.

We use SPOTROD (Béky et al. 2014) to produce light curves of transits with one or two spot-crossing events. In addition to the usual transit parameters, SPOTROD models the position and radius (in units of stellar radii) as well as the average spot-to-stellar flux contrast for each spot. We employ PyMultiNest (Buchner et al. 2014) to sample the parameter space and to calculate the Bayesian log-evidences $\mathrm{ln}(Z)$ through which we can compare the different models, as discussed in Section 6.1. This MultiNest implementation¹¹ of SPOTROD has also been used to study spots observed during transits of WASP-19b (Espinoza et al. 2019).

We use a second-order polynomial to detrend the white light curve from each night, and we then fit the combined spot and star model. We limit the uniform prior on spot size to R_spot/R_s < 0.08 or 0 < R_spot/R_s < 1 to probe for small and large spots. Using the effective temperature in Table 2 as a prior, the code then constructs and fits models for the flux contrast as a function of wavelength.

The optimal model parameters for each spot feature are presented in Table 4. The precision of the Transit 2 data only permits us to fit the spot feature and spectrum in a single, white-light bin. For this night, we find that our data are consistent with one- and two-spot solutions, as well as large and small spots. We opt for the simpler single-spot models and present the fit for both a single large and a single cool spot. The data for Transit 3 favor a single large spot versus multiple or smaller spots (${\rm{\Delta }}\mathrm{ln}(Z)\gt 2$). For this transit, we are able to fit the spot spectrum to the same bins as our transmission spectrum, as shown in Figure 7.

Zoom In Zoom Out Reset image size

Table 4.  Best-fit Parameters for the Spot-crossing Events in the Light Curves of Transits 2 and 3

	Rspot/Rs	fspot/fs	Ts − Tspot (K)	Tspot (K)
Transit 2a	${0.05}_{-0.01}^{+0.01}$	${0.49}_{-0.28}^{+0.20}$	880 ± 360	4670 ± 360
Transit 2b	${0.15}_{-0.07}^{+0.11}$	${0.79}_{-0.18}^{+0.11}$	410 ± 280	5140 ± 280
Transit 3	${0.27}_{-0.02}^{+0.03}$	${0.91}_{-0.01}^{+0.01}$	100 ± 5	5442 ± 50

Note. The data for Transit 2 are consistent with both a small (2a) and a large spot (2b). Median and 1σ confidence intervals are reported.

Download table as:  ASCIITypeset image

The large-spot model for Transit 3 seems inconsistent with smaller individual spots observed on the Sun; nevertheless, such wide features have been observed in transits before (e.g., Espinoza et al. 2019). While our evidence does not favor a two-spot model, it is possible that the large, low-contrast spot is in fact a more complex active region consisting of multiple spots or faculae, as has been observed on the Sun. Regardless of the actual structure of the feature, its average contrast is well described by a single spot with T_s − T_spot = 100 K, and we choose this model to calculate the correction in Section 4.7.

Later, in Section 6, we assume a spot temperature in order to model the effects of unocculted spots on the transit spectrum. The temperature we assume reflects small, cool spots resembling those on the Sun, but the results for Transit 3 in Table 4 appear to suggest the existence of large, warm features on the photosphere of WASP-4. In short, we do not believe that the spots discovered in our transit light curves are necessarily representative of the spots present elsewhere in the photosphere or below our detection threshold. We offer these arguments as justification:

• 1.  
  Sanchis-Ojeda et al. (2011) find evidence for a persistent small spot in multiple transit observations of WASP-4b. In this case, they find a spot with radius R_spot/R_s ≈ 0.05 and temperature contrast T_s − T_spot ≈ 600 K, suggesting that more Sun-like spots do exist on WASP-4.
• 2.  
  Cool spots have previously been discovered through Doppler imaging of active G and K dwarfs (e.g., O'Neal et al. 1998, 2004; Strassmeier 2009). However, this method of spot detection is subject to biases regarding spot size and temperature.
• 3.  
  The feature occulted during Transit 3 may indeed have a more complex structure consisting of cool spots and hot faculae that the precision of our data does not allow us to resolve.
• 4.  
  While the precise radius and temperature of the spot observed during Transit 2 are largely unconstrained, we can say with at least 1σ confidence that the spot is smaller and cooler than the feature observed during Transit 3. Therefore, even our data suggest that small spots may be present in the photosphere of WASP-4.

It is nevertheless possible that the distribution of spot sizes on WASP-4 favors larger and hot spots, in which case the spectral contamination signal would be less severe. However, our use of the contamination model is not meant to accurately account for the actual surface heterogeneity, but rather to demonstrate that this aspect should not be ignored in determining the planet's transmission spectrum.

To compare between our models, we compute their relative Bayes factors. A thorough introduction to the use of Bayes factors in astronomy is given by Trotta (2008), and their use in transit spectroscopy is well established (e.g., Benneke & Seager 2013). Here, we provide a brief overview.

Given a model M with parameters θ, the optimal values of θ to describe data values x are those that maximize Bayes' equation:

$$\begin{eqnarray*}&&P(\theta | x,M)=\displaystyle \frac{P(x| \theta ,M)P(\theta | M)}{P(x| M)},\end{eqnarray*}$$

where $P(x| \theta ,M)$ is the likelihood function and $P(\theta | M)$ is the prior distribution of θ. Most Bayesian model fitting algorithms seek to maximize the likelihood function, but the maximum likelihoods of two different models cannot be directly compared.

To compare two models, we first compute their Bayesian evidences, defined as

$$\begin{eqnarray*}&&{Z}_{M}\equiv \int P(x| \theta ,M)P(\theta | M)d\theta .\end{eqnarray*}$$

The evidence for a model is high if the data are well described by a large region of the prior parameter space. The Bayes factor between two models is defined as the ratio of their evidences:

$$\begin{eqnarray*}&&\mathrm{ln}(B)\equiv \mathrm{ln}({Z}_{2})-\mathrm{ln}({Z}_{1}).\end{eqnarray*}$$

In general, model 2 is weakly favored over model 1 if 0 < ln(B) < 2, and strongly favored if $\mathrm{ln}(B)\gt 5$. If $\mathrm{ln}(B)\lt 0$, then model 2 is disfavored to model 1. Complex models will be favored over simpler models only if the evidence justifies the additional parameters. However, a preference toward simpler models does not necessarily mean that the complex models are incorrect, but rather that they are not justified by the data at hand.

Our retrieval algorithm computes the Bayesian evidence parameter for each of the fitted models. In the following sections, we cite the Bayes factor $\mathrm{ln}(B)$ for each model with respect to the model for a uniform-opacity atmosphere with no contamination from the photosphere (i.e., a flat line). In general, we find that the more complex models have $\mathrm{ln}(B)\lt 0$, meaning that they are disfavored versus a flat line.

Assessing the level of TLSE contamination is critical to correctly interpreting high-precision transmission spectra of transiting exoplanets (Apai et al. 2018). This is particularly important for investigating features that could have a planetary or stellar origin, such the He 10833 Å line (Spake et al. 2018) or the Na and K alkali lines (Cauley et al. 2018). We have previously demonstrated this effect in the spectrum of GJ 1214b (Rackham et al. 2017), which has a featureless near-IR spectrum but a visible spectrum that is apparently offset below the near-IR transit depths. In this study, we found that for reasonable assumptions about the star's activity, the planetary transmission spectrum of GJ 1214b can be shown to be flat in the visible as well. We also detect the signature of stellar activity in one of our Magellan/IMACS transit data sets for WASP-19b (Espinoza et al. 2019).

WASP-4 is of slightly later spectral type than the Sun, and multiple spot-crossing events have been observed during previous transit observations of WASP-4b (Sanchis-Ojeda et al. 2011), as well as in two of our data sets. To assess the photometric variability, we retrieve four years of photometric monitoring data from the ASAS-SN online database¹² (Shappee et al. 2014; Kochanek et al. 2017). We remove outliers with a 5σ cut and create a Lomb–Scargle periodogram to search for periodic signals between 1 and 100 days, calculating the FAP as in Section 4.5. As shown in Figure 8, we find the largest peak at 22.4 days, and while the significance of the detection is only moderate, it corresponds well with the rotational period of 22.2 days that Sanchis-Ojeda et al. (2011) measure from consecutive spot-crossing events. This is to be expected if the photometric variability is dominated by the rotation of spots and faculae in and out of view. Modeling the variability as a sine curve and assuming a period of 22.4 days, we perform a least-squares fit to the phase-folded data to find a peak-to-trough variability amplitude of ∼6 mmag (0.6%). These results suggest that the surface of WASP-4 is moderately heterogeneous, so any analysis of the transmission spectrum of WASP-4b that ignores stellar contamination may be overly optimistic. Furthermore, since the spot-covering fraction cannot be determined from the variability alone (Rackham et al. 2018), the photosphere heterogeneity should be modeled alongside the atmosphere.

In the retrieval code described above, we characterize the heterogeneity with three parameters: T_occ, the effective temperature of the transit chord, T_het, the mean effective temperature of the heterogeneous features, and F_het, the fraction of the unocculted photosphere that is covered in spots. We then compute the effect on the observed transit depth following Section 4.7 and Equation (2), but replacing the area ratio ${({R}_{\mathrm{spot}}/{R}_{s})}^{2}$ with a covering fraction F_het.

We test more complex parameterizations that include multiple types of unocculted heterogeneities (e.g., Zhang et al. 2018), but we determine from the evidence that they are not warranted by the data.

Given the quality of our data, we find it necessary to fix some parameters of the heterogeneity model to reasonable values. In Table 5 we fix the spot temperature to match typical spots on the Sun (e.g., Solanki 2003, and references therein), and the transit chord temperature to the measured effective temperature of the photosphere.

Here, F_het is allowed to vary over all possible values and serves as a simple measure of the level of stellar contamination in the transmission spectrum. However, we note that F_het represents an enhancement over the effect of the large occulted spot in Transit 3, which has already been corrected for in Section 4.7.

In Figure 9 we present four of our atmosphere models with and without photosphere contamination. In Table 6 we list the Bayes factors for all 16 of the models tested.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 6.  Bayes Factors for the Full Suite of Atmosphere Models, Some of Which Are Shown in Figure 9

Model	$\mathrm{ln}{(B)}_{A}$	$\mathrm{ln}{(B)}_{A+P}$
Uniform opacity	0.0	−3.5
Clouds + Na	−1.8	−5.7
Clouds + K	−0.8	−5.0
Clouds + Na,K	−1.0	−5.4
Clouds + haze	−1.5	−5.4
Clouds + haze + Na	−1.9	−6.2
Clouds + haze + K	−1.9	−5.4
Clouds + haze + Na,K	−1.9	−6.2

Notes. The two columns represent atmosphere-only (A) and combined atmosphere–photosphere (A+P) models. Here, $\mathrm{ln}(B)$ is calculated relative to the flat-line model; since $\mathrm{ln}(B)\lt 0$ for the more complex models, they are all disfavored versus a flat line.

Download table as:  ASCIITypeset image

6.3.1. Favored Models

6.3.2. Constraints on Spot-covering Fraction

In Section 4.7 we detail a method for correcting the contamination due to the cool feature that was occulted by the planet during Transit 3. In this retrieval, however, we also include a varying parameter F_het and fixed parameter T_het to characterize the heterogeneity of the photosphere. We offer the following argument to justify our decision to visit the heterogeneity twice.

First, while the detailed structure of the occulted feature from Transit 3 is unclear, it is well described by a circular spot model with a constant temperature, and this is the same model that we use to correct for its contamination. By excluding this correction, we would not be leveraging all of the available information from our light curve. However, F_het must still be nonzero to describe the remaining unocculted features.

Second, we demonstrate that this correction is accounted for during the atmospheric retrieval. As a test, we fit the model of a uniform-opacity atmosphere with a heterogeneous photosphere (Figure 9, top right), this time fixing T_het = 5442 K to match the temperature of the occulted feature. When we fit the uncorrected combined spectrum, we find that the median of the posterior distribution for F_het is larger than for the corrected combined spectrum. The difference is consistent with a single spot of size R_spot/R_s ≈ 0.22. Similarly, when the retrieval is applied to the corrected and uncorrected spectra from Transit 3 only, the difference in the median value of F_het is consistent with a single spot of size R_spot/R_s ≈ 0.33. Both of these are similar to the spot size assumed in the correction, ${R}_{\mathrm{spot}}/{R}_{s}={0.27}_{-0.02}^{+0.03}$. This indicates that the effect of the correction is carried forward into the retrieval; we are not overcorrecting for the effect of the spot.

Huitson et al. (2017, hereafter H17) have previously published an optical transmission spectrum of WASP-4b from four combined transit observations with Gemini/GMOS, and they find that their data favor a uniform-opacity model. In this section, we assess the agreement of their results with our data from IMACS. In Figure 10 we compare the red and blue GMOS spectra with the combined spectrum from this work, and we find that the slopes of our spectra are generally in agreement.

Zoom In Zoom Out Reset image size

7.1.1. Absorption Features

May et al. (2018, hereafter M18) used Magellan/IMACS to study this target with the same grism as we used for Transits 1–3, and they also find no evidence for features in the transit spectrum. Figure 10 demonstrates that our results are in agreement.

7.2.1. Lack of a Spot-crossing Event

We combine the data from H17 and M18 with our own by subtracting the weighted mean from each spectrum, and then we repeat the retrieval analysis of Section 6. Table 7 displays the Bayes factor relative to a flat line for each of the 2x8 models. As before, we find that most atmosphere models are disfavored versus a uniform-opacity model by ln(B) = −1 to −2, but models including Na are slightly less disfavored than in Table 6. This indicates that there is more evidence for Na absorption in the combined data set than in ours alone, but still not enough to justify its inclusion.

Table 7.  Bayes Factors for the Full Suite of Atmosphere Models (Similar to Table 6) Where We Have Included Data from H17 and M18

Model	$\mathrm{ln}{(B)}_{A}$	$\mathrm{ln}{(B)}_{A+P}$
Uniform opacity	0.0	−3.8
Clouds + Na	−0.2	−5.5
Clouds + K	−1.3	−5.7
Clouds + Na,K	−0.7	−4.9
Clouds + haze	−1.7	−6.4
Clouds + haze + Na	−1.2	−6.0
Clouds + haze + K	−1.6	−6.3
Clouds + haze + Na,K	−1.6	−5.9

Download table as:  ASCIITypeset image

Models that include the photosphere are disfavored by $\mathrm{ln}(B)\sim -5.5$. In the case of a uniform-opacity atmosphere, the 95% upper limit on the spot-covering fraction is F_het < 1.8%, which is smaller than the value we report in Section 6.3.2. However, we caution that the red and blue GMOS spectra were observed during separate transits, which may mask the slope introduced by TLSE contamination.