Simultaneous Multiwavelength Variability Characterization of the Free-floating Planetary-mass Object PSO J318.5−22

The unbinned Spitzer light curve (30 s cadence) along with the best-fit sinusoid using a Levenberg–Marquardt least-squares minimization algorithm is presented in Figure 3. The sinusoidal model has four parameters: period (in hr), phase (in degrees), mean light curve value (since we have divided the raw light curve by the median flux over the whole observation, this should tend toward unity), and amplitude (in percent variation, peak to mean light-curve value). The best-fit model with Gaussian noise added is also shown and provides a good match to the observed light curve. We also plot the periodogram in Figure 3 of PSO J318.5−22 as well as a number of reference stars in the field to identify periodic variability. The 1% false-alarm probability (FAP, plotted in blue in the figure) is calculated from 1000 simulated light curves. These light curves are produced by randomly permuting the indices of reference star light curves (Radigan et al. 2014), producing light curves with Gaussian-distributed noise. The observed variability is reasonably well modeled with a sinusoidal model, although successive maxima and minima appear to be marginally increasing in flux. Thus we also considered sinusoidal + linear models.

Zoom In Zoom Out Reset image size

To fully explore the parameter space of both sinusoidal and sinusoidal + linear fits, we used the emcee Markov chain Monte Carlo (MCMC) package (Foreman-Mackey et al. 2013) to determine the full posterior probability distribution. We ran an MCMC chain using a ${\chi }^{2}$ likelihood function with 1000 walkers for 2000 steps. The first 100 steps of each chain were thrown out as part of the burn-in. Chains were checked by eye for convergence. We fit both a single-sinusoid model (results shown in Figure 4) and a sinusoid + linear model (shown in Figure 5). We adopt the 50% quantile value as the best value for amplitude and period. Best values of amplitude (peak-to-median value) and period as well as 68% confidence-interval and 95% confidence-interval errors are presented in Table 1. Both provide reasonably good fits—in Figure 6, we plot 1000 samples from our chain for both models. We calculated the Bayesian information criterion (BIC, Schwarz 1978) for the adopted best-value parameters for both the single-sinusoid and sinusoid + linear fit. The BIC is given by

$$\begin{eqnarray}&&\mathrm{BIC}=\mathrm{ln}(n)\,k-2\,{ln}({ \mathcal L }),\end{eqnarray} \tag{ 1 }$$

where n is the number of data points in the light curve (1700 for the Spitzer curve), k is the number of model parameters, and ${ \mathcal L }$ is the maximized value of the likelihood function of the model. We find $\mathrm{BIC}\simeq 12$ for the single sinusoid and $\mathrm{BIC}\simeq 19$ for the sinusoid + linear fit. The model with the lower value of BIC is preferred, thus we adopt the sinusoid-only model for further analysis and comparison to HST results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  MCMC Fit Results for Spitzer Variability Monitoring

Model	Parameter	Best Value	68%	95%
Sine	amplitude	1.69%	0.03%	0.07%
	period	8.61 hr	0.06 hr	0.11 hr
	phase	−105	24	48
Sine+slope	amplitude	1.81%	0.04%	0.07%
	period	8.63 hr	0.05 hr	0.10 hr
	phase	−84	22	43
	slope	0.05%	0.005%	0.01%

Note. Peak-to-median variability amplitudes are reported here.

Download table as:  ASCIITypeset image

Apai et al. (2017) recently modeled the rapid brightness evolution found in brown dwarf light curves using beating patterns between multiple planetary-scale wave surface features that move at slightly different velocities due to zonal wind speed variations. They considered a simple three-sinusoid model and also applied their Aeolus mapping package, with both spots and planetary-scale waves (Karalidi et al. 2016). We attempted to fit a similar three-sinusoid model to our Spitzer light curves, with the periods of the three sinusoids given as ${P}_{\mathrm{rot}}=({P}_{1}+{P}_{2})/2$ (wavenumber k = 1 waves) and ${P}_{3}={P}_{\mathrm{rot}}/2$ (k = 2 wave), where P_rot is the rotational period of PSO J318.5−22. In this case, the best-fitting model always reverted to a single sinusoid, with negligible amplitudes for the other two sinusoids relative to the uncertainties in the light curve. We only cover approximately two apparent rotational periods for PSO J318.5−22, thus it is unclear from our Spitzer light curve alone whether the nearly sinusoidal variation observed is the fundamental light curve of this object or if we have observed it during a period when multiple planetary-scale wave surface features happen to be in phase. This is qualitatively similar to the case of SIMP 0136 described in Apai et al. (2017), where the apparent periods of the two k = 1 waves are expected to vary only by ∼1%. However, Allers et al. (2016) find a maximum period for PSO J318.5−22 of 10.2 hr from $v\sin i$ measurements, ruling out the possibility of a double-peaked light curve, thus we expect the 8.6 ± 0.1 hr rotational period derived from our single-sinusoid MCMC fits to be accurate in this case.

Allers et al. (2016) constrain the inclination of PSO J318.5−22 to >29° based on their measured $v\sin i$, reasonable estimates of the radius of PSO J318.5−22 from on evolutionary models, and a lower limit on the period of ∼5 hr from Biller et al. (2015). With our high-precision measurement of the period from Spitzer observations, we can now directly measure the inclination for this object. Using Monte Carlo methods to account for uncertainties in $v\sin i$, radius, and period, we drew 30000 samples from the $v\sin i$ distribution found by Allers et al. (2016) and Gaussian distributions centered at radius = 1.4 ± 0.08 R_Jup (the mean and standard deviation of the radius values found in Allers et al. 2016) and period = 8.61 ± 0.06 hr. The resulting distributions for equatorial velocity (derived from radius and our measured rotation period), $v\sin i$, $\sin i$, and inclination are presented in Figure 7. We adopt a value for $\sin i$ of 0.83 ± 0.07 from the median and standard deviation of our $\sin i$ distribution. Propagating errors in the normal way, we find $i=56\buildrel{\circ}\over{.} 2\pm 7\buildrel{\circ}\over{.} 2$. Values of $\sin i\gt 1$ are unphysical and are a result of our uncertainties in measuring radius and $v\sin i$. We have tried to mitigate this issue in two manners: (1) discarding all values of $\sin i\gt 1$, we find an inclination of 561 ± 74 (median and standard deviation of remaining values), and (2) pinning all $\sin i$ values greater than 1 to 1 (as such a value does imply a high inclination), we find an inclination of $56\buildrel{\circ}\over{.} 2\pm 8\buildrel{\circ}\over{.} 1$. All methods provide consistent values for inclination.

Zoom In Zoom Out Reset image size

We adopted a similar MCMC method to interpret the HST light curves. Since the HST observation does not cover a full period, we fixed the period to 8.6 hr, as determined from the Spitzer light curve. The MCMC code was run for each of the bands from which we extracted light curves; posterior pdfs for each band are presented in the Appendix. Again, we adopt the 50%-quantile value as the best value for each parameter. For each synthesized light curve, a set of 100 samples drawn from the posterior pdf is overplotted on the unbinned HST light curve for that band in Figure 8. We find very little covariance between phase and amplitude in the posterior pdfs, suggesting that the fits are robust, at least for the purpose of estimating amplitudes and phase shifts. Amplitudes (peak to median value) and phases measured for each band from the MCMC fits, as well as phase offsets relative to the Spitzer light curve, are presented in Table 2. At the time of our observations (September 2016), the difference between BMJD and MJD was 7 minutes for HST and 5.6 minutes for Spitzer. In the HST phase offsets relative to Spitzer, we have not accounted for the 1.4 minute offset in timing between the two datasets, as this is considerably less than our 14 minute adopted cadence after binning for HST and is negligible over the 8.6-hour measured period of PSO J318.5-22.

Zoom In Zoom Out Reset image size

Table 2.  Peak-to-median Variability Amplitudes, HST Phases, and Phase Offsets Relative to the Spitzer Light Curve for Synthesized HST Band Light Curves

Band	Amplitude	68%	95%	HST Phase	68%	95%	Phase Offset	68%	95%
White Light	2.51%	0.11%	0.23%	1964	24	48	2069	31	61
(1.07–1.67 μm)
2MASS J	2.92%	0.16%	0.32%	1932	29	57	2037	34	68
2MASS H	2.28%	0.13%	0.25%	1981	29	57	2086	34	68
Water Band	2.38%	0.19%	0.37%	1999	39	77	2104	34	87
(1.34–1.44 μm)
Methane Band	2.20%	0.16%	0.31%	1985	36	72	2090	41	81
(1.60–1.67 μm)

Download table as:  ASCIITypeset image

We fit the median, orbit 1 (close to minimum), orbit 3 (bracketing maximum flux), and orbit 5 (close to minimum) spectra of PSO J318.5−22 with three sets of atmospheric models using very different treatments of cloud parameters: (1) the ExoREM models, which focus specifically on atmospheres with low surface gravity in the cloudy, clear, and partly cloudy cases (Baudino et al. 2015; B. Charnay et al. 2018, in preparation), (2) the BT-Settl models, which explore a wide range of dust species grain formation in the presence of hydrodynamical mixing (Allard et al. 2011), and (3) the thick-cloud models of Madhusudhan et al. (2011; henceforth M11).

5.4.1. ExoREM Model Fits

We compared the grid of ExoREM models to the median spectrum of PSO J318.5−22 through a ${\chi }^{2}$ minimization. Unlike most other 1D models, ExoREM enables the modeling of clear, patchy, and fully cloudy atmospheres, parameterized according to the f_c parameter, where f_c = 0 is a clear atmosphere and f_c = 1 is a fully cloud-covered atmosphere. As a consequence of the fits, we renormalized each synthetic spectrum by a dilution factor ${R}^{2}/{d}^{2}$ that minimizes the ${\chi }^{2}$, where R is the radius and d the distance to the target. The fit was performed independently for each value of f_c, and we compared the best ${\chi }^{2}$ a posteriori. We considered a fit with R left unconstrained, and another one with the radius varying in the interval R = 1.4 ± 0.08 R_Jup (derived from evolutionary model fits assuming an age of 23 ± 3 Myr, Allers et al. 2016). We adopt a distance of 22.2 ± 0.8 pc (parallax of 45.1 ± 1.7 mas) from Liu et al. (2016).

The data are systematically best represented by the ExoREM models with full cloud cover. When R is left unconstrained, we find a best fit for ${T}_{\mathrm{eff}}=1250\,{\rm{K}}$, log g = 3.4 dex, M/H = +0.5 dex, and R = 1.12 R_Jup. When an a priori range on the radius is given, the spectrum of PSO J318.5−22 is best reproduced for T_eff = 1150 K, log g = 3.3 dex, M/H = 0.0 dex, and R = 1.39 R_Jup. We show the solution in Figure 9 along with a ${\chi }^{2}$ map for the solar-metallicity models. The model reproduces the spectral slope and the object near-IR brightness simultaneously, but it fails to reproduce the detailed morphology of the H band and the strength of the water absorption at 1.3–1.5 μm. The variation in spectral features of PSO J318.5−22 are within the error bars of the median spectrum, thus the fitting solutions are identical when our various HST spectra of PSO J318.5−22 are considered.

Zoom In Zoom Out Reset image size

5.4.2. BT-Settl and M11 Model Fits

For the BT-Settl and M11 models, fits were performed with a Levenberg–Marquardt least-squares fitter $({numpy}.{optimize}.{curvefit})$, after binning the model spectra to the same resolution and sampling as the HST spectra. Best-fit models are shown in Figure 10. While both of these models produce spectra that are roughly qualitatively similar to our observed spectra, there are notable differences: the models do not reproduce the steepness of the observed spectral slope from 1.2 to 1.35 μm or from 1.4 to 1.7 μm. As the difference between the model spectra and our observed spectra were larger than the difference between observed spectra from different orbits, we found that a single model spectrum fit best all the single-orbit spectra as well as the median spectrum. In other words, we were unable to describe the changes seen in the observed spectra as a function of time by fitting different model spectra with varying T_eff and surface gravities. The 1D static models used here, however, are essentially averages over the entire visible surface area of the object and over multiple rotational periods, so it is is not surprising that we find similar fits for both single-orbit and median spectra.

Zoom In Zoom Out Reset image size

The best-fit BT-Settl model had T_eff = 1600 K and $\mathrm{log}(g)\,=3.5$, although a range of models with ${T}_{\mathrm{eff}}=1500\mbox{--}1700\,{\rm{K}}$ and $\mathrm{log}(g)=3\mbox{--}5$ fit the spectra nearly as well. The T_eff = 1600 K and $\mathrm{log}(g)=3.5$ best fit is driven by the fitting algorithm's attempt to fit the spectral slope in H band. In contrast, a T_eff = 1700 K, $\mathrm{log}(g)=5.0$ fits the J-band slope better, at the expense of a very poor H-band fit (similar to what was found for HR 8799de by Bonnefoy et al. 2016). These results are consistent with the best fits reported in Liu et al. (2013) for the IRTF SpeX spectrum (${T}_{\mathrm{eff}}=1400\mbox{--}1600\,{\rm{K}}$ and $\mathrm{log}(g)=4.0\mbox{--}4.5\,\mathrm{dex}$). As noted in Liu et al. (2013), this set of atmospheric model spectra fits yield T_eff values that are significantly higher than those obtained from evolutionary model fits to photometry. For the Madhusudhan et al. (2011) models, the best fits were obtained for model A (thick clouds) with 60 μm grains, ${T}_{\mathrm{eff}}\,=1100\mbox{--}1200\,{\rm{K}}$ and $\mathrm{log}(g)=3.75\mbox{--}4.25\,\mathrm{dex}$, consistent with the IRTF SpeX spectrum fits (model A, 60 μm grains,${T}_{\mathrm{eff}}=1100\,{\rm{K}}$ and $\mathrm{log}(g)=4.0\,\mathrm{dex}$) reported in Biller et al. (2015).

As with the ExoREM model fits, we renormalized each synthetic spectrum by a dilution factor ${R}^{2}/{d}^{2}$, which minimizes the ${\chi }^{2}$, where R is the radius and d the distance to the target. We again adopt a distance of 22.2 ± 0.8 pc (parallax of 45.1 ± 1.7 mas) from Liu et al. (2016). Thus, from the dilution factor obtained from the model fit, we estimate the radius R of PSO J318.5−22. This results in an unphysically small radius estimate of ∼0.7–0.8 R_Jup using the BT-Settl models for this young object with low surface gravity. Liu et al. (2013) obtained a similar result fitting this same model set to a low-resolution near-IR spectrum; in contrast, adopting an age of 23 ± 3 Myr, Allers et al. (2016) find T_eff = 1100–1200 K and radius values of 1.34–1.46 R_Jup using a variety of evolutionary model grids with and without clouds. The observed luminosity of a (very roughly blackbody) object is governed by the temperature and the radius; the high-temperature fit by the BT-Settl models has thus necessitated an unphysically small fit to the radius to produce the observed luminosity. For the M11 models, with a cooler fitted temperature of 1100–1200 K, we find radius estimates of ∼1.0–1.3 R_Jup, roughly consistent with our ExoREM model fits with radius left as a free parameter, but still somewhat smaller than estimates based on evolutionary models (Liu et al. 2013; Allers et al. 2016).

Over our 7 hr of HST monitoring, we only captured one clear extremum, a maximum in brightness that occurs in orbit 3. The minimum value of brightness measured during our time series occurred in orbit 1. However, orbit 1 was the most affected by the ramp effect. Orbit 5 is also near a minimum of the light curve and should not be affected as strongly by the ramp effect. The ratios of maximum and minimum spectra (taking both orbits 1 and 3 as potential minima) are plotted in Figure 11. The orbit 3 spectrum divided by the orbit 5 spectrum shows a monotonic and relatively small decrease in this ratio as a function of wavelength. In contrast, for the orbit 3 spectrum divided by the orbit 1 spectrum, the ratio of maximum to minimum spectral flux in the 1.4 μm water absorption feature is slightly smaller relative to that at adjacent shorter and longer wavelengths, superimposed on a monotonic small decrease in variability amplitude as a function of increasing wavelength. The suppression of variability in the water band relative to the adjacent continuum in orbit 1 appears to be robust—from Figure 8, all light curves except the water-band light curve display a significant deviation from the sinusoidal MCMC fits during orbit 1. This is likely not a result of the ramp effect, which should affect all wavelengths equally. For the L/T transition objects SIMP 0136 and 2M 2139, Apai et al. (2013) and Yang et al. (2016) found significantly smaller amplitudes in the water absorption feature relative to both the shorter and longer wavelength continuum, while for two mid-L dwarfs, Yang et al. (2015) found small but monotonic decreases of amplitude with increasing wavelength. Our results appear to be a hybrid of these cases, with different behavior observed in different orbits.

Zoom In Zoom Out Reset image size

As our HST observations did not cover a full period, we did not measure the full amplitude of variability. However, we did cover the majority of the period and have a robust period determination from the simultaneous Spitzer observations, so we can use the sinusoidal fits to the light curves to estimate amplitudes and phase shifts (relative to the Spitzer light curve) for the five broadband regions we have considered previously: the full white-light spectrum from 1.1 to 1.67 μm, 2MASS J, 2MASS H, water, and methane. As noted previously, this method appears to slightly underestimate the variability amplitude, as most of the synthesized light curves show some deviation from the sinusoidal fits during orbit 1. The measured amplitude and phase shift for each of the broadband regions we considered are plotted as a function of wavelength in Figure 12. Similar to the divided spectra for orbits 3 and 5, the amplitude appears to generally decrease as a function of increasing wavelength, with the sharpest break between J and the 1.4 μm water band.

Zoom In Zoom Out Reset image size

The mid-IR Spitzer [4.5 μm] light curve follows the same trend of decreasing amplitude with longer wavelength, with a peak-to-trough amplitude of ∼3.4% versus 4.4%–5.8% for the near-IR bands. For field brown dwarfs, the near-IR variability is generally found to have a significantly higher amplitude than the mid-IR variability when both are present. It is notable that the mid-IR variability amplitude for PSO J318.5−22 is so similar to its near-IR variability amplitude and is in fact one of the highest variability amplitudes ever measured in the mid-IR for a brown dwarf or planetary-mass object! The highest amplitude variable from Metchev et al. (2015) is the T6 dwarf 2MASS J22282889−4310262, with a 3.6 μm variability amplitude of 4.6% ± 0.2%, but most of the variables in their sample have amplitudes of <2%. This high-amplitude variability observed for PSO J318.5-22 in the mid-IR may be the effect of low surface gravity on the vertical structure of such an atmosphere. Low surface gravity allows cloud species to potentially extend up to lower pressures and higher altitudes compared to the high surface gravity case (Marley et al. 2012). In general, for brown dwarfs and free-floating planetary-mass objects, the photosphere in the mid-IR is at lower pressures and higher altitudes than the photosphere in the near-IR (see, e.g., Marley et al. 2012; Biller et al. 2013a). Thus, the extension of clouds up to lower pressure regions increases the chance of heterogeneous cloud opacity (and hence variability) at the low-pressure levels probed by mid-IR observations.

The viewing inclination of a given object will significantly affect its variability properties (see, e.g., Vos et al. 2017a). We presume that the surface features that generate variability are primarily equatorial, thus the same object observed at high inclination will appear to have a higher variability amplitude than if viewed at lower inclination. Additionally, Vos et al. (2017a) find that J-band variability is more affected by inclination angle than mid-IR variability. As J-band observations generally probe a deeper part of the atmosphere than mid-IR observations (Biller et al. 2013b; Yang et al. 2016; Vos et al. 2017a), Vos et al. (2017a) propose that this effect may be due to the increased atmospheric path length of J-band flux at lower inclinations. We measure a relatively high inclination for PSO J318.5−22 of 56° ± 8°, thus, we are observing close to the full amplitude in each band.

We measure a phase shift between the Spitzer light curve and HST light curves of ∼200°. At a lower significance, we find a ∼6°–7° phase shift between the J-band light curve and the other near-IR narrowband light curves. Similar results have been found for old, high-surface gravity brown dwarfs phase shifts between the near-IR and mid-IR light curves may be quite common for these objects. The first such phase shift was found by Buenzli et al. (2012) for the T6 object 2MASS J22282889−4310262, and Yang et al. (2016) recently found significant phase shifts between the mid- and near-IR for four brown dwarfs with high surface gravity (including 2MASS 2228). In three out of four cases from Yang et al. (2016), the measured phase shift between near-IR and mid-IR is nearly ∼180°, ranging from 150° to 210° (see their Figures 19–22), very similar to the ∼200° phase shift reported here for PSO J318.5−22. Phase shifts within the near-IR spectral bands are less common (Buenzli et al. 2012; Biller et al. 2013a; Yang et al. 2016), but have been reported now at multiple epochs for the T6 object 2MASS J22282889−4310262 (Buenzli et al. 2012; Yang et al. 2016). For PSO J318.5−22, all the near-IR bands we consider agree in phase at the 2σ level; at the 1σ level, 2MASS J is shifted by ∼6 degrees relative to the other near-IR bands. However, given that we have monitored less than one rotation period with the HST and have also assumed a sinusoidal light-curve shape, this phase shift is still within the errors expected from our sinusoidal model fitting. Observed phase shifts have generally been interpreted as different "top-of-atmosphere" locations at different wavelengths—near-IR generally probes deeper in the atmosphere than mid-IR (Buenzli et al. 2012; Biller et al. 2013a; Yang et al. 2016). In other words, the source of the inhomogeneity that drives the near-IR variability is located at a higher pressure level deeper in the atmosphere than the source of inhomogeneity driving mid-IR variability at a lower pressure level.

Following the method of Apai et al. (2013), we performed a principal component analysis (PCA) to determine how many spectral components drive the variability for PSO J318.5−22. We expect the observed variability to be driven by the rotation in and out of view of regions with differing spectra. The PCA identifies the smallest set of independent spectra that account for the majority of the observed variability. Taking each of the 49 spectra in our spectral sequence as a dimension and subtracting the mean value for each spectrum (as required for a PCA), we calculate the 49 × 49 covariance matrix between each spectrum using the ${numpy}.{cov}$ function in python. We then used ${numpy}.{linalg}.{eig}$ to determine eigenspectra and eigenvalues. Sorting on eigenvalues, the principal spectral component accounted for 99.6% of the observed variability, with the second component contributing 0.1% and the third component contributing 0.07%. This is similar to results for other variable objects—Apai et al. (2013) find that the principal spectral component accounted for 99.6% and 99.7% of the variability for 2M2139 and SIMP 0136, respectively. They argue that this implies that only two types of surface patches are required to explain the observed variability for these objects. However, subtracting the mean essentially means that we remove any gray variation—any shifting of the entire spectrum by a constant value. Figures 11 and 12 show that the observed variability does not possess a strong color component, as variability amplitude changes monotonically and slowly with wavelength over the 1.07–1.67 μm spectral range of the WFC3 G141 grism—for instance, the variability amplitude in the methane band is roughly 75% of that in the J band. Thus, removing the mean for the PCA actually removes most of the observed variability. Hence it is not surprising that the principal spectral component (which closely resembles the median spectrum) encompasses the vast majority of the power in the time series.

The key variables for constraining variability properties in both field brown dwarfs and young exoplanets are spectral type (which correlates with T_eff) and surface gravity. High-temperature objects (spectral types from early-L to mid-L) may have variability from cloud features as well as magnetic activity (e.g., star spots and aurorae); lower temperature objects (>L5) are assumed to have atmospheres too cool to produce magnetic activity and are presumed to have entirely cloud-based variability (Gelino et al. 2002). The variability of older field dwarfs with high surface gravity across the full L–T spectral has been studied in detail spectroscopically (Buenzli et al. 2012; Apai et al. 2013; Buenzli et al. 2015a, 2015b; Yang et al. 2015). However, only three young, low surface gravity have HST spectral variability monitoring to date—SIMP0136 (Apai et al. 2013), W0047 (Lew et al. 2016), and the observations presented here for PSO J318.5−22. Here we consider a number of drivers of variability suggested in the literature for objects with both high and low surface gravity (specifically, patchy cloud features, hot spots, high-level hazes, thermochemical instabilities, and magnetically driven aurorae), to determine if they can describe the variability properties observed for PSO J318.5−22.

Considering patchy salt and sulfide clouds as well as hot-spot models (heating at a specific pressure level) for objects with high surface gravity and ${T}_{\mathrm{eff}}\gt 375\,{\rm{K}}$, Morley et al. (2014) find that variability due to patchy clouds should drive high-amplitude variability across a wide spectral range, while variability due to hot spots should produce larger variability amplitudes within absorption features relative to continuum wavelengths. As noted in Section 6.2, for PSO J318.5−22, we find a smooth decrease of variability amplitude as a function of increasing wavelength across the 1.07–1.67 μm spectral range. For the 1.4 μm water absorption feature and the 1.6 μm methane absorption feature covered by the HST WFC3 grism, we find similar or slightly smaller variability amplitudes relative to the adjacent continuum. The lack of stronger variability within absorption features compared to continuum wavelengths for PSO J318.5−22 implies that inhomogeneous cloud features (thick and thin clouds, or a haze layer over a thick cloud surface) is likely to be a major driver of the observed variability.

The 1.4 μm water absorption feature can potentially provide a useful diagnostic of cloud-driven variability mechanisms. This feature is available from space with HST but is very difficult to observe from the ground. In the case of variability from high-level hazes (as observed for two mid-L field dwarfs by Yang et al. 2015), we expect similar variability amplitudes both with and within the water absorption feature, if the high-level hazes driving the variability are located at a very low pressure high altitude in the atmosphere where the water opacity is negligible. In the case of variability due to inhomogeneous thin and thick clouds, we expect the variability amplitude to be notably different in the water absorption feature relative to adjacent non-absorbed wavelengths. For instance, Apai et al. (2013) found variability to be suppressed at 1.4 μm relative to the adjacent continuum for two highly variable L/T transition brown dwarfs. However, as noted in Section 6.2 and Figure 11, PSO J318.5−22 displays both behaviors in different orbits! Thus, in this case, the amplitude inside and outside the 1.4 μm water absorption feature does not clearly distinguish between the cases of high-level hazes versus inhomogeneous thin and thick clouds.

However, recent work has suggested that clouds may not be necessary to model the spectra (Tremblin et al. 2016, 2017)—or variability—of objects of L and T spectral type. Tremblin et al. (2016) have recently produced cloud-free models of L- and T-type brown dwarf atmospheres, successfully modeling the red colors of L dwarfs as well as the T dwarf J-band brightening and reemergence of the FeH absorption feature using additional convection from thermochemical instabilities (in the CO/CH₄ transition in the case of the L/T boundary). They suggest that turbulence produced by CO or temperature fluctuations across the CO/CH₄ boundary may be a driver of the observed variability of brown dwarfs, and in particular, that the inhomogeneous top-of-atmosphere structure mapped via Doppler imaging for the L/T transition brown dwarf Luhman 16B (Crossfield et al. 2014) may be explained by inhomogeneities in CO versus CH₄ abundance or temperature. We consider whether this mechanism might drive variability for PSO J318.5−22. Variability due to abundance variations should drive increased variability amplitudes in the absorption features produced by the species in question. For this reason, we produced synthesized light curves in the 1.6 μm methane absorption feature—at least the portion of it that lies within the 1.07–1.67 μm spectral range of the HST WFC3 G141 grism. We do not find the variability amplitude in the methane band to be significantly enhanced or suppressed relative to the wider 2MASS H band (see Figure 12). We tentatively suggest that the observed variability is not driven by varying CH₄ abundance here—however, a full theoretical calculation of the expected variability amplitude in methane absorption features due to thermochemical instabilities is not yet available, thus, we await more quantitative theoretical predictions here.

While objects with lower temperature (>L5) are commonly assumed to lack magnetically driven variability (Gelino et al. 2002), Hallinan et al. (2015) suggest that this may not always be the case. Hallinan et al. (2015) find significant phase shifts between radio and various optical bands for the much hotter nearby M8.5 object LSR J1835+3259, which they interpret as auroral heating. In particular, electron beams from global auroral current systems feed energy from the magnetosphere into the atmosphere of this object. Hallinan et al. (2015) suggest that this mechanism may extend down even to very cool brown dwarfs, driving some of the more extreme examples of weather phenomena in brown dwarfs. They also note that radio emission has been detected from objects with spectral types as late as T6.5. From our current data set, we cannot confirm or refute if this is the case for PSO J318.5−22; radio and H-α monitoring would be necessary to do so. However, more generally, Miles-Páez et al. (2017) searched for H-α emission for a sample of eight L3-T2 field brown dwarfs, six of which have detections of photometric variability. The only H-α detection in this sample was from a non-variable T2 dwarf, suggesting that aurorae and other chromospheric activity do not commonly drive variability for L and T spectral type objects.

Assuming that variability may be cloud driven, what cloud species and cloud geometries are necessary to reproduce our observed amplitudes and phase shifts across the near- and mid-IR? To try to quantify the expected pressure level at which the photosphere is found at a given wavelength for PSO J318.5−22, we considered the best-fit model to our HST median spectrum to identify at what pressure level flux is being emitted at each wavelength. As our BT-Settl model fit in Section 5.4.2 yielded a higher temperature and smaller radius than is consistent with evolutionary model fits to the same object (Liu et al. 2013; Allers et al. 2016), we consider only our model fits using the M11 models and the ExoREM models for this analysis. We adopt the ExoREM fully cloudy model with ${{\rm{T}}}_{\mathrm{eff}}=1150\,{\rm{K}}$, log g = 3.3 dex, M/H = 0.0 dex, and R = 1.39 R_Jup. In Biller et al. (2015), we found that the SpeX spectrum for PSO J318.5−22 is best fit with the A60, 1100 K, solar metallicity, log(g) = 4 model from Madhusudhan et al. (2011) and found similar fits for our HST time-resolved spectra (see Section 5.4). As the SpeX spectrum fit is similar to what we find here and covers a wider wavelength range, we adopt this model fit for the M11 models. Photospheric pressures are provided with the publicly available M11 models; for the ExoREM models, a pressure spectrum was generated by combining the model spectrum with the pressure/temperature profile of the model. The flux at each wavelength was converted into the equivalent brightness temperature. We then interpolated using the corresponding pressure/temperature profile to obtain the photospheric pressure level. This method is correct if the source function varies linearly with optical depth. The resulting "pressure spectra" for both of these models are plotted in Figure 13, with the bandwidth for each of our light curves overplotted. While the photospheric pressure level versus wavelength varies between models, in both cases, the mid-IR flux is generated higher in the atmosphere than the near-IR flux.

Zoom In Zoom Out Reset image size

What cloud species are expected to dominate at the respective higher and lower pressures probed by near-IR versus mid-IR observations? In Figure 14, we plot the pressure/temperature profile for both our best-fit ExoREM cloudy model and an equivalent clear model with other parameters unchanged. Thick lines correspond to the photosphere (computed from 0.6 to 5 μm), and dashed lines are condensation temperatures for the different clouds present in the model. The presence of clouds (red curve) increases the temperature by around 200 K. In the cloudy case, the type of cloud that condenses varies according to the pressure/temperature profile and the condensation temperatures for different cloud species. Silicate and iron clouds form at around 1 bar and are optically thick up to around 0.3 bar at 1 μm. Thus, for these model atmospheres, silicate and iron clouds form below the photosphere pressures of the 1.4 μm water band and the 4.5 μm CO band. In contrast, Na₂S clouds form in the upper atmosphere at around 0.06 bar, so above the photosphere pressures of all molecular bands except the 4.5 μm CO band.

Zoom In Zoom Out Reset image size

Longitudinal variations in cloud thickness can potentially produce anticorrelated variability, in other words, ∼180° phase shifts between the near-IR and mid-IR light curves. In Figure 15 we plot wavelength versus brightness temperature, showing where in the spectrum the brightness temperature increases/decreases with clouds. In the cloudy case, the brightness temperature increases at longer wavelengths (e.g., 3 and 4.5 μm) by around 200 K relative to the clear case. The opposite is true at shorter wavelengths (∼1–2 μm), where the brightness temperature decreases by 200 K relative to the clear case. A hole in the cloud cover would thus produce a 180° phase shift between near-IR and mid-IR light curves, except for the 1.4 μm water band, which would be correlated with the mid-IR light curve, in contrast to the observations presented here for PSO J318.5−22. However, we do not expect any fully clear patches on this object (and indeed previous work suggests that this is the case for brown dwarfs in general, Apai et al. 2013), but rather longitudinal variations in the cloud thickness. The phase shift and the amplitude of light curves are then very dependent on the altitude, thickness, and placement of different cloud species.

Zoom In Zoom Out Reset image size

We consider a number of simple geometries for both silicate/iron clouds and sulfide clouds to model our observed amplitudes and phases. In Figure 16 we compute the light-curve amplitude that would be produced assuming (1) a spot with optically thinner silicate and iron cloud thickness, covering 10% of the surface, with homogeneous thick silicate and iron clouds over the rest of the surface; and (2) one hemisphere covered by high-altitude sulfide clouds and no sulfide clouds on the other hemisphere, with homogeneous silicate/iron clouds below the sulfide cloud layer altitude for both hemispheres. Case (2a) was computed assuming no horizontal heat redistribution between the less cloudy spot and the rest of the brown dwarf (solid line). Case (2b) was computed with very efficient heat redistribution (dashed line). For longitudinal variations (Case 1) in silicate and iron cloud thickness, we predict large variations in the amplitude within the 1.07–1.67 μm spectral range of the HST WFC3 G141 grism and, additionally, the 1.4 μm water-band light curve should be correlated with the 4.5 μm Spitzer Channel 2 light curve. For longitudinal variations in the sulfide cloud cover (Case 2a and b), the amplitude is quite constant in the HST bands, and for a case intermediate between efficient heat redistribution and no heat redistribution, the predicted amplitude and the phase shifts between different near-IR and mid-IR wavelengths could be compatible with our HST and Spitzer observations. This modeling remains very preliminary, but suggests that variations in the cover of high-altitude clouds could begin to explain the observations presented here. Na₂S clouds are a good candidate for such high-altitude clouds since they form at very low pressures high in the atmosphere (0.06 bar). Inhomogeneous Cr and MnS clouds also are potential candidates. An upper layer of silicate clouds could also produce this variability, but it would require a mechanism for forming or transporting cloud particles higher than the cloud deck. Cloud convection triggered by latent heat release (Tan & Showman 2017) or radiative heating (Freytag et al. 2010) may produce vertically extended clouds and a detached silicate haze layer.

Zoom In Zoom Out Reset image size

Our measured phase shifts between the HST bands and the Spitzer 4.5 μm light curve are in fact somewhat more than 180°, which is unsurprising, as the cloud geometry for PSO J318.5−22 is certainly more complicated than the simple geometries considered above. Modeling approaches that combine multiple 1D models (such as the one presented herein and Artigau et al. 2009; Apai et al. 2013; Karalidi et al. 2015, 2016) can reproduce correlated variability or 180° anticorrelated variability, but not other phase shifts. As demonstrated above, where the "top-of-atmosphere" occurs varies depending on wavelength and the specific opacity sources that dominate at different atmospheric levels. The observed "phase shifts" may simply be heterogeneous and uncorrelated structure at different altitudes, which is still modulated by the rotation period of the object in question. Most likely, full 3D models will be necessary to describe this structure, such as the C05BOLD model currently undergoing testing (F. Allard et al. 2018, in preparation), especially as rotation probably plays a significant role in the appearance and features of these atmospheres (Showman & Kaspi 2013).

Metchev et al. (2015) find increased mid-IR variability amplitudes for eight L3–L5 objects with low surface gravity with respect to the rest of their older survey sample with high surface gravity. We tentatively find that such a trend (in both the near- and mid-IR) may continue for mid- to late-L dwarfs with low surface gravity. Only four such objects have been surveyed in either the near- or mid-IR to date. Three out of the four have positive variability detections in both near- and mid-IR (Morales-Calderón et al. 2006; Biller et al. 2015; Lew et al. 2016; Vos et al. 2018); one is a non-detection in our ongoing SofI survey (Vos et al. 2018, in preparation). The three mid- to late-L dwarfs with low surface gravity and positive variability detections are PSO J318.5−22, W0047, and 2M2244. For our HST+Spitzer monitoring of PSO J318.5−22, we found peak-to-trough amplitudes of ∼3.4% for Spitzer Channel 2 and 4.4%–5.8% in the near-IR band (1.07–1.67 μm) covered by the WFC3 G141 prism. In the discovery epoch, Biller et al. (2015) found peak-to-trough variability amplitudes of 7%–10% in the J_S band and ∼3% in K_S, indicating evolution of the variability between the discovery epoch and our HST observations. The lower amplitude in K versus J during the discovery epoch is consistent with our finding in this work that the variability amplitude decreases with increasing wavelength across the 1.1–1.7 μm spectral range of the HST WFC3 grism. Lew et al. (2016) find a similarly high near-IR variability amplitude for W0047, with the relative variability amplitude decreasing from 11% at 1.1 μm to 6.5% at 1.7 μm. Vos et al. (2018) reported a mid-IR detection for this object with a relative variability amplitude of 1.07% ± 0.04%. Morales-Calderón et al. (2006) measured a Spitzer Channel 1 peak-to-peak variability amplitude of 8 mmag for 2M2244, an L6.5 AB Dor member (Vos et al. 2018). Variability in this object has recently been confirmed by Vos et al. (2018), who found a Spitzer Channel 1 peak-to-peak variability amplitude of 0.8 ± 0.2% and $\geqslant 3 \%$ variability amplitude in J band in a 4 hr long J-band UKIRT WFCAM observation of this object. All three of these objects have notably high near-IR amplitudes compared to field brown dwarfs with similar spectral types as well as planetary-mass objects with earlier spectral types. For instance, the detection of variability in the L5 planetary-mass object 2M1207b has a considerably lower near-IR amplitude of 1%–2% (Zhou et al. 2016). Of the three, PSO J318.5 also has a notably high mid-IR amplitude; mid-IR amplitudes for the other two objects are more in line with typical values for field brown dwarfs.

With such a small number of low surface gravity mid- to late-L variables to study, it is not clear whether these three objects are unusual or if objects with low surface gravity are inherently more variable then their counterparts with high surface gravity. The variability peak for field brown dwarfs appears to be at the L/T transition (Radigan et al. 2014; Radigan 2014), commonly attributed to the breakup or at least thinning of silicate clouds at this spectral type transition (Apai et al. 2013). It is hard to say if this is the case for objects with low surface gravity, with three high-amplitude near-IR detections for mid- to late-L objects with low surface gravity (Biller et al. 2015; Zhou et al. 2016; Vos et al. 2018), one high-amplitude detection in a young T2.5 object (Artigau et al. 2009; Gagné et al. 2017), and one tentative detection in a young T3.5 object (Naud et al. 2017). Predominantly early-L objects with low surface gravity have been surveyed to date (Vos et al. 2018, in preparation), largely because of the current scarcity of late-L, L/T transition, and T spectral type young objects with low surface gravity. Nonetheless, the few mid- to late-L objects surveyed to date appear to be notably variable, which is surprising given that late-L objects are expected to have thick (and probably homogeneous) cloud cover. If late-L spectral type young objects are as a class highly variable, this may draw into question the interpretation of high-amplitude variability as the breakup of silicate clouds between the L and T spectral type.

Mid- to late-L dwarfs with low surface gravity are particularly interesting because these objects are excellent proxies for several known giant exoplanet companions. The spectra of PSO J318.5−22 and W0047 are nearly identical to those of the inner two HR 8799 planets (Bonnefoy et al. 2016). De Rosa et al. (2016) find that the spectrum of the particularly red planet HIP 95086b (Rameau et al. 2013) closely matches that of 2M2244. The newly discovered exoplanet companion HIP65426b also has an L5–L7 spectral type (Chauvin et al. 2017). Given the significant variability of PSO J318.5−22, W0047, and 2M2244, we may expect exoplanet companions such as HR 8799bcde, HIP 95086b, and HIP 65426b to be similarly variable, although the viewing angle (likely pole-on for the HR 8799 system) may render that variability hard to detect.

Even if young planetary-mass objects have significant top-of-atmosphere inhomogeneities, we will only be able to detect such features if these objects are relatively rapid rotators (periods <20 hr). Many old, field brown dwarfs are rapid rotators (Zapatero Osorio et al. 2006). From conservation of angular momentum, one might expect young objects to be predominantly slower rotators compared to old, field brown dwarfs, as they have somewhat inflated radii (e.g., ∼1.4 R_Jup for PSO J318.5−22 Allers et al. 2016) compared to older objects (radius ∼1 R_Jup) and will be expected to spin up with age as they contract. At least preliminarily, however, there is a small cohort of young ($\leqslant 150\,\mathrm{Myr}$), planetary-mass objects with periods <20 hr, including PSO J318.5−22, as well as the bona fide exoplanet β Pic b and 2M1207b. In Figure 17 we plot estimated object mass versus measured equatorial velocity for these objects, solar system objects, and field brown dwarfs with measured periods from Vos et al. (2017a). Planetary-mass objects seem to encompass a similar range of equatorial velocities as older, field brown dwarfs, with both rapid rotators and notable slow rotators such as the young, 30–40 M_Jup brown dwarf companion GQ Lup b (Schwarz et al. 2016). However, statistics are still too sparse for a robust comparison to the brown dwarf population in general. Preliminary analyses do suggest that the rotation rate between free-floating and companion objects is similar: Bryan et al. (2017) recently measured the rotation rate for a number of companions with masses <20 M_Jup. Combining their measurements with others in the literature, they found no discernable difference in rotation speed between companions and free-floating objects with similar masses for a small sample of 11 objects.

Zoom In Zoom Out Reset image size