ACCESS: Confirmation of no potassium in the atmosphere of WASP-31b
Chima D. McGruder, Mercedes Lopez-Morales, Nestor Espinoza, Benjamin V. Rackham, Daniel Apai, Andres Jordan, David J. Osip, Munazza K. Alam, Alex Bixel, Jonathan J. Fortney, Gregory W. Henry, James Kirk, Nikole K. Lewis, Florian Rodler, Ian C. Weaver
DDraft version September 22, 2020
Typeset using L A TEX twocolumn style in AASTeX62
ACCESS: Confirmation of no potassium in the atmosphere of WASP-31b
Chima D. McGruder, ∗ Mercedes López-Morales, Néstor Espinoza, Benjamin V. Rackham, † Dániel Apai,
4, 5, 6
Andrés Jordán,
7, 8
David J. Osip, Munazza K. Alam, ∗ Alex Bixel,
4, 6
Jonathan J. Fortney, Gregory W. Henry, James Kirk, Nikole K. Lewis, Florian Rodler, andIan C. Weaver Center for Astrophysics | Harvard & Smithsonian, 60 Garden St, Cambridge, MA 02138, USA Space Telescope Science Institute (STScI), 3700 San Martin Dr, Baltimore, MD 21218 Department of Earth, Atmospheric and Planetary Sciences, and Kavli Institute for Astrophysics and Space Research, MassachusettsInstitute of Technology, Cambridge, MA 02139, USA Department of Astronomy/Steward Observatory, The University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA Lunar and Planetary Laboratory, The University of Arizona, 1640 E. Univ. Blvd, Tucson, AZ 85721 Earths in Other Solar Systems Team, NASA Nexus for Exoplanet System Science Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Av. Diagonal las Torres 2640, Peñalolén, Santiago, Chile Millennium Institute for Astrophysics, Santiago, Chile Las Campanas Observatory, Carnegie Institution of Washington, Colina el Pino, Casilla 601 La Serena, Chile Department of Astronomy and Astrophysics, University of California, Santa Cruz, USA Center of Excellence in Information Systems, Tennessee State University, Nashville, TN 37209, USA Department of Astronomy and Carl Sagan Institute, Cornell University, 122 Sciences Drive, 14853, Ithaca, NY, USA European Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago de Chile, Chile (Accepted 14 September 2020 by AJ)
ABSTRACTWe present a new optical (400–950 nm) transmission spectrum of the hot Jupiter WASP-31b(M=0.48 M J ; R= 1.54 R J ; P=3.41 days), obtained by combining four transits observations. Thesetransits were observed with IMACS on the Magellan Baade Telescope at Las Campanas Observatoryas part of the ACCESS project. We investigate the presence of clouds/hazes in the upper atmosphereof this planet as well as the contribution of stellar activity on the observed features. In addition, wesearch for absorption features of the alkali elements Na I and K I, with particular focus on K I, forwhich there have been two previously published disagreeing results. Observations with HST/STISdetected K I, whereas ground-based low- and high-resolution observations did not. We use equilib-rium and non-equilibrium chemistry retrievals to explore the planetary and stellar parameter spaceof the system with our optical data combined with existing near-IR observations. Our best-fit modelis that with a scattering slope consistent with a Rayleigh slope ( α = 5 . +2 . − . ), high-altitude cloudsat a log cloud top pressure of -3.6 +2 . − . bars, and possible muted H O features. We find that our ob-servations support other ground-based claims of no K I. Clouds are likely why signals like H O areextremely muted and Na or K cannot be detected. We then juxtapose our Magellan/IMACS trans-mission spectrum with existing VLT/FORS2, HST/WFC3, HST/STIS, and
Spitzer observations tofurther constrain the optical-to-infrared atmospheric features of the planet. We find that a steeperscattering slope ( α = 8 . ± . ) is anchored by STIS wavelengths blueward of 400 nm and only theoriginal STIS observations show significant potassium signal. Keywords: planets and satellites: atmospheres — planets and satellites: individual (WASP-31b) —stars: activity — stars: starspots — techniques: spectroscopic
Corresponding author: Chima D. [email protected] ∗ NSF Graduate Research Fellow †
51 Pegasi b Fellow a r X i v : . [ a s t r o - ph . E P ] S e p McGruder et al. INTRODUCTIONWe are at the forefront of an exoplanet revolution withover 4200 planets discovered to date. Of these planets,over half have their fundamental properties (i.e. orbitalparameters, masses and radii) characterized, which givea basic understanding of their formation and interiorcomposition . However, to obtain a more in-depth un-derstanding of how planets form and evolve, it is neces-sary to study their atmospheres in detail.To date, transmission spectroscopy is the most suc-cessful approach to probe exoplanets’ atmospheres.However, most of the instruments currently used fortransmission spectroscopy were not designed with thisscience in mind, and unconstrained instrumental sys-tematics are a major limitation to data precision. Inthe case of ground-based observations, effects relatedto the Earth’s atmosphere are also a non-trivial factor.Because each telescope/instrument combination has itsown unique systematics, there is no specific prescriptionto handle all of these limitations, which can lead to dif-ferent interpretations of the same data based solely onthe detrending process. This problem has caused sev-eral claims and counterclaims in this field (e.g. Espinozaet al. 2019 vs. Sedaghati et al. 2017a; Diamond-Loweet al. 2018 vs. Southworth et al. 2017). WASP-31b isone such planet that has had a detection of K I (Singet al. 2015) that was later refuted (Gibson et al. 2017,2019). To exacerbate the situation, stellar activity canmimic planetary features (see e.g., Pont et al. 2013;McCullough et al. 2014; Espinoza et al. 2019), and tel-lurics might mask signals in ground-based observations(Gibson et al. 2017). Thus, the best approach in un-derstanding the atmospheric transmission spectrum ofan exoplanet is to repeatedly observe the planet androbustly detrend the data.WASP-31 is a low-activity (3-5 Gyr), F6V dwarf, withlow metallicity ([Fe/H] = -0.19) and an effective tem-perature of about 6300 K. WASP-31 (V = 11.7) is also avisual binary with a V ∼ (cid:48)(cid:48) away. WASP-31b is a transiting low-density( ρ = 0.129 ρ J ), inflated (R = 1.54 R J ) hot Jupiter withan orbital radius of 0.047 AU, a period of 3.4 days, amass of 0.48 M J , and a zero bond albedo equilibriumtemperature of 1575 K (Anderson et al. 2011). Thisplanet is a particularly good target for transmissionspectroscopy because its lower gravity leads to a largeratmospheric scale height and stronger signals .WASP-31b’s transmission spectrum was observedwith the Space Telescope Imaging Spectrograph (STIS) see, e.g., https://exoplanets.nasa.gov/exoplanet-catalog/ and Wide Field Camera 3 (WFC3) instruments on theHubble Space Telescope (HST). These observationswere combined with Spitzer photometry to constructa high signal-to-noise transmission spectrum. Analysisof those data indicated a Rayleigh scattering signa-ture at short wavelengths, a cloud-deck at longer wave-lengths, a strong potassium absorption feature (4.2 σ confidence, with 78 Å wide bins), but no sodium (Singet al. 2015). This is mostly consistent with ground-basedobservations later made with FORS2 on the VLT, asidefrom the detection of K I (Gibson et al. 2017). Gibsonet al. re-analyzed the HST observations using Gaussianprocesses (GPs) to correct for systematics and foundthat the potassium detection significance decreased to2.5 σ . Additional high-resolution observations (R >80000) with the VLT’s UV-Visual Echelle Spectrograph(UVES) did not detect sodium or potassium, supportingthe idea that HST observations also might be suscepti-ble to spurious signals produced by systematics (Gibsonet al. 2019). However, because the detection of theK I feature from the ground can be contaminated byEarthâĂŹs O telluric features, it is worth further in-vestigating these results.Here we present four new transit observations ofWASP-31b with spectral resolutions 300 < R < 1300and wavelength coverage between 0.4–0.97 µ m, obtainedbetween 2013 and 2019 as part of ACCESS . ACCESSutilizes the multi-object spectrograph (MOS) IMACS(Dressler et al. 2011) on the 6.5-meter Magellan BaadeTelescope at Las Campanas Observatory, which allowsfor simultaneous collection of spectra from the exoplanethost star and many comparison stars to enable bettersystematic corrections. To date, ACCESS has producedtransmission spectra of WASP-6b (Jordán et al. 2013),GJ 1214b (Rackham et al. 2017), WASP-4b (Bixel et al.2019), WASP-19b (Espinoza et al. 2019), and WASP-43b (Weaver et al. 2020). To these we add here opticalmeasurements of WASP-31b, which we combine withpreviously published HST/WFC3 and Spitzer
IR datato produce a new transmission spectrum with coveragebetween 0.4 and 5.0 µ m.This paper is structured as follows: In Section 2, wepresent our Magellan/IMACS observations of WASP-31b in addition to photometric monitoring observationsto constrain stellar variability. In Section 3, we dis-cuss our transit light curve extraction process, binningschemes, and detrending techniques. We outline the pro-duction of our final transmission spectrum and present ACCESS is a multi-institutional collaboration aiming to pro-duce a large, homogeneous library of optical spectra of exoplanetatmospheres. http://project-access.space/
CCESS: WASP-31b
3a qualitative comparison of this spectrum to data fromthe literature in Section 4. In Section 5, we introducethe retrievals used to analyze the transmission spectraand report our best-fit results compared to best fits ob-tained from previous data. Interpretations of the re-trieval results in context of WASP-31b’s overall atmo-spheric properties are discussed in Section 6. We con-clude and summarize in Section 7. OBSERVATIONS2.1.
General Setup
We observed four transits of WASP-31b on UTdates UT130226, UT130425, UT140222, and UT190314(UTYYMMDD) with IMACS on the 6.5-meter BaadeMagellan Telescope at Las Campanas Observatory inChile. The large time-lapse between the first threetransits and the last, is because transit UT130226 couldnot be used for the transmission spectrum (as discussedin appendix B), which we only found to be the case laterin the analysis process.IMACS is a wide-field imager and spectrograph withtwo 8Kx8K CCD mosaic cameras at the f/2 and f/4foci (Dressler et al. 2011). Each camera has a multi-object spectroscopy (MOS) mode, allowing for customslit masks designs. The f/2 camera has a large, . (cid:48) field of view (FoV) and resolving powers up to R ∼ . (cid:48) x . (cid:48) )but higher spectral resolution (max R ∼ ), as described in detail in Rackham et al. (2017).The comparison stars selected for each observation aresummarized in Table 1.We designed masks with slits wide enough to preventslit losses in typical seeing conditions at LCO but nar- color distance, D, is defined in Rackham et al. (2017) as D = (cid:112) [( B − V ) c − ( B − V ) t ] + [( J − K ) c − ( J − K ) t ] where t and c represent the target and comparison star magnitudes, respec-tively. row enough to prevent contamination from nearby stars.This resulted in slit widths of (cid:48)(cid:48) and (cid:48)(cid:48) for the f/2and f/4 mode, respectively. The lengths of the slitswere (cid:48)(cid:48) for f/2 and (cid:48)(cid:48) for f/4 to adequately samplethe sky background. We used a 300 line/mm gratingat blaze angle of 4.3 ◦ for the f/4 observations and a300 line/mm grism at blaze angle of 17.5 ◦ for the f/2observation. These setups provided a wavelength cov-erage of 0.4–0.95 µ m and 0.45–0.97 µ m for the f/4 andf/2 observations, respectively. Finally, to reduce readouttime and improve the duty cycle of the observations, webinned the detectors by 2 × × spatial direc-tions) for the UT130266 transit. Transit UT130266 wasthe first observation done as part of ACCESS, while wewere still experimenting with settings. After this firsttry we realized that seeing conditions at Las CampanasObservatory are often very good (typically 0.7 arcsec orbetter), so 2 × × ∼
900 (FWHM dispersion of8.08Å).
Table 1.
Target and comparison star magnitudes and co-ordinates from the UCAC4 catalog (Zacharias et al. 2013).We used comparisons 2–4 for the white-light fit of epochUT130226, 1–4 for epochs UT130425 and UT140222, andstars 2, 5–7 for epoch UT190314. Comparison 1 for epochUT130226 was saturated in most bins and not used for thistransit.Star RA Dec B V J KWASP-31 11:17:45.4 -19:03:17.2 12.4 11.9 10.9 10.71 11:17:45.3 -19:07:25.9 11.8 11.1 10.0 9.62 11:17:32.0 -19:08:37.0 13.1 12.6 11.6 11.43 11:17:38.1 -18:59:19.6 13.7 13.0 11.9 11.54 11:17:44.3 -18:58:06.5 14.7 13.8 11.9 11.45 11:17:10.8 -19:14:26.0 12.5 11.9 11.0 10.76 11:16:35.7 -19:01:35.6 13.2 11.9 9.7 8.97 11:18:23.2 -18:58:45.8 13.3 12.8 11.9 11.6
IMACS Data Collection
The transits of WASP-31b last 2.6 hours, so to coverthe full transits and obtain enough baseline coverage toextract precise transit depths, we monitored the sys-tem during each transit epoch for about 5-hour contin-uous time blocks. We limited observations to epochs
McGruder et al. when WASP-31 was at airmasses below 1.7 during thefull observation window to minimize atmospheric re-fraction effects in the data. Exposure times were ad-justed from 27 to 60 seconds to provide maximum countrates of 25,000–35,000 in analog-to-digital units (ADU;gain= 0.56 e − /ADU for f/4 setup, 1.0 e − /ADU forf/2 setup), which was well within the CCD’s linearitylimit (Bixel et al. 2019). More details of specific transitobservations are given in Table 2.For each night we collected bias frames, high-SNRquartz lamp flats, and calibration arcs with HeNeArlamps, each with the same binning scheme as the sci-ence observations. From the biases we found that thelevels were essentially constant across the detector, so weadopted a constant bias level on each science frame us-ing the median of their corresponding overscan region.We did not take dark frames because dark counts arenegligible for our exposure times (Rackham et al. 2017).We took series of ∼
20, 1 sec flat frames before and afterthe science frames with the same configuration as thescience images. The exposure times of the flats wereset to give a count rate of about 35,000 ADU or underin order to stay within the linear regime of the CCDs(Bixel et al. 2019). However, during the data reduction(see Section 3.1), we found that flat-field corrections in-troduced additional noise in the data, similar to whatwas found with previous ACCESS observations (Rack-ham et al. 2017; Bixel et al. 2019; Espinoza et al. 2019;Weaver et al. 2020). Therefore, we decided not to applythe flat-field correction to any of our data sets. Finally,we took exposures with He, Ne, and Ar lamps using cal-ibration masks identical to the science masks but with0.5 (cid:48)(cid:48) slit widths to increase the spectral resolution andprevent saturation of the CCDs. These lamps were usedto wavelength-calibrate the spectra, as described in Sec-tion 3.1. 2.3.
Photometric Monitoring
We analyzed photometric time series of the host starWASP-31 to constrain potential stellar contributionsto the transmission spectrum of WASP-31b using ob-servations from the Tennessee State University (TSU)Celestron 14-inch (C14) Automated Imaging Telescope(AIT) at Fairborn Observatory (Henry 1999), the Tran-siting Exoplanet Survey Satellite (TESS, Ricker et al.2014), and the All-Sky Automated Survey for Super-novae (ASAS-SN, Kochanek et al. 2017). The combinedtime baseline of these observations cover all our transitepochs, as shown in Figure 2. However, the photometricprecision varies between surveys, so we decided to ana-lyze each observing season for each survey separately.2.3.1.
AIT Analysis
Figure 1.
Median extracted spectra of WASP-31 from tran-sits UT130226 (orange), UT130425 (green), UT140222 (red),and UT190314 (blue). The shaded regions of the same colorextending past the median lines are the 1- σ range of countsextracted for that night. The difference in counts betweenthe f/4 observations are because of the exposure times. Dif-ferent throughputs with the f/4 and f/2 setups explain thedifferent shape of the spectra from UT190314 and from thef/4 observations. Dark shaded regions indicate chip gaps; ascan be seen there are 3 gaps for f/4 observations and usabledata starts at .48 µ m for the f/2 mode. TSU’s AIT acquired 660 photometric observationsthrough a Cousins R filter spanning 6 observing sea-sons (2011-12 through 2016-17). Standard deviations ofa single observation from their seasonal means range be-tween 3.2 and 4.2 mmag. This is near the limit of mea-surement precision with the AIT, as determined fromthe comparison stars in the field. The seasonal-meandifferential magnitudes agree to a standard deviation ofonly 0.88 mmag, consistent with the absence of year-to-year variability and the low activity level of WASP-31. Periodogram analysis of each individual season re-sulted in the detection of significant periodicity only inthe 2012–13 season. We detected a small peak-to-troughamplitude of 4 mmag with a period of 8.60 days, corre-sponding to the star’s rotational period of 8.62 days asmeasured by (Sing et al. 2015). This suggests WASP-31is a relatively inactive star .2.3.2. TESS Analysis
WASP-31 was observed in TESS Sector 9 (2/28/19-3/26/19) for 27 days with 2-minute cadence observa-tions. To understand stellar variability, we used the out-of-transit Pre-search Data Conditioning Simple Aper- See Sing et al. (2015) for further details of AIT observations,reduction, and analysis.
CCESS: WASP-31b Table 2.
Observing log for WASP-31b data sets from
Magellan/IMACS . The resolutions were calculated from the FWHM,which was estimated using the second moment method (Markevich & Gertner 1989). The average seeing conditions at LCO are . (cid:48)(cid:48) − . (cid:48)(cid:48) (Thomas-Osip et al. 2008) and is the resolution limiting factor.Transit Date Obs. Start/ Airmass Exposure Readout + Frames binning camera min/max(UTC) End (UTC) Times (s) Overhead (s) resolution2013 Feb 27 02:58/8:20 1.02–1.3 40 29 283 2 × × × × ture Photometry (PDCSAP) light curve. After maskingthe transits of WASP-31b and binning to ∼ ASAS-SN Analysis
ASAS-SN monitors the entire visible sky (V-band fil-ter) with 24 telescopes all over the world. There are262 observations of WASP-31 from December 2013 toNovember 2018 with two ASAS-SN cameras (bg, be) atthe Cerro Tololo International Observatory, Chile andone camera (bc) at the Haleakala Observatory, Hawaii.The photometric precision of the ASAS-SN data is lessthan that of the AIT and TESS data, as seen in Figure 2.Thus, we cannot extract any evidence of stellar variabil-ity from the ASAS-SN light curves. Additionally, thereis little ASAS-SN coverage of our transits; as such, weused the AIT and TESS data to constrain WASP-31’slevel of activity. The AIT and TESS phase-folded lightcurves are shown in Figure 3.Our conclusion from the analysis of the AIT, TESSand ASAS-SN datasets is that WASP-31b is a typicalphotometrically quiet F-type star. However, there is ob-servational evidence that photometrically quiet F dwarfsare not completely inactive (e.g, López-Morales et al.2016). Furthermore, photometric monitoring alone isnot enough to fully characterize stellar activity levels, as axisymmetrically distributed active regions do not con-tribute to rotational variability (Rackham et al. 2019).Therefore, we still consider stellar activity in our re-trieval analyses in Section 5. DATA REDUCTION & LIGHT CURVEANALYSES3.1.
Reduction Pipeline
We reduced the raw IMACS data using the pipelinedescribed in previous ACCESS papers (Jordán et al.2013; Rackham et al. 2017; Espinoza et al. 2019; Bixelet al. 2019; Weaver et al. 2020). It performs standardbias and flat calibration, bad pixel and cosmic ray cor-rection, sky subtraction, spectrum extraction, and wave-length calibration. The pipeline produces sets of ex-tracted, wavelength-calibrated spectra for each epoch,for the target and each of the comparison stars in thefield.The wavelength calibration was done by first fittingLorentzian profiles to each spectral line on images takenwith the calibration masks. Next, using these pixel po-sitions and the known vacuum wavelengths, the wave-length solution for each spectrum was found by an itera-tive process in which a sixth order polynomial was fittedto the wavelengths as a function of pixel position, thedata point with the greatest deviation from the fit wasremoved, and the process was repeated until the rootmean square error value of the fit was less than 2 kms − ( ∼ ∼ µ m, which we excluded.The extracted spectra are then integrated across theentire wavelength range to produce white light curves orintegrated over a narrow wavelength range to produce McGruder et al.
Figure 2.
Photometric monitoring of WASP-31 with AIT (black), ASAS-SN (cyan), and TESS (blue). The times of the transitsobserved with Magellan are shown as dash dotted vertical lines. Given the differing observing bands of ASAS-SN (V), AIT (R),and TESS (650–1050 nm), we limit our comparison to relative magnitude changes. Thus, the values plotted are the differencesfrom the mean magnitude of each survey. All observations suggests that the star is relatively inactive, however, this cannotconfidently be constrained with photometry alone.
Figure 3.
Phase curves of the AIT and TESS observations, showing low-amplitude periodicity in periodogram analyses. AITdata (left) show a 8.6-day periodicity with an amplitude of 0.004 mag. For the TESS data (right), the transparent gray pointsare all ∼ the binned light curves. These integrated white lightcurves are shown in the top panel of Figure 4.3.2. Light Curve Detrending
We first modeled the light curves following the ap-proach of Jordán et al. (2013), using principal compo-nent analysis (PCA). For transit data (time-series mea-surements) we assume that the light curve L k ( t ) for agiven comparison star k is a linear combination of a setof signals s i ( t ) , which represent the different instrumen- tal and atmospheric effects on the light curves, i.e., L k ( t ) = K (cid:88) i =1 A k,i s i ( t ) , (1)where A k,i is the weight of each signal s i ( t ) in each lightcurve. CCESS: WASP-31b L k ( t ) as a row in a matrix D , weperform singular value decomposition . This returns theeigenvectors (cid:126)e i and eigenvalues λ i of the matrix. We ob-tain s i ( t ) from the product of the matrix of eigenvectorsand the original data matrix D . Additionally, the im-portance of each signal for reconstructing D is given bythe eigenvalues λ i (largest λ i , most important (cid:126)e i , etc.).This allows for the sequential fitting of those signals toour target star in order to find the simplest model thatpreserves the most important information present in thecomparison stars. For further details of PCA, we referthe reader to Jolliffe (2002), or Espinoza (2017) for itsusage in high signal-to-noise time-series photometry. Weimplemented PCA with the same procedure as Jordánet al. (2013) and Rackham et al. (2017) and found 2,4, 2, and 3 of the 4 PCA components to be the optimalnumber of components for each transit, chronologically .However, as can be seen in the middle panel of Fig-ure 4, using this PCA method alone was not sufficient toproperly correct all systematics. This is because thereare likely additional systematics that uniquely affectWASP-31 due to, e.g., color, detector, and instrumen-tal differences with respect to the comparison stars. Toaccount for this, we used Gaussian process (GP) regres-sion to further model the target light curve in conjunc-tion with PCA. A GP is a machine learning techniquefor non-parametric regression modeling (Rasmussen &Williams 2005), and has proven effective at detrendingradial velocity data (e.g, Aigrain et al. 2012) and trans-mission spectroscopy (e.g, Gibson et al. 2012). For ourGP analysis, we define a joint probability distributionof the form N [0 , Σ] , where the covariance function ( Σ )is defined as K SE ( x i , x j ) + σ w δ i,j . Here, σ w , δ i,j , and K SE ( x i , x j ) are a jitter term, the Kroenecker delta func-tion, and our kernel respectively.A multidimensional squared-exponential kernel wasused for our GP regression: K SE ( x i , x j ) = σ GP exp (cid:32) − D (cid:88) d =1 α d ( x d,i − x d,j ) (cid:33) , (2)where σ GP is the amplitude of the GP and α d are the in-verse (squared) length-scales of each of the componentsof the GP. The x i vectors here have components x d,i ,where each i denotes a time-stamp and where each dcorresponds to a different external parameter. The mo-tivation for using the SE kernel is that it is a relativelysmooth function and self-regulates the importance of the Cov t [ s i ( t ) s j ( t ) ]=0 assuming that the signals are uncorrelated. The optimal number of PCA components were determinedwith the k-fold cross-validation procedure (Hastie et al. 2007). external parameters ( x ) (Micchelli et al. 2006; Duvenaud2014).Because the joint probability of the data can be de-fined by a multivariate Gaussian using GPs, then thelog marginal likelihood is known and given by: ln L = − r T Σ − r −
12 ln | Σ | − N π ) . (3)Here N is the number of measurements of the data. Weused george (Foreman-Mackey et al. 2014) to eval-uate the marginalized likelihood of our dataset and PyMultiNest (Buchner et al. 2014) (a python wrapperof MultiNest (Feroz et al. 2009)) to sample the hyper-parameter space. In our analyses, we found including5 external parameters x useful to model the PCA cor-rected data with GPs. These time-dependent parame-ters were: i) the sky flux, ii) the variation of the FWHMof the spectra, iii) the airmass, iv) the position of thecentral pixel trace (perpendicular to the dispersion axis),and v) the drift of the wavelength solution.The first 3 parameters are proxies of how Earth’s at-mosphere is affecting the stellar spectra. Over all ex-posures the FWHM ranges from 4.3 to 17.1 pixels andthe airmass varied from 1.7 to 1.0 over the course ofa night for each transit. The last two parameters areproxies for how slight physical variations in the instru-ments/detector could affect the spectra. Ultimately,exactly how these parameters affect the spectra is notwell understood, but GPs are optimal at finding correla-tions between parameters and observables (Rasmussen& Williams 2005).The next task is to properly integrate the PCA, GP,and inherent transit model into one function defining ourlight curves. Our particular model used the magnitudeof the target and comparison stars and was also used byYan et al. (2020) to detrend ground-based spectroscopictransit data: M k ( t ) = c k + N k (cid:88) i =1 A k,i s i ( t ) − . log T ( t | φ ) + (cid:15), (4)where M k ( t )) is the (mean-subtracted) magnitude of thetarget star in the k th model, c k is a magnitude offset, N k is the number of PCA signals s i ( t ) , A i,k is the weight foreach signal in each models, T ( t | φ ) is the transit modelwith parameters φ , and (cid:15) is a stochastic componenthere modeled as a GP. https://github.com/dfm/george https://github.com/JohannesBuchner/PyMultiNest We used the Python package batman (Kreidberg 2015) to pro-duce the analytic transit model.
McGruder et al.
Figure 4. Top:
Raw white light curves for our Magellan/IMACS observations directly from ACCESS’s custom pipeline. Blackis WASP-31 and transparent red, green, blue, magenta, cyan, grey, and yellow are comparison stars 1-7, respectively. SeeTable 1 for more information about the stars.
Middle:
The PCA-corrected light curves are shown in green. The magentalines are GP fits on the light curves to correct for systematics using the 5 external parameters. This figure depicts that PCAcorrection alone is not sufficient to remove all systematics, and as such model averaging of PCA and GPs are justified.
Bottom:
Final white-light curve with the model-averaged PCA and GP systematics removed (blue points), along with the correspondingresiduals (red points). The best-fitting transit light curve given our model averaging procedure is depicted as a solid black line.The value of σ given in each panel is the standard deviation of the residuals in ppm. The middle panel is not an intermediatestep to arrive at the bottom panel, given that the bottom panel is obtained directly when simultaneously doing PCA and GPmodel averaging of the data in the first panel. As can be seen in the right column, tranist UT130226 has little pre-transit data,which might cause systematics in its transmission spectra. This is explored further in Appendix B. Though PCA provides a weighting scheme for the im-portance level of the signals s i ( t ) , the determination ofwhich components improves our light curve model is stillnon-trivial. In order to incorporate our ignorance ofthe functional form of nuisance signals into our modelfit, we used the model averaging technique outlined byGibson (2014). For our case, we averaged the posteriordistributions of each M k ( t )) , where the posteriors wereexplored using PyMultiNest . Our final model-averagedwhite light curves can be seen in the bottom row of Fig-ure 4. 3.3.
Light Curve Fitting
Using an improper limb-darkening (LD) law can intro-duce biases in the final light curve fit(Espinoza & Jordán2016), and as such we used the open source package ld-exosim to determine the most appropriate LD lawfor our transit fits. The square root law was found to bebest for all white light fits, given the data quality. Weused the same law for the spectroscopic fits, describedin Section 4. https://github.com/nespinoza/ld-exosim CCESS: WASP-31b σ -clipping to each of the white lightcurves to remove data points that deviated more than 3 σ from the model. This clipping was applied after obtain-ing an initial detrended light curve and correspondingbest-fit model with all available data and lead to 3, 5, 0,and 1 data points being removed from epochs UT130226,UT130425, UT140222, and UT190314, respectively.We then used the 3- σ -clipped light curves to find thebest transit parameters. The transit parameters used tofit the light curves were LD coefficients ( q , q ), planetorbital period ( P ), semi-major axis (in terms of star ra-dius; a/R s ), the planet-to-star radius ratio ( p ), impactparameter ( b ), mid-transit time ( t ), inclination ( i ), ec-centricity ( e ), and the argument of periapsis ( ω ). Wekept e and ω fixed to 0 and 90 ◦ (assuming a circularorbit), following the results of Anderson et al. (2011).We placed Gaussian priors on all other parameters us-ing the results of Anderson et al. (2011), then fit forthem following the procedure described in Section 3.2.The best-fit white light curve orbital parameters for eachepoch are summarized in Table 3. TRANSMISSION SPECTRATo produce the optical transmission spectrum ofWASP-31b, we generated transit light curves in thesame manner as described in Section 3, but using a setof narrower wavelength bins. We propagated the flagged3 σ outliers from the white-light analysis to each binnedlight curve. Our binning scheme was devised so it opti-mized spectro-photometric precision, while still lettingus probe for atmospheric features, such as a scatteringslope, and sodium and potassium lines. The average binwidths were 187 Å, where the lower-throughput, redderbins were made larger, and the regions around the Na Iand K I features were 60 Å wide centered on each dou-blet (5892.9 Å, 7681.2 Å). We found 60 Å to be enoughto encompass both doublet lines of each feature, provid-ing high signal from potential features while minimizingsuppression of signal from the surrounding continuum.This is also comparable to the 78 Å bin widths that Singet al. (2015) used to report a 4.2 σ K I detection.We employed the same analysis procedures as thewhite light curves (see Section 3.3) but with P , a/R s , b , t , and i held fixed to values obtained from the corre-sponding white light curve fits, leaving only q , q , and p as free parameters. The detrended light curves, theirbest-fit model, and residuals for each bin of each epochare shown in Figures 8–11 in Appendix A. Tables 6 and7 of Appendix A list the values of p obtained for eachbin. 4.1. Combined Transmission Spectrum
The average precision of each bin for each spectrumis 0.007 R p /R s . This is larger than the ∼ R p /R s expected signal from Na and K features for this planet(Sing et al. 2015). Therefore, we combined the trans-mission spectra from epochs UT130425, UT140222, andUT190314 to achieve higher precision. The UT130226transit was not used in our combined data set because itstransmission spectrum was inconsistent with the others,perhaps because of stellar activity or insufficient baselinebefore ingress. We discuss this further in Appendix B.Our spectra were combined by first applying an offsetto each spectrum equal to the difference between theirwhite light curve depths and the overall mean whitelight curve depth. We then combined the measurementsfor each spectroscopic bin by applying a weighted av-erage using each point’s R p /R s errorbars as weights.The combined transmission spectrum, shown as blackdiamonds in Figure 5, has an average precision per binof 0.0033 R p /R s . We emphasize that when combiningtransmission spectra of different epochs, as we did, in-formation about how the star and the terminator of theplanet’s atmosphere are behaving at a given epoch islost, leaving only persistent (average) behaviors in theplanet-star system.As can be seen in Table 3, the best-fit planet radiusvaries per epoch. These differences in our white lightcurve transit depths can be explained by the differentlocations of wavelength gaps in the detector for boththe comparison stars and WASP-31, as shown in Fig-ure 1, which causes the integrated flux to be slightlydifferent per epoch. Additionally, some epochs use dif-ferent comparison stars as shown in Table 1, which causedifferences in white light curve flux ratios when divid-ing the flux of the targets by the flux of the comparisonstars. We tested this hypothesis by using one matchingcomparison and a fixed wavelength range with no gapsfor all datasets and found consistent transit depths forall three epochs. However, this approach significantlyreduces the precision of the white light curves and thespectro-photometric bins, so we resorted back to apply-ing offsets between epochs in order to use the full wave-length coverages and all comparisons. This approach isjustified because the most important information fromthe transmission spectra is the relative changes in depthbetween wavelength bins. In addition, to address thisissue, we fit for offsets between our combined transmis-sion spectrum and other datasets in our atmosphericretrieval analysis, as described in Section 5.Figure 5 shows the individual transmission spectra foreach transit epoch, together with our combined trans-mission spectrum. Overlaid on the combined spectrumis a two-part fit with a scattering slope bluewards of0 McGruder et al.
Table 3.
Fitted white light curve values. These are calculated using all integrated flux and available comparisons, uniqueto each epoch. This causes variation in depths per night, but improves the precision of the light curve parameters needed forbinned light curve fits.parameter definition UT130226 UT130425 UT140222 UT190314 p planet radius/star radius . ± .
002 0 . ± .
002 0 . +0 . − . . ± . t − . e mid-transit (JD) . +0 . − . . ± . . +0 . − . . ± . P period (days) . ± . . ± . . ± . . ± . a/R s semi-major axis/star radius . +0 . − . . ± .
09 8 . ± .
09 8 . +0 . − . b impact parameter . +0 . − . . +0 . − . . ± .
006 0 . +0 . − . i inclination (degrees) . ± . . ± . . ± . . ± . q LD coeff 1 . ± . . +0 . − . . ± . . +0 . − . q LD coeff 2 . ± . . +0 . − . . +0 . − . . +0 . − . µ m and a featureless spectrum redwards of 0.55 µ m,similar to fits done by Sing et al. (2015) and Gibson et al.(2017). The fit is favored over a completely featurelessspectrum ( χ = 50.9) with a χ of 46.2.4.2. Previous Optical Transmission Spectra
All available low-resolution optical transmission spec-tra of WASP-31b are shown in Figure 6. When quali-tatively comparing the Magellan/IMACS data with theVLT/FORS2 (Gibson et al. 2017) and HST/STIS ob-servations (Sing et al. 2015), we see that all 3 spectracan be modeled with a scattering slope in the blue anda featureless spectrum beyond. The simple models inFigure 6 were created the same way as in Figure 5, justas a comparison with Sing et al. (2015)’s model. How-ever, in our retrieval analysis (Section 5.2), we do notuse a two part model, but rather obtain scattering slopeswith one continuous atmospheric model over the wholespectrum. This is also consistent with the findings ofBarstow et al. (2016), who found that a "split" spectrumwith the HST data is not significantly favored comparedto a continuous model. Like the FORS2 observations,our IMACS observations show no signs of potassium orsodium absorption. This is contrary to the detection ofpotassium with HST/STIS. Though the overall shapesof the spectra agree with each other (aside from potas-sium absorption), we elected not to combine previousstudies with our own in the subsequent analysis. Thiswas because we are uncertain of the biases introduced tothe data by the various detrending techniques for eachanalysis, and combining the spectra may lead to misin-terpretation of signals lost or produced by compoundedbiases. This rationale is supported by the analysis ofGibson et al. (2017), which found significantly differentresults for the same dataset depending on the detrend-ing method. However, we did combine the IMACS op-tical transmission spectrum with near-IR HST/WFC3 spectra and mid-IR (
Spitzer /IRAC) data (Sing et al.2015), because there is only one observation/analysisfor each and including them will help better characterizeWASP-31b’s atmospheric structure. In total, the com-bined transmission spectrum covers a wavelength rangefrom 0.4 µ m to 5.06 µ m, with average bin sizes of 187Åand 200Å for the optical(Magellan/IMACS) and near-IR (HST/WFC3) data. ATMOSPHERIC RETRIEVAL ANALYSESWe use atmospheric retrievals to provide a quanti-tative measurement of the existence of features in thetransmission spectra. Given that results can depend onthe priors and the specific retrievals used (e.g., Kirk et al.2019), we employed two independent retrievals to im-prove data interpretation:
PLATON (Zhang et al. 2019),and Exoretrievals (Espinoza et al. 2019).
PLATON and
Exoretrievals complement each otherin several ways:
Exoretrievals considers a smallnumber of user defined molecules, but
PLATON con-siders a large number (33) of pre-determined molecules.
Exoretrievals performs free retrievals, allowing forabundances of individual molecules to be retrievedwith non-equilibrium chemistry, whereas
PLATON as-sumes equilibrium chemistry throughout the atmo-sphere. Thus, using both retrievals gives a more thor-ough understanding of the planet. Additionally, findingagreeing results in overlapping parameters, strengthenssupport of accurate interpretation of the data, giventhat they use very different fundamental assumptions.In detail,
PLATON uses the planet’s limb tempera-ture ( T p ), atmospheric metallicity (Z/Z (cid:12) ), C/O ratio,cloud-top pressure ( P ), a scattering slope ( α ), andan offset term between instruments to characterize ob- https://github.com/ideasrule/platon CCESS: WASP-31b Figure 5.
All four transmission spectra of WASP-31b taken by Magellan/IMACS, before an offset is applied, (faded colors)and the weighted average of the last 3 spectra (black). The first transit, UT130226 (transparent orange), significantly deviatesfrom the other 3 spectra and has its highest ∆ R p /R s of 0.0385, but cannot be seen in this figure because emphasis is placed onthe other three consistent spectra. Transit UT130226 is likely an outlier because of insufficient baseline during observation orstellar activity (see Appendix B). The purple line is a two part fit with a scattering slope bluewards of 5500Å and a featurelessspectrum redwards of 5500Å, similar to fits done by (Sing et al. 2015) and (Gibson et al. 2017). The highlighted gray and yellowregions are centered around the Na I and K I lines, respectively. served transmission spectra using equilibrium chemistry. PLATON ’s sources of line lists for all 33 molecules con-sidered are shown in Table 2 of (Lupu et al. 2014),which mostly, but not exclusively, comes from HITRAN.
PLATON includes collision induced absorption, where allcoefficients come from HITRAN.
PLATON also simultane-ously fits for stellar activity using the occulted temper-ature of the star ( T star ), the temperature of unoccultedspots ( T spot ), and their area ( f spot ). PLATON uses thenested sampling algorithm dynesty (Speagle 2019) toexplore the parameter space. Exoretrievals uses a semi-analytical formalism tomodel the planet with an isothermal, isobaric, atmo-sphere with non-equilibrium chemistry (see Heng &Kitzmann 2017). The cross sections we consideredwith
Exoretrievals were CO, CO2, CH4, TiO, NH3,FeH, HCN, H2O, Na, and K where the line list of themolecules came from HITEMP (Rothman et al. 2010),ExoMol (Yurchenko et al. 2013; Tennyson et al. 2016),and TiO was synthesized using the latest (2012 ver-sion) line lists calculated by B. Plez, obtained fromVALD (Ryabchikova et al. 2015)). The Na and K cross-sections were determined from analytical Lorentzian-profiled doublet shapes used in MacDonald & Mad-husudhan (2017). Unlike
PLATON , Exoretrievals does https://github.com/joshspeagle/dynesty not include collision induced absorption. However, itdoes also account for stellar photospheric heterogeni-ties by modeling the occulted and unocculted regions ofthe stellar photosphere with PHOENIX models (Husseret al. 2013), following the approach of Rackham et al.(2018, 2019). The retrieval has three parameters tocharacterize the stellar surface: the temperature of theregion occulted by the planet ( T occ ), the temperatureof unocculted heterogeneities ( T het ), which can be hot-ter or colder than the mean photosphere, and the het-erogeneity covering fraction ( f het ). The planet’s atmo-sphere is characterized by a factor ( f ) multiplied to theplanetary radius to find a reference radius at which theatmosphere is optically thick (i.e. R = f R p ), the cor-responding pressure where the atmosphere is opticallythick ( P ), the limb atmospheric temperature (T p ), mix-ing ratios of specified species (i.e. Na, K, and H O),offset terms for different instruments, and two parame-ters to characterize the cloud/haze extinction ( κ ). Thetwo cloud/haze parameters are a Rayleigh-enhancementfactor ( a ) and a cloud/haze power law ( γ ), and theyquantify extinction with: κ ( λ ) = aσ ( λ/λ ) γ , (5)where σ = 5 . × − m and λ = 350nm (Sedaghatiet al. 2017b; MacDonald & Madhusudhan 2017). Thecode uses nested sampling with PyMultiNest to mea-sure the evidence for each of the retrieved models. The2
McGruder et al.
Figure 6.
Transmission spectra of WASP-31b taken by Magellan/IMACS (excluding epoch UT130226, black), VLT/FORS2(purple, Gibson et al. 2017), HST/STIS (cyan, Sing et al. 2015), and combined FORS2 and reanalyzed STIS data (magenta,Gibson et al. 2017), with vertical offsets for clarity. Each spectrum has a two-part fit of a scattering slope blueward of 5500Åand a featureless spectrum redward of 5500 Å, excluding the possible K signal centered at 7682.5 Å for the HST/STIS data,overlaid as dashed lines. All three spectra have signs of a blue scattering slope and featureless optical spectrum. The highlightedgray and yellow regions are centered around the Na I and K I lines, respectively. priors we used for both retrievals are shown in Table 4.Because
Exoretrievals allows for individual moleculesto be tested, we used the results of
Exoretrievals todetermine the presence of water, potassium, or sodium.5.1.
Combined Magellan Data
We run
PLATON and
Exoretrievals on the combinedIMACS, HST/WFC3, and
Spitzer data. For each re-trieval run, we fit for an offset between the optical dataand the IR data but no offset between the
Spitzer andthe HST/WFC3. The
Spitzer and HST/WFC3 data wastaken as a single IR dataset, where Sing et al. (2015) an-alyzed them both using the same ephemeris and systemparameters. Using the retrievals we tested for differentscenarios: 1) a featureless exoplanet atmosphere, whereall potential signals are covered by high-altitude clouds,2) a spectrum with any signal detected in the spectrumcoming from active regions on the surface of the star,3) an atmosphere with scatterers, and 4) an atmospherewith scatterers, and stellar activity contamination. Ad-ditionally, we tested
Exoretrievals ’ models with thesame 4 scenarios, but including only H O, Na, or K, in-dividually and simultaneously in order to isolate theiratmospheric effects.Table 5 summarizes the Bayesian evidences of eachscenario with respect to a featureless spectrum ( ∆ ln Z ).Following Trotta (2008) and Benneke & Seager (2013),the frequentist significance of those ∆ ln Z values scalesas: | ∆ ln Z | of 0 to 2.5 is inconclusive with < 2.7 σ sup- port for the higher evidence model, | ∆ ln Z | of 2.5 to 5corresponds to a moderately significant detection of 2.7 σ to 3.6 σ , and | ∆ ln Z | ≥ Exoretrievals could also fit for CO, CO , CH , TiO,NH , FeH, and HCN. Given the temperature(Fortneyet al. 2008; Evans et al. 2018), wavelength range, andlikely high altitude clouds (see discussion in section 6)we do not expect to detect any of those 7 species. As aprecaution, we did run Exoretrievals with the 3 casesdiscussed above (featureless, activity, scatters) and in-cluding each of the 7 species individually. As expected,there was insufficient evidence to support any of thosemodels.The best model fit from
Exoretrievals was a flatspectrum. However, the | ∆ ln Z | between the featurelessmodel and the majority of the others was below 2.5, im-plying that those models cannot be distinguished from aflat spectrum. For PLATON the best fit was a model witha near-Rayleigh scattering slope (5.3 +2 . − . ), log metallic-ity of 1.3 +1 . − . , and C/O = 0.52 +0 . − . , but again the sup-port for this model is indistinguishable from the others( ∆ ln Z = 1.7). The undefined features (i.e lack of Naor K) and nominal detection of a scattering slope foundhere is consistent with the results of Gibson et al. (2017).Gibson et al. found no significant detection of K or Naand had only a tentative detection of a scattering slopewith their FORS2 data. It was only when including thebluer HST/STIS data between 0.3 and 1.0 µ m that they CCESS: WASP-31b Table 4.
The priors for
Exoretrievals and
PLATON . These priors were set to allow for a wide parameter space to be surveyed,but contained within physical regimes. Not all parameters were included in each model fit (see Tables 5). We used 10000 livepoints for each
Exoretrievals run and 5000 for
PLATON . For further description of the parameters of
Exoretrievals , pleaserefer to the Appendix of Espinoza et al. (2019).
Exoretrievals PLATON parameter prior function prior bounds parameter prior function prior bounds cloud top log-uniform -8 to 2 cloud top log-uniform -3 to 7pressure (P , bars) pressure (P , pascals)planetary atmospheric uniform 800 to planetary atmospheric uniform 800 totemperature (T p ) 1900K temperature (T p ) 1900Kstellar temperature) uniform 6000.0 to stellar temperature gaussian mean=6300K,(T occ ) 6300.0K (T star ) std=150Kstellar heterogeneities uniform 2300.0 to stellar heterogeneities uniform 5900.0 totemperature (T het ) 7000.0K temperature (T spot ) 6700.0Kspot covering uniform 0 to .8 spot covering uniform 0 to .8fraction (f het ) fraction (f spot )offset (depth) gaussian mean=0, offset (depth) uniform -8400 tostd=1000ppm 8400ppmcloud/haze amplitude (a) uniform -30 to 30 scattering factor log-uniform -10 to 10cloud/haze power law ( γ ) a uniform -10 to 4 scattering slope ( α ) b uniform -4 to 10cloud absorbing log-uniform -80 to 80 metallicity (Z/Z odot ) log-uniform -1 to 3cross-section ( σ cloud )trace molecules’ log-uniform -30 to 0 C/O uniform 0.05 to 2mixing ratiosreference radius factor ( f ) uniform 0.8 to 1.2 reference radius (R p ) uniform 1.35 to 1.6R j a This is the exponent of the scattering slope power law, where − is a Rayleigh scattering slope. b This is the wavelength dependence of scattering, with 4 being Rayleigh. confidently detect a scattering slope. Unfortunately ourMagellan data does not reach that far into the blue.Though the highest model with
Exoretrievals wasthat of a featureless atmosphere, in order to com-pare results from both retrievals, we examine the
Exoretrievals scattering model with the highest evi-dence ( ∆ ln Z = -1.77). This model is still indistinguish-able from the featureless model. We find that all over-lapping parameters agree with each other within theiruncertainties. The planet’s terminator temperature ofT p = 1223 +429 − K with
Exoretrievals agrees with thecalculated equilibrium temperature (1575K Andersonet al. 2011), and the T p of 1219 +262 − K retrieved with
PLATON is nearly 1- σ away from the calculated equilib-rium temperature. A low retrieved terminator temper-ature is not uncommon, and an inherent bias producedwhen using 1-D transmission spectrum models for a 2-Dproblem (Pluriel et al. 2020; MacDonald et al. 2020).However, using equation 1 of López-Morales & Seager(2007) to calculate the equilibrium temperature of the planet, we find the retrieved temperatures are still con-sistent if the planet has a bond albedo above 0.3 . Thismight be a slightly high albedo for a hot Jupiter butis completely reasonable relative to the solar system’sjovians (Rowe et al. 2008; Mallonn et al. 2019).The log cloud top pressure of -2.9 +3 . − . bars with Exoretrievals and -3.6 +2 . − . bars with PLATON , whilenot well constrained, implies that there are high-altitudeclouds. High-altitude clouds provide a good explanationfor the data, given that features such as Na and K arenot present, the possible water feature in the IR datais strongly muted, and the retrieved parameters havelarge uncertainties. But again, it cannot be confidentlyclaimed that this is the case because the lack of con-vergence of the cloud top pressure. The γ of -3.8 +3 . − . obtained with Exoretrievals and scattering slope of5.3 +2 . − . with PLATON are within 1 σ of a Rayleigh scatter- assuming efficient heat distribution, f =1/4, and using theretrieved system parameters McGruder et al. ing slope, which was detected at some level by previousstudies (Sing et al. 2015; Gibson et al. 2017).The highest ∆ ln Z transmission spectrum model re-trieved is that of a scattering slope with PLATON , andis overlayed on top of our combined Magellan/IMACSdata in Figure 7, with corner plot results shown in Fig-ure 12.Previous detailed atmospheric studies on WASP-31b’satmosphere with data introduced in Sing et al. (2015)are consistent with our retrieval analysis using IMACSdata introduced in this manuscript. Specifically, Wake-ford et al. (2017) fit models to the HST/STIS WASP-31bdata and found that clouds made of enstatite (MgSiO )and iron act as a gray opacity source from log pres-sures of -4.1 to -1.7 in WASP-31. This agrees with thecloud deck pressure, and muted features we retrieved. Infact, given the high uncertainties of our retrieved scat-tering slope, a near gray (little color dependent scatter-ing) cloud deck is also supported. Barstow et al. (2016)fit 3600 atmospheric models to the HST/STIS WASP-31 data, and found a best fit for a gray cloud deck at100mbars, but there were a few high evidence modelswith a Rayleigh slope and cloud top pressure as low as0.01mbars. Pinhas & Madhusudhan (2017) also foundmultiple indistinguishable fits for the HST/STIS WASP-31b data with their models, but had a best fit slope of-5.52+/-1.27 which agrees with our best fits. These find-ings in conjunction with ours emphasizes the difficultyin converging on the exact picture of WASP-31b’s at-mosphere, due to cloud formation.5.2. Previous Data
In order to test the consistency of our IMACS ob-servations and previous observations of WASP-31b, weran the same retrieval analysis against the three sepa-rate optical data in the literature, which were the orig-inal HST/STIS (Sing et al. 2015), VLT/FORS2, andthe combined FORS2 and reanalyzed STIS data (Gib-son et al. 2017). We combined each of these subsets ofdata with the WFC3 and
Spitzer
IR data.In summary, the VLT/FORS2 retrieval analysis pro-duced nearly identical results to our IMACS analysis.That is, inconclusive support for a flat spectrum with
Exoretrievals and a slightly favored, indistinguishablefrom a flat line ( ∆ ln Z = 0.26), near-Rayleigh opticalscattering slope with Platon . Given the similar wave-length coverage of IMACS (0.4–0.96 µ m) and FORS2(0.39–0.85 µ m) and their precision (0.0033 R p /R s ), weexpect to see this agreement.We detected potassium ( ∆ ln Z =3.8) and a scatteringslope with Exoretrievals when performing retrievalanalysis with the HST/STIS data. Interestingly, the best model with
Platon was one with strong support( ∆ ln Z =20.9) towards only stellar activity, which islikely a portrayal of a degeneracy between scatteringslope and stellar activity (Rackham et al. 2017).The evidences obtained with Exoretrievals , utilizingthe combined FORS2 and reanalyzed STIS data subset,favored the model with water and a scattering slope at ∆ ln Z =2.5. This data subset also moderately favored ascattering slope with Platon at ∆ ln Z =4.2.With both the STIS and combined FORS2 and reana-lyzed STIS data subsets, the retrieved scattering slopeswere higher (> 2 σ ) than Rayleigh. This is likely dueto the fact that these data subsets include bluer (0.3–0.4 µ m) wavelengths, which constrain the high scatteringslope. We tested this by running the combined FORS2and reanalyzed STIS data excluding the first two data-points (centered at 3350 and 3859Å) against Platon ,and obtained a shallower slope that was consistent withRayleigh scattering at 1 σ . RETRIEVAL INTERPRETATION6.1.
Optical Scattering Slope
As discussed in Sections 5.1 and 5.2, the Bayesian ev-idences of the FORS2 or IMACS datasets were unableto distinguish a flat spectrum from that with a scatter-ing slope. The models that fit for a scattering slope withthese data did retrieve slopes that were within Rayleigh,but those retrieved slopes had large uncertainties. Incontrast, retrieval runs utilizing the STIS or combinedFORS2 and reanalyzed STIS data subsets moderatelydetected a scattering slope. These slopes were aboveRayleigh (| γ | = α = 8.3) and better constrained withapproximate error bars of ± µ m.6.2. Sodium & Potassium Features
Of the four available WASP-31b optical datasets, onlythe STIS data had a Bayesian evidence that favoredpotassium absorption. As Gibson et al. (2017) imply, Platon is only capable of calculating the spectrum from 0.3 to30.0 µ m. Therefore, we had to shift the central wavelength of thebluest datapoint by 50Å. To do so, we did a linear interpolationbetween the three bluest bins, using the linear method of the scipy.interpolate package (Virtanen et al. 2020). We found a45ppm change in depth. CCESS: WASP-31b Table 5. ∆ ln Z for various Exoretrievals (left) and
PLATON (right) models relative to a featureless spectrum with the subsetof data that included the combined Magellan/IMACS (excluding transit UT130226), HST/WFC3, and
Spitzer
IR data. For
Exoretrievals no model is significantly preferred, but a flat spectrum has the highest evidence. For
PLATON a spectrum withan optical scattering slope is slightly favored, but indistinguishable from the other models.
Exoretrievals PLATON
Model: featureless H O Na K H O + K + Na Model: featureless . − . − . − . − . featureless . scatterers − − − − . − . − . − . scattering . activity − . − . − . − . − . stellar activity . scatterers & activity − − − − . − . − . − . Both . Figure 7.
Combined transmission spectrum of the three usable Magellan/IMACS optical transits, HST/WFC3, and
Spitzer
IRtransits. With this data the best-fit model (highest ∆ ln Z ) was one utilizing Platon and obtained a scattering slope consistentwith Rayleigh (5.3 +2 . − . ), log metallicity (Z/Z (cid:12) ) of 1.3 +1 . − . , C/O ratio of 0.52 +0 . − . , planet limb temperature of 1220 +260 − K, andlog cloud top pressure of -3.6 +2 . − . bars. The solid black line is the retrieval’s best fit (highest evidence), the shaded region is the1- σ confidence interval, and the colored dots with error bars are the data for various instruments. McGruder et al. this is likely produced as a side effect of the detrendingtechnique used. Another cause may be that the starwas active during the HST observations, which Pontet al. (2013), Rackham et al. (2017), and Espinoza et al.(2019) have shown could mimic transmission spectralsignals. As introduced in Section 5.2, the Bayesian ev-idence using
Platon strongly favored a model charac-terizing the STIS transmission spectrum with stellar ac-tivity. This is supported by the the AIT photometricmonitoring that only showed signs of activity during the2012-2013 season, which coincides with the STIS obser-vations taken in June of 2012. However, Keck/HIRESobservations of WASP-31’s Ca II H & K emission linessuggest that the star was extremely quiet (logR’ HK =-5.225) when observed (2002-2009) (Sing et al. 2016).Thus, it is possible, while not certain, that stellar activ-ity could contribute to the cause of this unique trans-mission spectrum. No matter the cause of the potas-sium signal in the original STIS data, given that thereare two independent low-resolution transmission spectraand a high-resolution transmission spectrum that showno signs of K, we conclude that WASP-31b shows nooptical atomic absorption.It could also be argued that because the K feature(7681Å) is near prominent tellurics ( ∼ ∼ Water Features
None of the fits in Sections 5.1 and 5.2, aside from theFORS2 and reanalyzed STIS data subset, showed signif-icant detection of H O. However, in general the modelsthat include H O are more favored (see Table 5), sug-gesting that there are signs of water features in the IRdata. This is consistent with the findings of Sing et al.(2015) and Stevenson (2016), where Stevenson (2016)found a weak H O – J index (0.86 ± +0 . − . , which is within uncertainties of the re- trieved C/O ratios using the other data subsets, suggestan oxygen-rich atmosphere. The C/O ratio is relativelyunconstrained, because the only observations that couldconstrain CO features are two, wide wavelength datapoints from Spitzer . None-the-less this ratio and therelatively low temperature retrieved (T p = 1220 +260 − K)imply that water would be a large absorber in WASP-31b’s atmosphere (Madhusudhan 2012). Thus, the lackof significant detection of water is likely due to high-altitude clouds muting the water features.6.4.
Final Interpretation
The possible extremely-muted water features, lack ofoptical absorption features from K and Na, and retrievedlow pressures with all data subsets (see Sections 5.1 &5.2) point towards WASP-31b possessing high-altitudeclouds. Furthermore, comparing our spectra and re-trieval results with atmospheric modeling by Wakeford& Sing (2015) and Pinhas & Madhusudhan (2017) thedata suggests: 1) Their likely is a large number of par-ticles of modal size between ∼ × − and .25 µ m loftedhigh in WASP-31b’s atmosphere. 2) The compositionof these particles are predominately enstatite and othersilicates. 3) These high opacity structures are clouds,which is an accumulation of particles that condense inequilibrium or near equilibrium chemistry; rather thanhazes, which form from photochemistry or other non-equilibrium chemical processes. Given that this planetlikely has clouds obscuring the optical and near-IR fea-tures, spectral observations in the IR (6-30 µ m) couldclarify the cloud composition and particle size distribu-tion. As suggested by Wakeford & Sing (2015) and Pin-has & Madhusudhan (2017), JWST observations wouldbe optimal for such clarity. SUMMARY & CONCLUSIONWe observed four transit epochs of WASP-31b be-tween 2013 and 2019 with the Magellan/IMACS spec-trograph. From these observations we derived opti-cal transmission spectra of the planet between 0.4 µ mand 0.97 µ m. We excluded the first transit (UT130226),which had inconsistent retrieved parameters, and com-bined the other three transits to produce a transmissionspectrum with an average precision of 0.0033 R p /R s per187Å bin widths. We combined our transmission spec-trum with literature HST/WFC3 and Spitzer observa-tions and modelled the full datasets using two indepen-dent retrieval codes. We find signs of a Rayleigh scat-tering slope and water features, but neither are signifi-cantly detected. We also find no Na or K signals likelydue to a high-altitude cloud deck. Our
PLATON scatter-ing slope model had the highest ∆ ln Z , and retrieved CCESS: WASP-31b
17a Rayleigh scattering slope of 5.3 +2 . − . , log metallicity of1.3 +1 . − . , C/O of 0.52 +0 . − . , and cloud top pressure of -3.6 +2 . − . bars.In addition, we ran the retrievals against data from theliterature (Sing et al. 2015; Gibson et al. 2017), whichall showed weak signs of water absorption likely mutedby high-altitude clouds. The retrieved scattering slopesfor the datasets including the bluer STIS observationswere steeper than Rayleigh by over 2 σ . This could meanthat WASP-31b’s scattering slopes is beyond Rayleigh,but we acknowledge there could be other causes of thisretrieved slope. Of all available data on WASP-31b, onlythe original STIS transit showed strong K absorption.From this work we find that the atmosphere of WASP-31b is likely to be significantly covered by high-altitudeclouds, explaining high uncertainties of parameters andlow significance of feature detection.It is imperative to understand which planets formhigh-altitude clouds or hazes, and how. Understandingthis will be key in the JWST era. We need a large li-brary of thoroughly analyzed optical to IR transmissionspectra for this. In particular, the optical componentallows for absolute abundances of molecules retrieved inthe IR by constraining the level at which the molecularfeatures are muted by clouds. Such a library would al-low for better comparative studies on observed spectralfeatures and other system parameters (e.g., stellar irra-diation levels, spectral type, mass, radius, planet densityetc.) to find correlations, and will drastically enhancethe field of exoplanetology. Uniform datasets for a widerange of planets, such as the ones ACCESS is building,will be crucial for such analyses.We appreciate the anonymous referee for helpful com-ments and feedback. The results reported herein ben- efited from support, collaborations and information ex-change within NASA’s Nexus for Exoplanet System Sci-ence (NExSS), a research coordination network spon-sored by NASA’s Science Mission Directorate. This pa-per includes data gathered with the 6.5 meter Magel-lan Telescopes located at Las Campanas Observatory,Chile. We thank the staff at the Magellan Telescopesand Las Campanas Observatory for their ongoing inputand support to make the ACCESS observations pre-sented in this work possible. We also appreciate thesupport from the NSF Graduate Research Fellowship(GRFP), grant number DGE1745303. B.V.R. thanksthe Heising-Simons Foundation for support. G.W.H. ac-knowledges long-term support from NASA, NSF, Ten-nessee State University, and the State of Tennesseethrough its Centers of Excellence program. A.J. ac-knowledges support from FONDECYT project 1171208,and from ANID âĂŞ Millennium Science Initiative âĂŞICN12_009. C.D.M. appreciates the help Josh Speagle,a colleague and friend, provided in better understand-ing the Bayesian statistics behind the detrending andsampling techniques used. Software:
Astropy (Astropy Collaboration et al.2013), corner (Foreman-Mackey 2016), Matplotlib(Hunter 2007), NumPy (Oliphant 2006), Multinest(Feroz et al. 2009), PyMultiNest (Buchner et al. 2014),SciPy (Jones et al. 2001), batman (Kreidberg 2015),george (Foreman-Mackey et al. 2014) dynesty (Speagle2019), PLATON (Zhang et al. 2019), ld-exosim (Es-pinoza & Jordán 2016)
Facilities:
Magellan:Baade, Smithsonian InstitutionHigh Performance Cluster (SI/HPC)REFERENCES
Aigrain, S., Pont, F., & Zucker, S. 2012, MNRAS, 419, 3147Anderson, D. R., Collier Cameron, A., Hellier, C., et al.2011, A&A, 531, A60Astropy Collaboration, Robitaille, T. P., Tollerud, E. J.,et al. 2013, A&A, 558, A33Barstow, J. K., Aigrain, S., Irwin, P. G. J., & Sing, D. K.2016, The Astrophysical Journal, 834, 50Benneke, B., & Seager, S. 2013, ApJ, 778, 153Berta, Z. K., Charbonneau, D., Bean, J., et al. 2011, ApJ,736, 12Bixel, A., Rackham, B. V., Apai, D., et al. 2019, AJ, 157, 68Buchner, J., Georgakakis, A., Nandra, K., et al. 2014,A&A, 564, A125 Charbonneau, D., Berta, Z. K., Irwin, J., et al. 2009,Nature, 462, 891Chen, G., Pallé, E., Welbanks, L., et al. 2018, A&A, 616,A145Colón, K. D., Ford, E. B., Redfield, S., et al. 2012,MNRAS, 419, 2233Diamond-Lowe, H., Berta-Thompson, Z., Charbonneau, D.,& Kempton, E. M. R. 2018, AJ, 156, 42Dressler, A., Bigelow, B., Hare, T., et al. 2011, PASP, 123,288Duvenaud, D. 2014, PhD thesis, University of CambridgeEspinoza, N. 2017, PhD thesis, Pontificia UniversidadCatólica de Chile McGruder et al.
Espinoza, N., & Jordán, A. 2016, MNRAS, 457, 3573Espinoza, N., Rackham, B. V., Jordán, A., et al. 2019,MNRAS, 482, 2065Evans, T. M., Sing, D. K., Goyal, J. M., et al. 2018, AJ,156, 283Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS,398, 1601Foreman-Mackey, D. 2016, The Journal of Open SourceSoftware, 24Foreman-Mackey, D., Hoyer, S., Bernhard, J., & Angus, R.2014, George: George (V0.2.0)Fortney, J. J., Lodders, K., Marley, M. S., & Freedman,R. S. 2008, ApJ, 678, 1419Gao, P., Thorngren, D. P., Lee, G. K. H., et al. 2020,Nature Astronomy, arXiv:2005.11939 [astro-ph.EP]
Gibson, N. P. 2014, Monthly Notices of the RoyalAstronomical Society, 445, 3401Gibson, N. P., Aigrain, S., Roberts, S., et al. 2012,MNRAS, 419, 2683Gibson, N. P., de Mooij, E. J. W., Evans, T. M., et al.2019, MNRAS, 482, 606Gibson, N. P., Nikolov, N., Sing, D. K., et al. 2017,MNRAS, 467, 4591Gillon, M., Jehin, E., Lederer, S. M., et al. 2016, Nature,533, 221Hastie, T., Tibshirani, R., & Friedman, J. 2007, TheElements of Statistical Learning, 2nd edn.(Springer-Verlag)Heng, K., & Kitzmann, D. 2017, MNRAS, 470, 2972Henry, G. W. 1999, PASP, 111, 845Hunter, J. D. 2007, Computing In Science & Engineering,9, 90Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013,A&A, 553, A6Jolliffe, I. 2002, Principal Component Analysis (New York:Springer-Verlag)Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy:Open source scientific tools for PythonJordán, A., Espinoza, N., Rabus, M., et al. 2013, ApJ, 778,184Kirk, J., López-Morales, M., Wheatley, P. J., et al. 2019,AJ, 158, 144Kochanek, C. S., Shappee, B. J., Stanek, K. Z., et al. 2017,PASP, 129, 104502Kreidberg, L. 2015, PASP, 127, 1161Lomb, N. R. 1976, Ap&SS, 39, 447López-Morales, M., & Seager, S. 2007, ApJL, 667, L191López-Morales, M., Haywood, R. D., Coughlin, J. L., et al.2016, AJ, 152, 204 Luger, R., Sestovic, M., Kruse, E., et al. 2017, NatureAstronomy, 1, 0129Lupu, R. E., Zahnle, K., Marley, M. S., et al. 2014, ApJ,784, 27MacDonald, R. J., Goyal, J. M., & Lewis, N. K. 2020,ApJL, 893, L43MacDonald, R. J., & Madhusudhan, N. 2017, MonthlyNotices of the Royal Astronomical Society, 469, 1979MacDonald, R. J., & Madhusudhan, N. 2017, MNRAS, 469,1979Madhusudhan, N. 2012, ApJ, 758, 36Mallonn, M., Köhler, J., Alexoudi, X., et al. 2019, A&A,624, A62Mallonn, M., Herrero, E., Juvan, I. G., et al. 2018, A&A,614, A35Markevich, N., & Gertner, I. 1989, Nuclear Instrumentsand Methods in Physics Research A, 283, 72McCullough, P. R., Crouzet, N., Deming, D., &Madhusudhan, N. 2014, ApJ, 791, 55Micchelli, C. A., Xu, Y., & Zhang, H. 2006, J. Mach. Learn.Res., 7, 2651Morley, C. V., Fortney, J. J., Marley, M. S., et al. 2012,The Astrophysical Journal, 756, 172Oliphant, T. E. 2006, Guide to NumPy, Provo, UTPinhas, A., & Madhusudhan, N. 2017, MNRAS, 471, 4355Pluriel, W., Zingales, T., Leconte, J., & Parmentier, V.2020, A&A, 636, A66Pont, F., Sing, D. K., Gibson, N. P., et al. 2013, MNRAS,432, 2917Rackham, B., Espinoza, N., Apai, D., et al. 2017, ApJ, 834,151Rackham, B. V., Apai, D., & Giampapa, M. S. 2018, ApJ,853, 122—. 2019, AJ, 157, 96Rasmussen, C., & Williams, C. 2005, Gaussian Processesfor Machine Learning (Massachusetts Institute ofTechnology: Adaptive Computation and MachineLearning)Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, inProc. SPIE, Vol. 9143, Space Telescopes andInstrumentation 2014: Optical, Infrared, and MillimeterWave, 914320Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010,JQSRT, 111, 2139Rowe, J. F., Matthews, J. M., Seager, S., et al. 2008, ApJ,689, 1345Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015,PhyS, 90, 054005Scargle, J. D. 1982, ApJ, 263, 835Seager, S., & Sasselov, D. D. 2000, ApJ, 537, 916
CCESS: WASP-31b Sedaghati, E., Boffin, H. M. J., MacDonald, R. J., et al.2017a, Nature, 549, 238—. 2017b, Nature, 549, 238Sing, D. K., Désert, J. M., Fortney, J. J., et al. 2011, A&A,527, A73Sing, D. K., Wakeford, H. R., Showman, A. P., et al. 2015,MNRAS, 446, 2428Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature,529, 59Southworth, J., Mancini, L., Madhusudhan, N., et al. 2017,AJ, 153, 191Speagle, J. S. 2019, arXiv e-prints, arXiv:1904.02180Stevenson, K. B. 2016, ApJL, 817, L16Tennyson, J., Yurchenko, S. N., Al-Refaie, A. F., et al.2016, Journal of Molecular Spectroscopy, 327, 73 Thomas-Osip, J. E., Prieto, G., Johns, M., & Phillips,M. M. 2008, in Society of Photo-Optical InstrumentationEngineers (SPIE) Conference Series, Vol. 7012,Ground-based and Airborne Telescopes II, 70121UTrotta, R. 2008, Contemporary Physics, 49, 71Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020,Nature Methods, 17, 261Wakeford, H. R., & Sing, D. K. 2015, A&A, 573, A122Wakeford, H. R., Stevenson, K. B., Lewis, N. K., et al.2017, ApJ, 835, L12Weaver, I. C., López-Morales, M., Espinoza, N., et al. 2020,AJ, 159, 13Yan, F., Espinoza, N., Molaverdikhani, K., et al. 2020,arXiv e-prints, arXiv:2007.15485Yurchenko, S. N., Tennyson, J., Barber, R. J., & Thiel, W.2013, Journal of Molecular Spectroscopy, 291, 69Zacharias, N., Finch, C. T., Girard, T. M., et al. 2013, AJ,145, 44Zhang, M., Chachan, Y., Kempton, E. M. R., & Knutson,H. A. 2019, PASP, 131, 034501 McGruder et al.
APPENDIX A. LIGHT CURVES
CCESS: WASP-31b Table 6.
Magellan/IMACS optical transmission spectrum (planet radius/star radius, p ) for each of our four transit epochs:UT130226, UT130425, UT140222, and UT190314, respectively. The data were produced implementing the reduction anddetrending processes discussed in Section 3. These depths do not include the offsets for combining each spectra. Gaps inthe spectra (see Figure 1, prevented a few bins from exactly overlapping in wavelength space. Those bins were still weightedaveraged together (x and y direction), with the resulting bin composing of the full wavelength width of all bins. In this table,each bin for the combined transmission spectrum is separated by a vertical line. Wavelength UT130226 UT130425 UT140222 UT190314range ( Å ) R p /R s R p /R s R p /R s R p /R s . − . . +0 . − . . ± . . +0 . − . − − −− . − . . +0 . − . . +0 . − . . +0 . − . − − −− . − . . +0 . − . . +0 . − . . +0 . − . − − −− . − . . +0 . − . . +0 . − . . +0 . − . − − −− . − . . +0 . − . . +0 . − . − − −− − − −− . − . − − −− − − −− . +0 . − . − − −− . − . . +0 . − . . +0 . − . − − −− − − −− . − . . +0 . − . . ± . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . ± . . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − . . ± . . +0 . − . . +0 . − . . ± . . − . . +0 . − . . +0 . − . . +0 . − . . ± . . − . . +0 . − . . ± . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . − − −− − − −− . − . − − −− − − −− − − −− . +0 . − . . − . − − −− − − −− . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . − − −− − − −− . − . . +0 . − . . +0 . − . . ± . . +0 . − . . − . . +0 . − . . ± . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . ± . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − .
25 0 . +0 . − . . ± . . +0 . − . . +0 . − . . − .
25 0 . +0 . − . . ± . . +0 . − . . +0 . − . . − . − − −− − − −− . +0 . − . − − −− . − . . +0 . − . . ± . − − −− . +0 . − . . − . − − −− − − −− − − −− . +0 . − . . − . − − −− − − −− . +0 . − . − − −− . − . . +0 . − . . +0 . − . . +0 . − . . +0 . − . . − . . +0 . − . . +0 . − . . +0 . − . . ± . . − . . +0 . − . . +0 . − . . +0 . − . − − −− . − . − − −− − − −− − − −− . +0 . − . McGruder et al.
Figure 8.
Detrended, binned light curves of transit UT130226 (left) with residuals of the detrended data from the best-fittransit model on the right. The first axis on the left panel displays the mean wavelength of that bin in microns, and the secondaxis displays the light curves’ normalized flux, which were incrementally offset by 0.01 for visual purposes. The axis on theright panel displays the standard deviation ( σ ) of the residuals in ppm. Each bin was detrended with model averaging of PCAand GP machine learning techniques. These light curves encompass the spectrum from 4018Å to 9450Å, aside from the spacescaused by chip gaps. CCESS: WASP-31b Figure 9.
Same as Figure 8 but for transit UT130425. McGruder et al.
Figure 10.
Same as Figure 8 but for transit UT140222.
CCESS: WASP-31b Figure 11.
Same as Figure 8 but for transit UT190314. Because of the IMACS chip gaps, this particular dataset fullyencompasses the spectrum from 4918Å to 9670Å. McGruder et al.
Table 7.
Combined transmission spectrum from transit epochs UT130425, UT140222, and UT190314. The spectra werecombined following the procedure outlined in Section 4.1. In the retrieval analysis, because the retrievals do not take asymmetricwavelength errors, the wavelength range of the few overlapping bins were re-centered based on the weighted mean wavelengthof the bin.
Wavelength range (Å) Mean (Å) R p /R s . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
50 0 . +0 . − . . − . .
65 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
30 0 . +0 . − . . − . .
90 0 . +0 . − . . − . .
45 0 . +0 . − . . − . .
78 0 . +0 . − . . − . .
62 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − .
25 7553 .
62 0 . +0 . − . . − .
25 7681 .
25 0 . +0 . − . . − . .
01 0 . +0 . − . . − . .
69 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
00 0 . +0 . − . . − . .
70 0 . +0 . − . CCESS: WASP-31b B. TRANSIT UT130226As can be seen in Figure 5, the transit depths of epoch UT130226’s transmission spectrum varied widely over ourwavelength coverage. Additionally, as can be seen in Figure 4, there was little baseline before ingress during thisobservation, which could cause improper detrending of systematics. We tested the impact of our insufficient baselineby comparing the transit duration of epoch UT130226 with those of other epochs when all transit parameters, asidefrom p , were fixed. We also compared the transit depth obtained for this epoch between this test case (with fixedparameters) and our nominal case (with other transit parameters free). For both cases, the transit durations and depthsagreed within 1 σ , suggesting that the baseline was sufficient. This is not conclusive support of sufficient baseline, butreasoning to explore if the transmission spectrum is indeed caused by the system.We used our retrieval framework to analyze just IMACS transit UT130226 and the IR data to further explore thecause of the outlier transmission spectrum . It is possible the transit can be explained by unique physical scenarios inthe system during the single epoch, such as the host star being relatively active or the planet having an exceptionallylarge amount of scattering agents lofted high in the atmosphere. We tested the plausibility of these scenarios byrunning a variety of retrievals. Table 8 shows ∆ ln Z for these retrievals compared to a featureless spectrum. Theevidences for both PLATON and
Exoretrievals strongly support stellar activity and a large scattering slope ( α =9.61 +0 . − . and γ = -9.50 +0 . − . ).The retrieved large scattering slopes are inconsistent with the other Magellan/IMACS and VLT/FORS2 observations(discussed in Sections 5). Additionally, the best-fit terminator temperatures of 1805 +66 − K for
Exoretrievals and1846 +35 − K for
PLATON were high relative to the other IMACS and FORS2 data ( ∼ Exoretrievals the best-fit stellar activity model was that where T het = 5720 ±
310 Kand f het = 0.4 +0 . − . , and for PLATON it was where T spot = 5950 +240 − K and f spot = 0.31 +0 . − . . These high spot coveringfractions imply a completely contrary level of activity from other activity indicators. Thus, we have no concreteexplanation of this transmission spectrum’s abnormal features, and decided to exclude this transit from our combineddata in order to prevent biasing our atmospheric interpretations. Table 8.
Same as Table 5 but with the subset of data that included only the UT130226 Magellan/IMACS and IR data. Theretrievals with a strong scattering slope ( γ ∼ − . , α ∼ . ) and stellar activity (T het ∼ f het ∼ spot ∼ f spot ∼ Exoretrievals and
PLATON evidences, respectively.
Exoretrievals PLATON
Model: featureless H O Na K H O + K + Na Model: featureless . − . − . − . − . featureless . scatterers − − − .
53 14 .
13 14 .
00 13 . scattering . activity .
08 10 .
04 3 .
02 3 .
40 9 . stellar activity . scatterers & activity − − − .
11 15 .
08 14 .
72 15 . Both . C. CORNER PLOTS The signal of one Magellan/IMACS transit (0.007 R p /R s ) is not enough to detect features like Na or K ( ∼ R p /R s , Sing et al.2015), but sufficient to interpret the large slope in the spectrum. McGruder et al.
Figure 12.
Corner plot of the best-fit retrieval on the combined Magellan/IMACS, HST/WFC3, and
Spitzer data. This isusing the