MOPSS. II. Extreme Optical Scattering Slope for the Inflated Super-Neptune HATS-8b

IMACS is a wide-field imager and optical multi-object spectrograph located on the 6.5 m Magellan Baade Telescope at Las Campanas Observatory. We use the f/2 camera on IMACS and custom-made observing masks to simultaneously observe our target and a number of calibration stars. The IMACS f/2 CCD consists of eight separate rectangular chips in a 4 × 2 arrangement. There is a gap of ∼57 pixels between the short edges of the chips, and ∼92 pixels between the long edges of the chips.

We use the 300 lines/mm grism at a blaze angle of 175, resulting in theoretical wavelength coverage from 3900 to 9000 Å, but practically, from ∼4600 to ∼8000 Å. The shortened wavelength coverage on the blue end is due to the use of a blocking filter for light below 4550 Å. The throughput at these wavelengths is low, so we choose this addition in order to eliminate spectral contamination from 7800 to 9000 Å from the grism's second order. On our first night we did not use this filter, so we do not use wavelengths past 7800 Å in our analysis. Further, the gap between chips occurs at ∼7800 Å for our target star, with the main O₂ absorption in our atmosphere at ∼7600 Å. In our analysis, we skip any bins that contain either of these features.

As mentioned above, we use a custom-made observing mask designed to maximize the number of calibrator stars, all aligned in a way to allow the full wavelength span to fall across the chips. The main target is placed on a central chip (see Figure 1). Calibrators were selected based on their magnitudes and spectral type if available, or their color in the common case that a spectral type was not determined. We select those stars that are most similar in color and are within ∼0.75 mag of the main target when possible. Each calibrator and the main target is centered in a wide slit (15'') with large lengths (20'') so that we are not concerned with slit losses and to improve background subtraction. With this setup our resolution is limited by the seeing each night (see Section 2.2 for seeing data). Our observational efficiency was limited by the readout time for the CCD (82 s in fast mode and 1 × 1 binning, 29 s in fast mode and 2 × 2 binning).

Zoom In Zoom Out Reset image size

HATS-8b is an ∼1300 K hot Neptune orbiting a G-type star with a V mag of 14.03. Discovered in 2015 (Bayliss et al. 2015), HATS-8b was immediately noted to be an ideal target for transmission spectroscopy. However, as of writing, there are no previous transmission spectra observations in the literature.

HATS-8b was observed on the nights of 2017 July 23 in 1 × 1, binning with exposure times of 300 s, and 72 exposures and 2017 August 11 in 2 × 2 binning, with exposure times of 180 s and 92 exposures. This corresponds to an observational efficiency of 78% and 86%, respectively. Seeing was between ∼08 and 15 throughout the first night, with the seeing improving as the night went on. Our wide slits mean we were not concerned with slit losses, even at the 15 seeing. We were not impacted by clouds or wind throughout the night. During the second night, seeing was around 06–07 throughout the night with no cloud coverage. In addition to our main target, we observed an additional seven calibrator stars which are listed in Table 1. We reserve one of the calibrators (star #3 in Figure 1 and Table 1) to serve as a check that our calibrators have a constant flux throughout the night. This allows us to ensure that our pipeline is correctly accounting for airmass, seeing, and other instrumental effects throughout the night.

Our reduction pipeline was developed in Python and is described in detail in MOPSS1. It follows the traditional techniques for spectral reduction. We follow the process outlined in Nikolov et al. (2014) for cosmic-ray removal. We use SpectRes (Carnall 2017) for fast flux-conserving binning when down sampling the spectra to lower resolutions.

Wavelength calibration is done using a second mask with small, 1'' square slits at the locations of each object. We use a HeNeAr lamp and take several frames at the beginning of the night. By measuring shifts in the spectral direction as described in MOPSS1, we can shift the spectrum back to the wavelength frame of the first exposure, which should most exactly match the pixel locations of our arc frames. This paper includes improvements to our treatment of biases introduced due to atmospheric and instrumental effects, as described in the following sections.

Because our red-noise model is best applied to normalized data, we first correct for extinction due to airmass. This process requires fitting for extinction as a function of wavelength following the relationship

$$\begin{eqnarray}&&{m}_{\mathrm{observed}}(\lambda ,t)={m}_{\mathrm{true}}(\lambda ,t)+k(\lambda )Z(t)\end{eqnarray} \tag{ 1 }$$

where m_observed is the observed magnitude at a given time and wavelength, m_true is the true magnitude of the object, k(λ) is the wavelength-dependent extinction coefficient, and Z(t) is the airmass as a function of time. We can use the observed magnitude and airmass from our first exposure and the definition of magnitudes to write Equation (1) as

$$\begin{eqnarray}&&{F}_{\mathrm{observed}}(t)={F}_{{Z}_{0}}{10}^{-k(\lambda )(Z(t)-{Z}_{0})/2.5}\end{eqnarray} \tag{ 2 }$$

where F_observed(t) is the number of counts detected at time t, ${F}_{Z0}$ is the number of counts at time = 0, k(λ) is the wavelength-dependent extinction, Z(t) is the airmass at time t, and Z₀ is the airmass at t = 0.

We fit for the average k(λ) across all objects except the target (the transit event will bias the fit) and then remove the trend described in Equation (2) from all objects. Figure 2 shows an example of the extinction curve for the white light binning of our check star on night one.

Zoom In Zoom Out Reset image size

Typically, this trend is addressed simply through the division of the target star by the calibrator star(s). However, because we wish to model the correlated noise component in each star independently, we choose to remove the airmass trend in the manner described here. We find that treating each star independently in our noise model provides a greater precision in our final results.

We expect the dominate source of correlated noise in our data to be caused by variations in seeing and the spectra shifting on the chips throughout the night. Though IMACS has very stable pointing, our night-long coverage of the object results in the objects shifting on the chips throughout the night at a measurable level of a few pixels. We extract these shifts during our data reduction process so that we can decorrelate the data against these shifts. The seeing is calculated by converting the spatial FWHM measured when flattening the 2D spectra to 1D, multiplied by the detector's plate scale of 02/pixel.

In addition to the spectral and spatial shifts and seeing, we decorrelate the data against the background counts as well. We use a linear combination of all of the above parameters fit to the out-of-transit times to predict in-transit data. Our model is described as

$$\begin{eqnarray}&&f={\boldsymbol{X}}\beta \end{eqnarray} \tag{ 3 }$$

where f is the data, X is a matrix containing the spatial and spectral positions relative to time = 0, seeing, background counts, and a column of unity values to account for a constant offset, and β is an array containing the relative importance of each term. By inverting Equation (3),

$$\begin{eqnarray}&&\beta ={\left({{\boldsymbol{X}}}^{T}{\boldsymbol{X}}\right)}^{-1}f,\end{eqnarray} \tag{ 4 }$$

one can solve for β during the out-of-transit times and use the results to calculate the predicted flux during transit under the current conditions (chip location, seeing, background).

We perform this fit for each object and bin independently. Though relative pixel shifts are the same for all objects, due to differences in pixel response, a shift of 1 pixel in the spatial direction may mean an increase in counted photons for one object or wavelength bin, but a decrease for another object or wavelength bin. Figure 3 shows an example of our noise model as applied to the white light data for our check star on night two. We are able to explain the majority of the trends and scatter with this method. Any trends left over we attribute to instrumental effects we are unaware of or do not include. Theses trends are small in comparison to what we fit for here and are removed through dividing our target by a master calibrator star.

Zoom In Zoom Out Reset image size

To generate our light curves, we first create our master calibrator which consists of all reference stars except the one chosen as the check star. The master calibrator allows us to correct for instrumental effects not captured in our noise model. After applying the master calibrator, a small systematic generally remains in the light curve. Though this trend is smaller than in previous work due to our airmass and red-noise corrections, we still find it necessary to include a baseline fit. For night one, we fit this as a third-order polynomial in time, as detailed in MOPSS1. For night two, we use a second-order polynomial in time due to there being relatively few data points before transit to constrain the fit.

3.4.1. White Light Curves

The white light curves for HATS-8b are generated from 4200 to 8000 Å on night one and 4600–9000 Å on night two. The difference is due to a blocking filter we implemented on night two as described in Section 2.2. The white light curves are used to fit for orbital parameters in the same manner as in MOPSS1. We use Batman (Kreidberg 2015) with a quadratic limb-darkening function, as well as emcee (Foreman-Mackey et al. 2013) to fit for center of transit, period/semimajor axis, inclination of orbit, white-light planet radius, and white-light limb darkening. All orbital parameters are fit with a Gaussian prior defined by the values and errors given in Table 2. We find that our Markov Chain Monte Carlo (MCMC) chains do not converge as easily if we have both the period and semimajor axis as free parameters, and so we maintain the relation between period and semimajor axis by only fitting for period while simultaneously calculating a corresponding semimajor axis based on the stellar mass and radius and their errors cited in Table 2.

Table 2.  HATS-8b: Stellar and Orbital Parameters

Stellar Parameter	Value
Mass, (M⊙)	1.056 ± 0.037
Radius, (R⊙)	${1.086}_{-0.059}^{+0.149}$
Teff, (K)	5679 ± 50
[Fe/H], (dex)	0.210 ± 0.080
ξ, (km s−1)	2.00 ± 0.50
(microturbulent velocity)
$\mathrm{log}g$, (cm s−1)	4.386 ± 0.071
Planet Parameter	Value
Mass, (MJup)	0.138 ± 0.019
Radius, (RJup)	${0.873}_{-0.075}^{+0.123}$
Teq, (K)	${1324}_{-38}^{+79}$
Period, (days)	3.583893 ± 0.000010
Semimajor axis, (au)	0.04667 ± 0.00055
Eccentricity	<0.376
Inclination, (deg)	${87.8}_{-1.8}^{+1.2}$

Note. All values have been adopted from Bayliss et al. (2015). Our fits use an eccentricity of 0.0 due to the lack of a strong constraint from previous work.

Download table as:  ASCIITypeset image

In all light curve fits we check for convergence by calculating the autocorrelation time (τ) of our chains for each parameter. We use the autocorrelation function suggested by the emcee package from Goodman & Weare (2010). We run our chains long enough that τ approaches a consistent value regardless of chain length. We find that for all bin sizes (white and binned in wavelength space) all parameters converge to τ ∼ 100 steps fairly quickly. Our chains are run at least 500× longer than τ.

For all fits we run an initial short chain and use the resulting reduced χ² to rescale the uncertainties on our light curves such that we obtain a reduced χ² of ∼1 on our full production runs. This rescaling accounts for any underestimated uncertainties and produces more reliable uncertainties on the fit transit depth.

3.4.2. Binned Light Curves

We bin the data from each night into bins of width 400, 200, 100, 50, 40, and 20 Å in order to capture the overall shape of the transit at high precision, while also searching for absorption features that may be averaged over in wider binnings. For each resultant light curve, we use emcee to fit for R_p/R_S and the two quadratic limb-darkening parameters. We do not fit for orbital parameters as a function of wavelength as they should not have a wavelength dependence. We use the orbital parameters we retrieve from our white light emcee runs. Figure 4 shows the 400 Å binned light curves and their fits from our MCMC runs.

Zoom In Zoom Out Reset image size

Our initial quadratic limb-darkening parameters are interpolated between the Johnson filter values from Claret & Bloemen (2011). While it is understood that these may not be completely accurate, our use of them as a starting value does not affect our results since we use a noninformative flat prior which allows our walkers to fully explore the parameter space without being biased to our initial guess. We confirm this by also starting white light chains at values of 0.0 for both limb-darkening parameters and the planet radius, which return the same results as those that begin at the Claret & Bloemen values. Because the unconstrained runs take longer to converge, we choose to use the starting guesses from Claret & Bloemen and a more informed initial guess for the planet radius. For this planet, we find that our returned limb-darkening values agree well with those from Claret & Bloemen (2011). Figure 5 shows the theoretical quadratic limb-darkening coefficient values and our fit values for both nights.

Zoom In Zoom Out Reset image size

We fit our data from each night independently. We present results here for the 400 Å binnings, which demonstrate the overall scattering slope. We find no definitive evidence of sodium or potassium absorption at any of our finer binnings. Figure 6 shows our transmission spectra and scattering slopes. Table 3 lists our results for both nights.

Zoom In Zoom Out Reset image size

Table 3.  HATS-8b: MCMC Fit Results

	2018 July 23
Bin Center	Rp/Rstar	Δ Rp/Rstar	q0	Δ q0	q1	Δ q1	RMSE (ppm)	Photon (ppm)	RMSE/Photon	${\chi }_{\mathrm{red}}^{2}$
4800 Å	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
5200 Å	0.0845	0.0054	0.611	0.065	0.17	0.02	1139.35	345.41	3.3	1.017
5600 Å	0.0741	0.0052	0.635	0.058	0.187	0.025	1040.31	310.21	3.35	1.027
6000 Å	0.0719	0.0053	0.466	0.057	0.195	0.026	984.31	295.99	3.33	1.022
6400 Å	0.0719	0.0056	0.367	0.053	0.269	0.028	974.4	292.16	3.34	1.012
6800 Å	0.0694	0.0034	0.493	0.046	0.209	0.027	629.04	291.02	2.16	1.012
7200 Å	0.0623	0.0056	0.392	0.041	0.213	0.028	780.43	302.87	2.58	1.014
7600 Å	0.062	0.0063	0.382	0.043	0.25	0.03	957.18	335.13	2.86	1.037
8000 Å	0.0542	0.0047	0.375	0.038	0.278	0.03	933.42	363.73	2.57	1.053
8400 Å	0.0509	0.0061	0.372	0.037	0.26	0.029	1383.42	375.9	3.68	1.054

2018 August 28
Bin Center	Rp/Rstar	Δ Rp/Rstar	q0	Δ q0	q1	Δ q1	RMSE (ppm)	Photon (ppm)	RMSE/Photon	${\chi }_{\mathrm{red}}^{2}$
4800 Å	0.085	0.0038	0.571	0.061	0.13	0.015	747.54	426.19	1.75	1.31
5200 Å	0.0831	0.0046	0.487	0.056	0.15	0.02	821.12	350.88	2.34	1.018
5600 Å	0.0823	0.0025	0.434	0.051	0.183	0.024	655.54	317.37	2.07	1.02
6000 Å	0.0778	0.0026	0.406	0.052	0.181	0.026	731.19	302.65	2.42	1.022
6400 Å	0.0742	0.0031	0.362	0.052	0.205	0.029	724.4	298.44	2.43	1.022
6800 Å	0.0745	0.0033	0.347	0.049	0.216	0.028	785.47	297.77	2.64	1.016
7200 Å	0.0775	0.0029	0.341	0.046	0.212	0.028	615.89	313.4	1.97	1.019
7600 Å	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
8000 Å	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯
8400 Å	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

Note. Here q0 is the first quadratic limb-darkening parameter and q1 is the second quadratic limb-darkening parameter. Night one does not have values reported past Earth's 7600 O₂ absorption feature due to second-order contamination. For both nights, the bin centered on this absorption feature is discarded. The column labeled reduced mean squared error (RMSE) contains the RMSE of the light curve fit in parts per million, the column labeled photon contains the photon noise level in parts per million, and the column labeled RMSE/photon contains a value describing how close we are to photon limited observations.

Download table as:  ASCIITypeset image

Our data covers only optical wavelengths, which can only constrain combinations of atmospheric properties by itself. Due to the lack of absorption features, we report only measurements of the strength of the scattering slope. The optical slope we detect is typically a signature of scattering, given by

$$\begin{eqnarray}&&\displaystyle \frac{1}{{R}_{{\rm{S}}}}\displaystyle \frac{{{dR}}_{{\rm{p}}}}{d\mathrm{ln}\lambda }=\alpha \displaystyle \frac{1}{{R}_{{\rm{S}}}}\displaystyle \frac{{k}_{{\rm{B}}}{T}_{\mathrm{limb}}}{\mu g}\end{eqnarray} \tag{ 5 }$$

first described in Lecavelier Des Etangs et al. (2008), where R_S is the stellar radius, R_p is the planetary radius, α is the power of the wavelength dependence in the scattering cross section (−4 for Rayleigh scattering), k_B is the Boltzmann constant, T_limb is the temperature at the limb of the planet, μ is the mean molecular weight of the atmosphere, and g is the planetary gravity. 1/R_S is included on both sides of the equation because transmission spectra are typically plotted as the planet radius relative to the stellar radius.

Our main unknowns in Equation (5) are the atmospheric mean molecular weight (μ) and the limb temperature (T_limb). Typical values for μ range from 2.2 for Jupiter to 2.5–2.7 for Neptune. If the atmosphere of the planet is clear, the gas should produce Rayleigh scattering and α will be equal to −4 as on Earth. However, if there exist particles in the atmosphere large enough to produce Mie scattering, α can take on a wider range of values.

Our measured transmission slope on night one, 2017 July 23, is (1/R_S) (dR_p/d ln λ) = −0.142 ± 0.024, ∼27 times stronger than nominal Rayleigh scattering, defined as the right-hand side of Equation (5) with α = −4, T_limb = T_eq, μ = 2.3, and g = 464.51 cm s⁻² (calculated from the planetary parameters given in Table 2). With these parameters, there is no reasonable combination of μ and limb temperature to explain the slope at α = −4.0. At the equilibrium temperature of the planet and a mean molecular weight of 2.3, we can match the slope with an extremely strong wavelength dependence of λ⁻⁹⁵ for scattering.

On night two, the measured slope is $(1/{R}_{{\rm{S}}})({{dR}}_{{\rm{p}}}/d\mathrm{ln}\lambda )=-0.057\pm 0.024$, approximately 40% of that from night one. Prior to searching parameter space for a combination of μ, T_limb, and α that matches this slope, we first discuss the effect unocculted starspots may have on the measured scattering slopes.

An important consideration is how uncertainties on planetary properties may affect the nominal Rayleigh slope, which in turn affects our expectations. Taking our parameters and uncertainties from Table 2, we note that R_planet has the highest relative uncertainty. When estimating the Rayleigh slope, planet radius matters only for the planet gravity. This results in a factor of ${R}_{\mathrm{planet}}^{2}$. If the true planet radius sits at the edge of the upper error bar, so that R_planet = 0.996 × R_Jupiter, the nominal Rayleigh slope will increase by a factor of 1.3. This is not a significant change when compared to the factor of ∼27 we are looking to explain and so we do not explore this option further.

In addition, HATS-8b is a low-density, low gravity exoplanet, with ρ = 0.259 g cm⁻³, 2.6× lower than Saturn. It is possible that the extended atmosphere requires a different treatment for atmospheric models. In particular, Exo-Transmit (Kempton et al. 2017) does not account for variations in scale height with altitude, a potentially important factor for these low-density, inflated planets.

As noted above, we also search for the presence of sodium absorption (5890 Å) and potassium absorption (7665 Å) with our narrow bins. We find that the presence or lack of potassium absorption could not be well established due to the deep O₂ telluric line at ∼7600 Å diluting any signal. We do not find strong evidence of sodium absorption.

Although the star HATS-8 is known to be an inactive star from exoplanet discovery observations (Bayliss et al. 2015) and we do not detect an occulted spot during either of the observed transit events, we cannot rule out that the host star has some amount of spot coverage. Particularly because unocculted spots can appear as a scattering slope in transmission spectrum, and our results suggest a stronger than expected slope with differences between the two nights. If the rotational period of the host star is approximately twice the time between observations, we would be viewing the planet against a different side of the star and different spot coverage.

Following the approach of May et al. (2018), we use the formulation of Louden et al. (2017) to investigate the influence unocculted spots may have on our results. Here we assume that night one is heavily influenced by spots, while night two is relatively uninfluenced by spots because of the stronger slope measured on night one. The measured transit depth, δ_m(λ), is described as a function of the true transit depth, δ(λ), the spot coverage fraction, η, and the relative flux from a spot, F_λ(spot), compared to the stellar photosphere, F_λ(star).

$$\begin{eqnarray}&&{\delta }_{m}(\lambda )=\delta (\lambda )\displaystyle \frac{1}{1-\eta \left(1-\tfrac{{F}_{\lambda }\left(\mathrm{spot}\right)}{{F}_{\lambda }\left(\mathrm{star}\right)}\right)}.\end{eqnarray} \tag{ 6 }$$

The overall dimming of the star due to unocculted spots can be written as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{\lambda }(\mathrm{star},\mathrm{corrected}) & = & \eta {f}_{\lambda }{F}_{\lambda }(\mathrm{star})+(1-\eta ){F}_{\lambda }(\mathrm{star})\\ & = & [1-\eta (1-{f}_{\lambda })]{F}_{\lambda }(\mathrm{star})\end{array}\end{eqnarray} \tag{ 7 }$$

with f_λ as the spot contrast given as F_λ(spot)/F_λ(star). When the star dims due to unocculted spots, the transit depth looks relatively larger compared to transits with no unocculted spots. Because the spot's blackbody spectra peaks at longer (redder) wavelengths, while the star peaks at relatively shorter (bluer) wavelengths, the spot contrast level is higher at short wavelengths, and lower at red wavelengths. This is therefore a wavelength-dependent effect, with blue wavelengths more heavily affected than red, resulting in a star-induced slope being injected into the data.

Following Equation (6), we calculate the required δ_m(λ)/δ(λ) to match the observed night one slope to both the nominal Rayleigh slope and the observed slope on night two. Figure 7 shows the corrected/expected slopes (where a ratio of 1 can fully explain the slope) at a variety of spot coverage fractions (η) and ΔT = T_star–T_spot.

Zoom In Zoom Out Reset image size

The required spot coverage fraction to fully explain our night one slope versus nominal Rayleigh slope at a ΔT of 1500 K is ∼80%, much higher than one would expect. HATS-8 is a G star with an age of 5.1 ± 1.7 Gyr (Bayliss et al. 2015), comparable to the Sun. On magnetically active low-mass stars, Jackson & Jeffries (2013) suggest that the spot coverage could be as high at 40% with T_spot/T_star = 0.7 (similar to our value). Further, long term spot-coverage studies by Alekseev & Kozhevnikova (2018) find spot coverages up to 42% for 13 active G and K stars. In agreement with these literature values for G stars, we can explain the difference between our two nights with a spot coverage fraction of 40% and a ΔT of 1500 K. We note, however, that HATS-8 is a slowly rotating G dwarf and aside from metallicity is very similar to our Sun (Bayliss et al. 2015). Typical spot coverage fractions for inactive stars such as the Sun range from  0.05% to 0.5% at solar minimum and solar maximum, respectively. (Bogdan et al. 1988; Solanki & Unruh 2004). Further, if HATS-8 had 40% spot coverage on night one, it is rather unlikely that the observed transit did not exhibit a spot crossing event unless the orbit and rotational axis of the star are very well aligned and the spots are strongly bound to latitudes away from the path of the transit. In addition, a 40% spot coverage fraction on night one is even more unlikely because our two observations are only one month apart. Therefore, because our two slopes are consistent within 1.3σ, we instead move forward by computing a weighted average of the slopes of our two nights. Our weighted slope is then $(1/{R}_{{\rm{S}}})({{dR}}_{{\rm{p}}}/d\mathrm{ln}\lambda )=-0.091\pm 0.017$.

With this slope, we explore the range of atmospheric parameters that can explain the observed slope in a clear atmosphere. In Figure 8 we show a range of possible parameters to satisfy the observed combined scattering slope in our data. Based on the unreasonable parameter values shown in Figure 8, we require a strong wavelength dependence of approximately λ⁻⁵⁰ and a low mean molecular weight atmosphere (μ ≲ 2.0) to explain our slope. This is inconsistent with a clear atmosphere, where α must equal −4. Further, because a μ < 2.0 suggests a large fraction of the atmosphere is atomic rather than molecular (unphysical at these temperatures), we choose to exclude solutions for such low mean molecular weight values.

Zoom In Zoom Out Reset image size

As described below, we prefer an explanation that involves a portion of the slope due to clouds (see Section 4.4), and the remainder resulting from differences between our expected and realized values of mean molecular weight and the limb temperature.

Naturally, the next consideration is that the atmosphere of HATS-8b is not clear, and instead contains scattering particulates of some kind. The well-known hot Jupiter HD 189733b is a natural comparison here, due to its strong optical slope and similar equilibrium temperature (1200–1400 K). A detailed analysis of the optical slope for HD 189733b is presented in Pont et al. (2013). As demonstrated in Sing et al. (2016), HD 189733b has the strongest scattering slope among hot Jupiters observed with the Space Telescope Imaging Spectrograph (STIS). Of importance, however, is that Pont et al. (2013) state that the optical slope for HD 189733b can be explained by an α = −4 Rayleigh scattering slope at a temperature of ∼1300 K for particle sizes below 0.1 μm. This is not the case for HATS-8b, as the slope is steeper than α = −4 even after considering the possibility of an extreme level of unocculted starspots. Therefore we note that the optical scattering slope of HATS-8b is the strongest detected scattering slope to date.

Pinhas & Madhusudhan (2017) explore the effects of a variety of species on the wavelength dependence (α) of the scattering cross section. In their work, Pinhas & Madhusudhan show that in some cases, the presence of cloud particles can result in a slope as strong as α = −13.

Of the species Pinhas & Madhusudhan discuss, they suggest that transmission spectra, which show strong optical slopes, are likely to contain sulphide clouds—Na₂S, MnS, or ZnS. These three species have condensation temperatures of 1176 K, 1139 K, and 700 K, respectively, at 10⁻³ bar (Wakeford & Sing 2015). It is probable that Na₂S and Mns could condense out in the nightside atmosphere of HATS-8b, which has an equilibrium temperature of 1321 K for an albedo of 0. If a large amount of scattering particles are present, we would expect a higher albedo, and therefore lower temperature that allows the condensation of these species. Cloud formation along the limbs is an expected (Parmentier et al. 2016; Roman & Rauscher 2019) consequence of atmosphere dynamics.

Specifically, Pinhas & Madhusudhan find that MnS produces the steepest slopes at a modal particle size of 10⁻² μm for a cloud scale height (H_c) equal to the gas scale of k_B T_eq/(μ g). Under these conditions, one can attain a scattering slope with α = −13 for a hot Jupiter of g = 24.79 m s⁻², μ = 2.3, and a solar abundance of MnS. The slope models in Pinhas & Madhusudhan (2017) depend on the cloud scale height (a smaller H_c brings the slope closer to the standard Rayleigh slope), the modal particle size based on their assumed distribution (smaller particles cause steeper slopes), the reference pressure (only varies the absolute depth), the grain abundance, and the molecular abundance of the species. Pinhas & Madhusudhan state that the slopes are insensitive to the planet gravity, because the same value is assumed for both the the gas scale height and the cloud scale height.

Additional condensates at these wavelengths and temperatures may include SiO₂, Al₂O₃, Fe₂O₃, Mg₂SiO₄, and MgSiO₃ with condensation temperatures of 1725 K, 1677 K, 1566 K, 1354 K, and 1316 K, respectively, at 10⁻³ bar (Wakeford & Sing 2015). Of these, Pinhas & Madhusudhan find that MgSiO₃ can result in the strongest scattering slope, with an α ∼ −5 for small grains (a modal particle size of 10⁻² μm) if H_c = H. However, this is not much different from a nominal Rayleigh slope of α = −4, and insufficient for the scenario we explore here. We note that this large slope may also be a result of a yet to be identified photochemical haze.

In Figure 9 we show the values of limb temperature and mean molecular weight which, combined with α = −13 for MnS clouds and small particle sizes, could explain our data. We plot the same levels as Figure 8 for comparison.

Zoom In Zoom Out Reset image size

We find that we can explain at most 40% of the observed slope through the inclusion of MnS clouds, but only for a warm, low mean molecular weight atmosphere. However, because a μ < 2.0 suggests a large fraction of the atmosphere is atomic rather than molecular, for reasonable values of μ ≈ 2.0, we can at most explain ∼20% of the slope at the equilibrium temperature of the planet.

Finally, in the limiting case that nothing is changing on the star itself, there must be a different explanation for the varying scattering slope we measure. Our two slopes differ by a factor of 2.5, which, following the methods of Pinhas & Madhusudhan (2017), suggests that the cloud scale height is up to 2.5 times larger on night 1. In order to get the strongest scattering slope (α = −13) we assume H_c = H, which only gets us within a factor of 2 of matching our data. Invoking an increase in the vertical extent of the cloud to enhance the scattering would require vigorous vertical mixing in the atmosphere, but the change between the two observations would require that the mixing also be time variable, which has not been seen in hot-Jupiter atmospheric simulations, to our knowledge. Therefore, if the particulates in the atmosphere of HATS-8b are composed of small MnS particles, it is not likely that we can explain the difference between the two nights as only an increase in cloud coverage on the terminator.