Ground-based Transmission Spectroscopy with VLT FORS2: Evidence for Faculae and Clouds in the Optical Spectrum of the Warm Saturn WASP-110b

We observed two primary transits of WASP-110b on UT 2017 June 30th and September 25th with FORS2 (Appenzeller et al. 1998) attached to the Unit Telescope 1 (UT1, Antu) of the VLT at the European Southern Observatory (ESO) on Cerro Paranal in Chile (Figure 1) as part of Large Program 199.C-0467 (PI: Nikolov). We used a similar observing setup and strategy to our earlier VLT studies (Nikolov et al. 2016, 2018; Gibson et al. 2017; Carter et al. 2020; Wilson et al. 2020).

Zoom In Zoom Out Reset image size

During the two transits, we collected data in the multiobject spectroscopy mode with a mask consisting of two broad slits centered on the target and on a reference star of similar brightness. The reference star (known as 2MASS 20234536-4403456) is a bright source in the FORS2 field-of-view Vmag = 12.999 ± 0.013 and is located at an angular separation of $2^{\prime} .9$ eastward from the target and a distance of 604 ± 14 pc from the Earth. The reference is of higher temperature (T_eff = 5912 ± 126 K) with magnitude differences (target minus reference) from the PPMXL and Gaia DR2 catalogs: Δ B = −0.36, Δ G = −0.30, Δ V = −0.24, Δ R = −0.35, Δ J = −0.46, Δ H = −0.52, and Δ K = −0.55 (Roeser et al. 2010; Gaia Collaboration et al. 2018). We used broad slits spanning 22'' along the dispersion and ∼120'' along the spatial (perpendicular) axis to minimize slit light losses due to seeing variations and guiding imperfections. We observed both transits with the same slit mask and the red detector (Massachusetts Institute of Technology), which is a mosaic of two chips. We positioned the instrument field of view such that each detector imaged one source. To improve the duty cycle, we made use of the fastest available readout mode (200 kHz, ∼30 s). During both nights, we ensured that the Longitudinal Atmospheric Dispersion Corrector is in its neutral (park) position, i.e., inactive.

In addition to the bright reference star, a fainter star is located at ∼5'' eastward from WASP-110. The star (known as 6678937482610261632 in the Gaia DR2 catalog) is an A-type dwarf (T_eff = 9617 ± 128 K, M ∼ 2.280 M_⊙, where M_⊙ is the mass of the Sun, R ∼ 1.890 R_⊙, where R_⊙ is the radius of the Sun, $\mathrm{log}g\sim 4.2$, and luminosity L ∼ 21 L_⊙, where L_⊙ is the luminosity of the Sun) of an apparent brightness of Vmag = 14.992 ± 0.046 and Gmag = 14.9710 ± 0.0006 at a distance of 4640 ± 1237 pc, and is not physically associated with the WASP-110 system (Gaia Collaboration et al. 2018). The A-type dwarf is resolved from WASP-110 in the acquisition images and the spectroscopy time series of both nights.

During the first night we monitored the flux of WASP-110 and the comparison star under thin cirrus conditions. We used the dispersive element GRIS600B (hereafter blue), which covers the spectral range from 3600–6200 Å at a resolving power R = λ/Δλ ≈ 600. The field-of-view rose from an airmass 1.47 to 1.06 and at the end of the observations the field of view was at an airmass of 1.12. The seeing gradually improved from 05 to 04, as measured from a Gaussian fit to the spatial profile of the stellar spectra. We collected a total of 87 exposures for 5h with integration times of 80 and 200 s, for the first 11 and remaining 76 exposures.

During the second night, we collected data under photometric conditions. We made use of the dispersive element GRIS600RI (hereafter red), which covers the range from 5400–8200 Å, in combination with the GG435 blocking filter to isolate the first order. The field of view rose from an airmass of 1.063 to 1.060 and the observations ended with the target at airmass of 1.853. The seeing varied between 03 and 06. We monitored WASP-110 and the reference star for 4^h49^m with ∼8 minutes interruption at UT 3:33 due to a computer operational issue and collected a total of 154 spectra with integration times between 80 and 100 s.

WASP-110 was observed by the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) on Camera 1 CCD 2 during the Extended Mission between 2020 July 4 and 30 (Sector 27) with a 10 minutes cadence in the full frame images (FFIs). A light curve was produced using tools available in the Lightkurve v1.9 package (Lightkurve Collaboration et al. 2018) and applied to poststamps 15 × 15 pixel ($5\buildrel{\,\prime}\over{.} 25\times 5\buildrel{\,\prime}\over{.} 25$) wide. They were extracted from FFIs with the TESSCut ¹⁵ online tool (Brasseur et al. 2019) with WASP-110 centered in a frame. Fluxes were obtained employing a quadratic aperture which was three pixels wide without one corner pixel in order to avoid blending with a nearby bright star. The sky background level was determined with an algorithm based on standard-deviation thresholding. Trends caused by systematic effects or stellar variation were removed with the Savitzky–Golay filter with a window width of 12 hr. Prior to this step, in-transit and in-occultation data points were masked out using a trial ephemeris.

The final light curve is plotted in the upper panel in Figure 2. Six consecutive transits of WASP-110 b were observed. For further analysis, the individual transit light curves were extracted with time margins equal to ±2.5 times a transit duration. Their photometric noise rates (Fulton et al. 2011) were found to be between 4.0 and 6.0 parts per thousand (ppth) of the normalized flux per minute of observation. Individual light curves are displayed in the middle and bottom panels in Figure 2.

Zoom In Zoom Out Reset image size

We divided the raw flux of the target by the raw flux of the reference star to remove the variations of Earth's atmospheric transparency, Figure 1. We modeled the transit and systematic effects of the white-light curves by treating the data as a Gaussian process ¹⁶ and assuming quadratic limb darkening for the star (Gibson et al. 2012). The transit parameters: midtime T_mid, orbital inclination i, normalized semimajor axis a/R_* (where R_* is the radius of the star), the planet-to-star radius ratio R_p /R_*, and the two limb-darkening coefficients u₁ and u₂ were allowed to vary in the fit to each of the two white-light curves, while the orbital period was held fixed to the previously determined value of Anderson et al. (2014).

Under the Gaussian process assumption, the data likelihood is a multivariate normal distribution with a mean function μ describing the deterministic transit signal and a covariance matrix K that accounts for stochastic correlations (i.e., poorly constrained systematics) in the light curves:

$$\begin{eqnarray}&&p({\boldsymbol{f}}| {\boldsymbol{\theta }},\gamma )={ \mathcal N }(\mu ,K),\end{eqnarray} \tag{ 1 }$$

where p is the probability density function, f and θ are vectors containing the flux measurements and mean function parameters, respectively, γ is a function containing the covariance parameters, and ${ \mathcal N }$ is a multivariate normal distribution. The mean function μ is defined as:

$$\begin{eqnarray}&&\mu ({\boldsymbol{t}},\hat{{\boldsymbol{t}}};{c}_{0},{c}_{1},{\boldsymbol{\theta }})=[{c}_{0}+{c}_{1}\hat{{\boldsymbol{t}}}]\,T({\boldsymbol{t}},{\boldsymbol{\theta }}),\end{eqnarray} \tag{ 2 }$$

where t is a vector of all central exposure timestamps in Julian Date, $\hat{{\boldsymbol{t}}}$ is a vector containing all standardized times (that is, with subtracted mean exposure time and divided by the standard deviation), c₀ and c₁ describe a linear baseline trend, T( θ ) is an analytical expression describing the transit, and θ = (i, a/R_*, T_mid, R_p /R_*, u₁, u₂). We made use of the analytical formulae of Mandel & Agol (2002) and Kreidberg (2015).

We computed the theoretical values of the coefficients of the quadratic limb darkening law using a three-dimensional (3D) stellar atmosphere model grid (Magic et al. 2015), factored by the throughputs of the blue and red grisms. In these calculations, we adopted the closest match to the effective temperature, surface gravity, and metallicity of the exoplanet host star found in Anderson et al. (2014). The choice of a quadratic versus a more complex law (such as a four-parameter nonlinear law) was motivated by the study of Espinoza & Jordán (2016), in which the two-parameter law has been shown to introduce negligible bias on the measured properties of transiting systems similar to WASP-110. The quadratic law also requires a much shorter computational time to compute the model transit light curve. The theoretical limb darkening coefficients were employed as priors in our light-curve fits.

Similar to our earlier VLT studies, we defined the covariance matrix as $K={\sigma }_{i}^{2}{\delta }_{{ij}}+{k}_{{ij}}$, where σ_i contains the photon noise uncertainties, δ_ij is the Kronecker delta function, and k_ij is a covariance function. The white noise term σ_w was assumed to have the same value for all data points and was allowed to freely vary. We chose to use the Matérn ν = 3/2 kernel with time t for the covariance function. Our kernel choice is motivated by the study of Gibson et al. (2013b), where the Matérn ν = 3/2 kernel is empirically motivated using simulated data, and is the first to use this kernel for light-curve analysis. We also experimented by adding additional terms, including the spectral dispersion and cross-dispersion drifts x and y as input variables, and the full width at half maximum (FWHM) measured from the cross-dispersion profiles of the two-dimensional spectra and the speed of the rotation angle z. As with the linear time term, we also standardized the input parameters before the light-curve fitting. We chose to use the time for both observations instead of combinations of other terms because a GP of time resulted in well-behaved residuals. The covariance function was defined as:

$$\begin{eqnarray}&&{k}_{{ij}}={A}^{2}(1+\sqrt{3}{D}_{{ij}})\exp -\sqrt{3}{D}_{{ij}},\end{eqnarray} \tag{ 3 }$$

where A is the characteristic correlation amplitude and

$$\begin{eqnarray}&&{D}_{{ij}}=\sqrt{\displaystyle \frac{{\left({\hat{{\boldsymbol{t}}}}_{i}-{\hat{{\boldsymbol{t}}}}_{j}\right)}^{2}}{{\tau }_{t}^{2}}},\end{eqnarray} \tag{ 4 }$$

where τ_t is the correlation length scale and the hatted variables are standardized. The parameters X = (c₀, c₁, T_mid, i, a/R_*, R_p /R_*, u₁, u₂) and Y = (A, τ_t ) were allowed to vary and fixed the orbital period P to its literature value. Uniform priors were adopted for X and log-uniform priors for Y .

We marginalized the posterior distribution p(θ, γ∣f) ∝ p(f∣θ, γ)p(θ, γ) using the Markov chain Monte Carlo software package emcee from Foreman-Mackey et al. (2013). We identified the maximum likelihood solution using the Levenberg–Marquardt least-squares algorithm (Markwardt 2009) and initialized three groups of 150 walkers close to that maximum. The first two groups were run for 350 samples and the third one for 4500 samples. To ensure faster convergence, we resampled the positions of the walkers in a narrow space around the position of the best walker from the first run before running for the second group. This helps prevent some of the walkers starting in a low-likelihood area of parameter space, which can require more computational time to converge. Figure 1 and Table 1 show transit models for each of the two observations computed using the marginalized posterior distributions and the fitted parameters, respectively. We find residual dispersion of 378 and 223 parts per million for the blue and red light curves, respectively.

Table 1. System Parameters

Parameter	Value
Orbital period (day)	3.7783977 (fixed)
Orbital eccentricity	0 (fixed)
GRIS600B, λ (Å)	4013−6173
Tmid (MJD)	${57934.23238}_{-0.00028}^{+0.00026}$
i (°)	${88.00}_{-0.37}^{+0.45}$
a/R*	${10.73}_{-0.35}^{+0.30}$
Rp /R*	${0.1390}_{-0.0064}^{+0.0063}$
u1	${0.52}_{-0.24}^{+0.22}$
u2	${0.61}_{-0.47}^{+0.50}$
$\mathrm{ln}A$	$-{13.0}_{-1.5}^{+3.3}$
$\mathrm{ln}{\eta }_{\mathrm{time}}$	$-{0.9}_{-1.5}^{+2.7}$
c0	${1.0000}_{-0.0018}^{+0.0011}$
c1	$-{0.00032}_{-0.0007}^{+0.0010}$
GRIS600RI, λ (Å)	5293−8333
Tmid (MJD)	${58021.13912}_{-0.00023}^{+0.00025}$
i (°)	${88.52}_{-0.46}^{+0.67}$
a/R*	${11.19}_{-0.20}^{+0.21}$
Rp /R*	${0.1389}_{-0.0044}^{+0.0041}$
u1	${0.18}_{-0.22}^{+0.25}$
u2	${0.74}_{-0.44}^{+0.40}$
$\mathrm{ln}A$	$-{12.4}_{-1.3}^{+3.3}$
$\mathrm{ln}{\eta }_{\mathrm{time}}$	$-{0.7}_{-1.2}^{+2.4}$
c0	${1.0014}_{-0.0027}^{+0.0013}$
c1	${0.00097}_{-0.00080}^{+0.0017}$
Weighted mean:
i (°)	${88.16}_{-0.32}^{+0.37}$
a/R*	${11.07}_{-0.17}^{+0.18}$
GRIS600B	(fixed i, a/R*, Tmid)
Rp /R*	${0.1435}_{-0.0036}^{+0.0047}$
u1	${0.66}_{-0.16}^{+0.14}$
u2	${0.19}_{-0.19}^{+0.20}$
GRIS600RI	(fixed i, a/R*, Tmid)
Rp /R*	${0.1396}_{-0.0037}^{+0.0033}$
u1	${0.12}_{-0.17}^{+0.20}$
u2	${0.68}_{-0.24}^{+0.20}$

Download table as:  ASCIITypeset image

We computed the weighted mean values of i and a/R_* and repeated the fits allowing only R_p /R_*, u₁ and u₂ to vary. We held fixed i and a/R_* to the weighted mean values and the two central times, T_mid, were held fixed to the values determined from the first fit. We report our results in Table 1.

The TESS light curves were used to refine the transit ephemeris for WASP-110b. In addition, we used a follow-up photometric time series from the discovery paper (Anderson et al. 2014). This early epoch observation was acquired on 2013 August 16 with the 1.2 m Euler-Swiss telescope equipped in the EulerCam photometer. In total, 121 measurements were collected in a Cousins I filter with a median cadence of 137 s. The photometric noise rate was found to be 1.15 ppth per minute of observation. The original timestamps in HJD_UTC were converted into BJD_TDB using an online applet ¹⁷ (Eastman et al. 2010). For timing purposes, we do not correct the TESS transit light curves for the third-light contamination. Thus, we used fluxes prior the dilution correction applied in Anderson et al. (2014).

The midtransit times were determined by fitting a trial transit model with the Transit Analysis Package (TAP; Gazak et al. 2012). The TESS and EulerCam light curves were modeled simultaneously with the orbital inclination, the semimajor axis scaled in star radii, and the ratio of planet-to-star radii linked together. The limb darkening effect was approximated with a quadratic law, the coefficients of which were bilinearly extrapolated from tables of Claret & Bloemen (2011). The coefficients for the TESS passband were calculated as averages of R, I, and Sloan Digital Sky Survey $z^{\prime}$ values. Their values were allowed to vary under the Gaussian penalties with a conservative width of 0.1. Possible trends in the time domain were accounted for using second order polynomials for each transit light curve. The best-fitting values and their uncertainties were calculated as the medians and 15.9 and 84.1 percentiles of the marginalized posteriori probability distributions generated from 10 MCMC chains, each 10⁶ steps long with 10% burn-in phase. The results are given in Table 2.

We converted to BJDs the central transit times from our VLT observations and combined these with the BJD TDB central times from our Euler and TESS light-curve analysis and computed an updated transit ephemeris. We fitted a linear function of the orbital period (P) and transit epoch (E): T_C(E) = T₀ + EP. We find a period of P = 3.77840121 ± 0.00000082 (days) and a midtransit time of T_C = 2457934.73871 ± 0.00015 (days). The updated ephemeris is in agreement with the ephemeris reported by Anderson et al. (2014) and does not indicate the presence of transit timing variations Figure 3.

Zoom In Zoom Out Reset image size

To produce spectroscopic light curves, we summed the flux of the target and reference star in bands with variable widths from 80–253 Å, Table 3. The blue and red grisms overlap in the wavelength range from 5300–6200 Å and each covers the sodium D lines at 5890 and 5896 Å. We choose common bins within the overlapping region to allow a direct comparison when combining the blue and red transmission spectra. We enlarged the first three bins in the blue grism and merged two pairs of bins, covering the O₂ A and B telluric bands from 7594–7621 Å and from 6867–6884 Å, respectively. In this way, we increased the signal-to-noise ratio of the broader bands and produced a total of 46 light curves.

We established wavelength-independent, i.e., common-mode, systematic correction factors for each night, similar to our previous studies with FORS2 (Nikolov et al. 2016, 2018; Gibson et al. 2017; Carter et al. 2020; Wilson et al. 2020). To obtain the correction factors, we divided the white-light curve from each night by the transit model computed with the weighted-mean system parameters (Section 4.1). The common-mode factors for each night are shown in Figure 1.

We modeled the spectroscopic light curves with a function accounting for the transit and systematics simultaneously. Before fitting, we divided each spectroscopic light curve by the wavelength-independent systematic correction factors for that night. In the fits, we allowed only the relative planet radius R_p /R_* and the linear limb-darkening coefficient u₁ to vary. We obtained theoretical values for u₁ and u₂ following the same approach as for the white light curves. We also fitted for both limb-darkening coefficients and found that the uncertainty of u₂ is large and consistent with the theoretical prediction. We interpret this as an indication for insufficient constraining power of the data for the quadratic coefficient. Given the transmission spectra did not substantially change we chose to fix u₂ and to fit only for u₁.

Zoom In Zoom Out Reset image size

We accounted for the systematics using a low-order polynomial (up to a second degree with no cross terms) of dispersion and cross-dispersion drift, air mass, FWHM, and the rate of change of the rotator angle. We produced all possible combinations of detrending variables and performed separate fits including each combination within the systematics function. This approach has been prefered as opposed to GP regression, as the common-mode-corrected VLT spectroscopic light curves exhibit a lower level of systematic effects. Our choice also rests on results from a VLT comparative follow up of WASP-39b with atmospheric features detected with the Hubble Space Telescope (Nikolov et al. 2016). We computed the Akaike Information Criterion for each fit and estimated the statistical weight of the model depending on the number of degrees of freedom (Akaike 1974). We chose to rely on the Akaike Information Criterion instead of other information criteria, for example the Bayesian Information Criterion (Schwarz 1978), because the Akaike Information Criterion selects more complex models, resulting in more conservative error estimates. We marginalized the resulting relative radii and linear limb-darkening coefficient following Gibson (2014). We found that for most light curves systematics models parameterized with linear air mass, dispersion drift, and FWHM terms resulted in the highest statistical weight.

Before fitting each light curve, we set the spectrophotometric uncertainties of each band to values that include photon and readout noise. To determine the best-fit models, we used a Levenberg–Marquardt least-squares algorithm and rescaled the uncertainties of the fitted parameters using the dispersion of the residuals (Markwardt 2009). Residual outliers larger than 3σ, typically 2–3 per light curve, were excluded from the analysis. Correlated residual red noise was accounted following the methodology of Pont et al. (2006) by modeling the binned variance with the relation ${\sigma }^{2}={\sigma }_{w}^{2}/N+{\sigma }_{r}^{2}$ relation, where σ_w is the uncorrelated white noise component, N is the number of measurements in the bin, and σ_r is the red noise component. We found white noise dispersion in the range from about 450–900 parts per million. For the red noise, we found a dispersion in the range from about 40–70 parts per million. The resulting time series are shown in Figures 10 and 11.

Activity and rotation are known to complicate the interpretation of transmission spectra and in particular when combining multi-instrument, multi-epoch data sets (Huitson et al. 2013; McCullough et al. 2014; Pinhas et al. 2018). To assess the level of stellar activity and photometric variability, we inspected the cores of the Ca ii H&K lines and photometry time series from archival data. We analyzed five MPG/ESO 2.2 m telescope high-resolution spectra obtained with the FEROS spectrograph (ESO program: 099.A-9010, PI P. Sarkis) at a resolution of R = λ/δλ = 48,000, Table 4. The spectra have been obtained using a fiber with a diameter of 2'' and reduced with ESO's FEROS pipeline to phase three products, i.e., extracted one-dimensional spectra with a wavelength solution. A comparison between the tabulated and observed wavelengths of the Ca ii H&K line cores indicates an offset (Figure 5). This can be attributed to the low signal to noise of the individual spectra, which reach ∼5 in the blue wavelengths, but can also be an indication for the presence of emission in the core of the lines. The activity index, S_HK, was estimated from the spectrum using the method described by Vaughan et al. (1978) and transformed to the $\mathrm{log}R{{\prime} }_{\mathrm{HK}}$ system using the relations from Noyes et al. (1984). We measure $\mathrm{log}R{{\prime} }_{\mathrm{HK}}$ = −${4.9}_{-0.24}^{+0.09}$, which is consistent with low-to-moderate activity of the host star.

Zoom In Zoom Out Reset image size

We analyzed an archival light curve consisting of 250 photometric data points from the Ohio State University's All-sky Automated Survey for Supernovae (ASAS-SN) Photometry Database (Shappee et al. 2014; Jayasinghe et al. 2019). The data has been obtained during the period from 2014 May to December and from March to December in the years from 2015–2018, respectively (Figure 12). We performed Lomb–Scargle periodogram analysis with trial frequencies ranging between 0.005 and 1 day⁻¹, corresponding to a period range between 1 and 200 days. We find a rotational modulation with a period of 13.74 days with a false alarm probability of ∼0.25. We phase folded the photometry and fitted a sinusoid to measure the amplitude of the flux variation. We found an amplitude of 0.5% ± 0.1% and residual scatter of 12.8 parts-per-thousand (ppt), Figure 4. Our result agrees with the 4 mmag (0.4%) upper limit for flux variability of WASP-110 reported by Anderson et al. (2014). Anderson et al. (2014) also found that WASP-110 is an evolved 8.6 ± 3.5 Gyr G9 star, which implies a gyrochronologic rotation comparable or exceeding the equatorial rotation period of the Sun i.e., ≳25 days. Anderson et al. (2014) report $v\sin ({i}_{* })=0.2\pm 0.6$ km s⁻¹, which translates to a rotation period P_*(i_* = 90°) ≳ 34 days. While this is a discrepancy with our finding, a detailed comparison is hampered by the unconstrained stellar axial tilt, which can significantly decrease the actual rotation period. It is worth mentioning that empirical relations of the rotation periods to the chromospheric activity level $\mathrm{log}R{{\prime} }_{\mathrm{HK}}$ predict ∼22 days for WASP-110 (Suárez Mascareño et al. 2015, their Figure 13). Given the spread of ∼10 days in the empirical relationship for the rotation period for G-type stars with −$0.1\lt \left[\mathrm{Fe}/{\rm{H}}\right]\lt 0.1$, the predicted estimate is consistent with our measured period. This is also the case for $v\sin ({i}_{* }\sim 4^\circ )$ of Anderson et al. (2014).

The measured transmission spectrum of WASP-110b is presented in Figures 6 and 7. To combine the transmission spectra from the blue and red grism observations, we computed the weighted mean of both data sets within the overlapping wavelength region. We find a marginally significant ∼1.5σ radius difference from the two observations of ΔR_p /R_* = (33 ± 22) × 10⁻⁴. This corresponds to a depth variation of ∼11 ppm, which is within the 0.5% ± 0.1% stellar variation of the ASAS-SN photometry. The blue and red pairs of radius measurements show an excellent agreement within the overlapping region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The spectrum is marginally flat without an increased absorption from sodium or potassium, which are the main absorbers expected in the optical transmission spectra of irradiated gas giants at T_eq ∼ 1100 K similar to WASP-110b. An increased absorption is present in the wavelength range from ∼0.63 to ∼0.75 μm, spanning ∼3 pressure scale heights. To interpret the spectrum we turn to model comparison and retrieval analysis.

We compared the transmission spectrum with clear, cloudy, and hazy synthetic spectra with solar abundances from the model presented in Fortney et al. (2008, 2010). These spectra include a self-consistent treatment of radiative transfer and chemical equilibrium of neutral and ionic species. Chemical mixing ratios and opacities were computed assuming solar metallicity and local chemical equilibrium, accounting for condensation and thermal ionization but not photoionization (Lodders 1999, 2002; Freedman et al. 2008, 2014).

Similar to our previous studies, the cloudy and hazy models were computed using a simplified treatment of the scattering and absorption to simulate the effect of small particle haze aerosols and large particle cloud condensates at optical and near-infrared wavelengths. In the case of haze, Rayleigh scattering opacity $(\sigma ={\sigma }_{0}{(\lambda /{\lambda }_{0})}^{-4})$ has been assumed with a cross section that was 1000× the cross section of molecular hydrogen gas (σ₀ = 5.31 × 10⁻²⁷ cm² at λ₀ = 3500 Å (Lodders 1999). To account for the effects of a flat cloud deck, we included a wavelength-independent cross section that was a factor of 100× the cross section of molecular hydrogen gas at λ₀ = 3500 Å.

We obtained the average values of the models within the wavelength bins of the observed transmission spectrum and fitted these theoretical values to the data with a single parameter accounting for their vertical offset. The χ² and Bayesian Information Criterion statistic quantities were computed for each synthetic spectrum with the number of degrees of freedom for each model determined by ν = N − m, where N is the number of measurements and m is the number of free parameters in the fit. We find the cloudy spectra to best describe our observations (χ² = 61) followed by the haze and cloud-free simulated spectra (Figure 7).

We additionally considered inverse methods to constrain the atmospheric properties of WASP-110b via atmospheric retrievals using the AURA code (Pinhas et al. 2018). AURA computes model-transmission spectra using line-by-line radiative transfer in a plane parallel atmosphere in hydrostatic equilibrium. The forward model-generating component is interfaced with the PyMultiNest Nested Sampling package (Buchner et al. 2014), which allows for Bayesian parameter estimation and model comparison. We describe the terminator pressure–temperature profile using the 6-parameter prescription of Madhusudhan & Seager (2009). Our model atmosphere includes opacities arising from H₂–H₂ and H₂–He collision-induced absorption, as well as absorption from H₂O (Rothman et al. 2010) and Na and K (Welbanks et al. 2019); the latter two accounting for the effects of pressure broadening. We avoid including high-temperature species in our atmospheric model (e.g., TiO, VO and other metal hydrides/oxides, ions, etc.), given the relatively low equilibrium temperature of WASP-110 b. Inhomogeneous coverage of clouds and hazes is incorporated in our model using the prescription of Pinhas et al. (2018), treating clouds as gray opacity below the cloud deck and hazes as a modification to Rayleigh scattering above the cloud deck.

AURA also incorporates the effects of stellar heterogeneity, being the first retrieval code demonstrated to successfully include such effects (Pinhas et al. 2018). The stellar heterogeneities are described using three parameters: δ, the fraction of the projected stellar disk covered by cooler starspots/hotter faculae; T_het, the temperature of the starspots/faculae; and T_phot, the temperature of the immaculate photosphere. We construct the spectral contributions of the two stellar surface components by interpolating through a grid of PHOENIX spectra (Husser et al. 2013), fixing the stellar metallicity and gravity to literature values. Hotter faculae/cooler starspots result in a larger/smaller apparent stellar radius toward the blue end of the spectrum, which in turn results in lower/higher observed transit depths at shorter wavelengths.

We find that our retrieval produces a good fit to the data, as shown in Figure 8. We successfully obtain constraints for the three stellar heterogeneity parameters, finding $\delta ={0.20}_{-0.09}^{+0.13}$, ${T}_{\mathrm{het}}={5643}_{-134}^{+195}$ K and ${T}_{\mathrm{phot}}={5387}_{-81}^{+85}$ K. We show the retrieved posterior distribution for these three parameters in Figure 9. Our retrieved T_het value is higher than both the retrieved T_phot and literature T_eff values, indicating the presence of hotter faculae. We determined the detection significance (D.S.) for stellar heterogeneity by carrying out a Bayesian model comparison with the findings of a second retrieval that excludes the effects of stellar heterogeneity. We find a D.S. of 3.1σ, indicating a strong detection. This stems from the overall trend in the data of increasing transit depths at longer wavelengths, which are attributed to unocculted stellar faculae rather than opacity sources in the planetary atmosphere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We additionally consider separate retrievals for the grism 600B and 600RI data sets, without applying any offsets to either. We find that the constraints on the cloud and heterogeneity parameters are governed primarily by the 600RI data set, with the 600B data set providing only weak constraints. The constraints from both data sets are consistent with each other and with the joint data set to within the 1σ uncertainties. We therefore find that the two data sets can be combined for a joint retrieval with a single set of heterogeneity parameters without significant loss of information.

Our retrieval is unable to detect any of the chemical species considered, including H₂O. We additionally obtain weak constraints for a high-altitude cloud deck pressure, $\mathrm{log}({P}_{\mathrm{cloud}})\,=-{3.4}_{-1.3}^{+1.8}$ (in bar), while also finding that our retrieval shows a preference for near-full cloud coverage fractions. As before, we conducted a third retrieval, this time excluding clouds and hazes to determine their D.S.. We find that our retrievals detect clouds and hazes with a 2.1σ significance, suggesting only a tentative detection. Clouds are only tentatively detected because a featureless spectrum can be caused by either clouds or low abundances of alkali metals. The overall posterior distribution is shown in Figure 13.

We also conduct further retrievals to explore how more elaborate model considerations affect our findings. In one case, we additionally consider several chemical species with significant optical opacities, such as TiO, VO and AlO (Fortney et al. 2008; Gandhi & Madhusudhan 2019), even though they are not expected in thermochemical equilibrium at the low atmospheric temperatures of WASP-110 b. While we find nominal indications of these species, none has a D.S. greater than 2σ to constituting a strong detection. We also run retrievals with an offset applied to either the 600RI or 600B data sets as an additional retrieved parameter. In this case, the D.S. of the heterogeneity is reduced to ∼1–2σ, depending on which of the two data sets is being offset, while clouds and hazes are detected at a D.S. of ∼2σ.