BROWN DWARF PHOTOSPHERES ARE PATCHY: A HUBBLE SPACE TELESCOPE NEAR-INFRARED SPECTROSCOPIC SURVEY FINDS FREQUENT LOW-LEVEL VARIABILITY

The initial target sample was selected from the Dwarf Archives,⁷ a compendium of over 1000 brown dwarfs, to uniformly span the spectral subtypes between L5 and T6, sampling effective temperatures between ≈1700–800 K and very different stages of cloud evolution. No special selection for or against known young low-gravity objects were made. The selection was unbiased in the sense that prior knowledge about variability was not a selection criterion except for the exclusion of the three known variable T dwarfs previously observed with HST. We excluded objects that were known resolved binaries (two known binaries were mistakenly included, one of them in the observed sample) and objects with known Two Micron All Sky Survey (2MASS) sources within 20''  in order to minimize the risk for overlapping spectra. We prioritized the targets by their J-band brightness, but ensured that the selected sources were evenly distributed in color–magnitude space and spectral type. Because of the random selection of the subset of targets that were actually observed from this sample, the final spectral type distribution of the observed sample is not entirely uniform.

We binned the subset of 22 targets that were observed in the SNAP survey into four spectral bins that correspond to distinct evolutionary stages: mid-L (L5–L7), late-L (L7.5–L9.5), early-T (T0–T3), and mid-T (T3.5–T6). The early-T dwarfs, which host the strongest known near-IR variables (Artigau et al. 2009; Radigan et al. 2012; Apai et al. 2013; Gillon et al. 2013), are the smallest sample with only three objects. The mid-L and mid-T sample has seven objects each, the late-Ls five objects. The target properties are summarized in Table 1.

Observations were taken between 2011 October and 2012 October with HST/WFC3 in the IR channel. We used the G141 grism which provides slitless spectra for wavelengths 1.1–1.7 μm. The detector is a Teledyne HgCdTe with 1024 × 1024 pixels. The pixel size is 013. All observations were executed in the 256 × 256 subarray mode with a field of view of approximately 30'' × 30''. In this mode, for our exposure times all images acquired in a visit could be stored without the need for a WFC3 buffer dump, ensuring maximum observing efficiency. The first order of the spectrum is fully captured on the subarray, while the zeroth and second orders are not recorded. The spectrum has a dispersion of 4.65 nm pixel⁻¹ and spans ≈140 pixels. For wavelength calibration, we obtained a direct image at the beginning of the visit through the F132N or F127N narrowband filter that provides an accurate measurement of the location of the source on the detector.

We used the SPARS25 readout mode for all targets. In this mode, a sequence of non-destructive reads is taken in one exposure. At the beginning of each exposure, a 0 s read and a 0.27 s read are taken. Then, a number of reads of 22.34 s are obtained. The time for each read is fixed, but the number of reads can be chosen between 1 and 15 depending on the brightness of the target. We set the number of reads to 1, 2, 5, or 10 such that the maximum number of counts acquired by the brightest pixel was between 5000 and 15,000 counts. This ensured that there was no significant image persistence, which can become relevant if a pixel is exposed above about half-well (≈40,000 counts). The number of exposures varied between 9 and 67 depending on the number of reads (effective exposure time) and duration of the visit. The spectroscopic time series lasted between 32 and 45 minutes and the cadence (exposure time plus overhead) was between 41 and 242 s. The observation details for each source are summarized in Table 2.

We kept the location of the spectra on the detector fixed and did not dither during the observations. This avoided issues with pixel-to-pixel sensitivity variations that cannot be corrected to sub-percent precision by flat fielding. The positional stability was better than 0.06 pixels for all targets as verified by cross-correlation of the images.

For four targets, parts of the spectrum could not be used: for 2M0310, the proper motion given in Faherty et al. (2009) was found to be incorrect.⁸ Due to the resulting location of the spectrum near the upper edge of the subarray, wavelengths from 1.18 to 1.26 μm and >1.65 μm were affected by a large number of bad pixels and had to be cut from the spectrum. For 2M0817, a mistake with the proper motion value in the observing preparation file resulted in the spectrum being located across the right edge of the subarray, and wavelengths >1.61 μm were not read out. For target 2M1624, the first- and second-order spectra of a background star that was missed in the target selection overlap partially with the target spectrum for λ < 1.18 μm and λ > 1.65 μm. For target 2M1711, the spectrum of a very faint background star which is not in the 2MASS Point Source Catalogue partially overlaps for λ ≳ 1.65 μm.

We reduced the data that is output by the standard WFC3 pipeline with custom IDL routines and the PyRAF software package aXe,⁹ a tool developed for extracting and calibrating slitless spectroscopic data. The WFC3 pipeline calfw3 delivers two-dimensional spectral images that are zero-read and dark subtracted, corrected for non-linearity and gain and include flags for bad pixels.

We started from the .flt files that already include the combined images from all subreads of an exposure. Unlike in Apai et al. (2013), we do not extract the individual subreads because the signal-to-noise ratio (S/N) of the spectra in the individual subreads is too low. We therefore do not apply the correction for the small flux loss within one exposure as a function of subreads. This correction would slightly increase the average flux of each exposure by a common factor and would not introduce any relative change in the time series.

Cosmic ray hits were identified as >5σ outliers for a given pixel compared to the same pixel in the nearest eight frames in the time series and replaced by the median. We corrected bad pixels that were flagged by the pipeline by interpolating over nearest neighbors in the same row. Like Berta et al. (2012) we found that only flag numbers 4 (dead pixel), 32 (unstable pixel), 256 (saturated), and 512 (bad in flat field) impact the flux in a pixel in a significant way, we therefore only corrected bad pixels with those flags.

Because the aXe pipeline cannot handle subarray images, we embedded the frames into full-frame images and flagged the extra pixels with a data quality flag to exclude them from further processing. With the axeprep routine we subtracted the background that is determined by scaling a master-sky frame. We then used the axecore routine that flat fields the frames, performs wavelength calibration, extracts the two-dimensional spectra, and flux-calibrates with the G141 sensitivity curve. We chose the extraction width to be 6 FWHM. For smaller values, we found that for some objects an artificial variability slope could be introduced that increased for narrower extraction widths. We did not find significant differences when using 6–8 FWHM, therefore we used 6, which minimizes the noise. The FWHM is measured on the direct image at the beginning of the observations.

We calculated the error in each pixel as a combination of the photon noise of the source, the error from the sky subtraction and the readout noise. In the final step we corrected a ramp effect that we previously identified and discussed in Apai et al. (2013). There, we corrected the effect by using an analytical function fitted to data from a non-variable star. Because our current program includes several additional non-variable sources, we refine the ramp correction as elaborated in the next section.

Time series observations with the WFC3/IR channel have shown the presence of a ramp effect where the measured flux increases strongly at the beginning of an orbit and then flattens out. The effect is present for several readout modes, but characteristics differ. The ramp was first analyzed for the RAPID readout mode by Berta et al. (2012), and characterized in more detail for different subarray sizes and sampling modes in RAPID and SPARS10 readout mode by Swain et al. (2013) and Mandell et al. (2013). They find a correlation of the ramp features with the length of the WFC buffer dump. Here we discuss the ramp for time series data taken in the SPARS25 mode and the 256 × 256 subarray, which allows observations of faint targets without intra-orbit buffer dumps.

The flux increase corresponds to ≈2% during the first orbit of a visit, and only ≈0.7% for subsequent orbits in the same visit (Apai et al. 2013). Because all data in this program were taken in single orbits, we focus here only on the ramp characteristics of the first orbit.

In a first step, we corrected the data from our targets with the ramp derived from the non-variable star data from program GO12314 (Apai et al. 2013). There, we had already found that the ramp depends neither on wavelength nor count rate. Figure 1 shows the ramp of the star for different wavelength regions. It agrees very well for all wavelengths. Berta et al. (2012) also found their ramp to be achromatic. We therefore derived the ramp for each target by integrating over the full spectral range from 1.12 to 1.66 μm to maximize the S/N. We excluded the four sources with missing wavelength segments or overlap with a background star in order to avoid other systematics. It was immediately evident that six of the remaining 18 objects showed clearly different behavior than the non-variable star, and different for each object. This indicates inherent time-variability for these objects on top of the detector systematic. We therefore excluded these objects from further analysis for the ramp. From the other 12 objects and the non-variable star, we calculated a moving average with time steps of 30 s and bin width of 180 s. Subtracting this ramp from the data, we found small trends for three additional sources, again with different trends for each source. We therefore removed those as well for the final ramp correction. For the remaining nine objects (2M0243, 2M0421, 2M0539, 2M0801, 2M0908, 2M0909, 2M1039, 2M1515, and 2M1632), the agreement of the ramp with the non-variable star was very good for all times after ≈180 s after the beginning of the observations (Figure 2). There is a large scatter in the first 3 minutes where the ramp is steepest. We therefore disregard the first 180 s of each time series from all further analysis in order not to introduce artificial trends at the beginning of the time series. Beyond 43 minutes, only a few objects have data. From there, we extend the ramp as a horizontal line. Comparing with the previous 10 minutes, this is valid to ≈0.1% level.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The final ramp is shown in Figure 2. We find evidence for small fluctuations that deviate from the analytical fit made in Apai et al. (2013). If we create the ramp by averaging different subsets of these objects, the differences are <0.1%. After dividing the data by the derived average ramp, we find a standard deviation of 0.115% for the residuals. Because of the small-scale fluctuations, there may be correlated errors on the order of 0.1% over the span of 5–10 minutes. However, the overall distribution of residuals is reasonably close to Gaussian (Figure 3) with the standard deviation similar to the average random error of individual points (0.090% ± 0.015%). We add an error of 0.1% in quadrature to account for the uncertainties in the ramp correction. For points shortward of 5 minutes and longward of 39 minutes, where the ramp correction is less certain, we add an error of 0.2% in quadrature.

Zoom In Zoom Out Reset image size

We tested whether the ramp or the residuals after the correction for the different objects correlate with the positional shifts of 0.01–0.06 pixels perpendicular to the spectral trace within one orbit. We do not find any correlation and conclude that the ramp and the trends after the ramp correction are not caused by positional instabilities as it is the case for warm Spitzer photometry (e.g., Heinze et al. 2013; Lewis et al. 2013). We note, however, that for some objects, small shifts (≲ 0.1 pixels) parallel to the spectral trace can result in small wavelength shifts. This can produce an artificial variability trend if the spectra are integrated over a narrow wavelength interval with a strong gradient in counts, in particular at >1.66 μm where the grism sensitivity drops strongly. We cannot reliably fit for the shifts because their influence on the other parts of the spectrum is too small. Similarly, Mandell et al. (2013) used only the spectral edges to fit for these shifts in their data of transiting planets. We therefore cannot reliably disentangle the effect of the shifts from potential true spectral variability that occurs only in that region. We exclude these outermost wavelength regions from our analysis and check for all detected variables that wavelength shifts cannot be responsible for the observed trends.

Finally, we used the data of a background star from a newer HST program (13176, PI: Apai) that used the same observing mode to verify our reduction and calibration. We reduced and analyzed it in the same way as all the data in our program. For all wavelength bins we find the star to be non-variable to very good precision, as expected. The results are shown in Figure 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

View figure set (23 panels)

Tarball of figure set images

Because the error bars in each point are of similar order as the measured flux changes, the random noise can influence the slope that is ultimately measured. To determine the probability density function for the true slope a_t given the measured slope a_m and error bars, we use Bayes' theorem, that states

$$\begin{equation} p(a_t | a_m) = \frac{ p(a_m | a_t) p (a_t) }{ \int _0^\infty p(a_m | a_t) p(a_t) \mathrm{d} a_t }, \end{equation} \tag{ 1 }$$

where p(a_t|a_m) is the posterior probability distribution, p(a_m|a_t) the likelihood function, and p(a_t) the prior probability density function for variability slopes.

We perform a Monte Carlo simulation to determine the likelihood function for the slope for each source and each wavelength range. For a range of slopes, we add random Gaussian noise with the appropriate standard deviation and measure the new slopes. Repeating this 50,000 times, we measure how often the measured slope occurs, where we accept the slope to be equal if it is within a 1σ interval of the measured slope. For an smaller acceptance interval, the likelihood function would be slightly narrower. The resulting likelihood function, normalized to the maximum likelihood, is shown in Figures 7–10 as the solid black line.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Choosing an appropriate prior distribution is not obvious because this is what we ultimately want to measure in this study. Using a flat prior (which would equate the posterior probability distribution with the likelihood function) is inappropriate in this case because a_t does not have an upper bound. We know from earlier studies that strong variability is rare. We therefore adopt an exponentially declining function as prior $p(a_t) \propto e^{-ba_t}$ where we use two different values of b to explore a more optimistic and a more pessimistic case and study the resulting differences. The priors are plotted as dashed blue or red lines in Figures 7–10, and the resulting posterior probability distribution as solid blue or red lines (both normalized to their maximum value). Because none and very small variability slopes are favored for this prior, which may not necessarily be true, the derived amplitudes may be slightly underestimated. However, the likelihood function corresponds to the posterior probability for a flat prior with a large upper bound, which can be used for comparison. In any case, the choice of the prior has only minimal influence on the solid detections of variability, where the likelihood function is very small for slopes near zero. However, the measured large slopes in wavelength regions with large error bars, most notably for the water band of 2M1624, are overwhelmed by the prior. Our short duration observations are therefore not able to draw firm conclusions about variability in deep absorption bands, where the S/N is intrinsically lower.

From the resulting posterior probability distributions, we compute the 68% and 95% credible intervals (highest density) for the slopes. These Bayesian confidence intervals denote the range in which the parameter lies with 68% or 95% probability. They do no correspond to frequentist confidence intervals because our prior is not uninformative. These intervals, together with the maximum likelihood and maximum probability for the slope, are tabulated in Table 4 for the different wavelength bands for all sources. All slopes are given in units of % relative flux change per hour in order to be able to compare between sources with different visit durations. We classify a detection of variability as significant if the lower boundary of the 95% confidence interval level for the strong prior (L95, P2 in Table 3) is >0.3% hr⁻¹, and as tentative if the lower boundary of the 68% confidence interval for the strong prior (L68, P2 in Table 3 is >0.3% hr⁻¹). An exception to this rule is elaborated in the next section for a case where the strong prior overwhelms even a strong likelihood function. The slopes and confidence intervals for sources with confident or tentative variability are listed in Table 3. In the following discussion of the variability, we generally refer to the slope a as column Max P2 in Table 3 unless otherwise noted.

As an alternative, we also calculate χ² values. We derive the χ² probability p(null) of the null hypothesis that the data can be explained by a flat line. We also provide the reduced χ² value for the best-fit linear model. On average, the reduced χ² is 0.82, indicating that our error bars may be slightly too large. This is mainly due to the additional 0.1% error that we add to all points for the uncertainty of the ramp correction, and 0.2% for points at the beginning and end of the time series. For all confident variables from our Bayesian analysis, p(null) < 0.01 in at least one wavelength band, most often significantly lower. For the bright objects and the broader filters, where the ramp correction error dominates the total error, the p(null) value is overestimated. On the other hand, for a few non-variable cases with larger scatter than explained by the error bars, p(null) may be underestimated. We do not consider the results from the χ² analysis for our further discussion, but provide them in Tables 3 and 4 for comparison.

Figures 7–10 show the observations as well as the likelihood, prior distribution, and posterior distribution for the true slope for each target for selected wavelength regions with a significant signal for confident and tentative variability. Figures 4.2–4.23 available in the online journal include all selected wavelength bands for all sources and Table 4 all confidence intervals including upper limits for non-variables.

Out of seven objects with spectral types between L5 and L7, we find two which show significant variability. The most obvious is 2M1750−00 (L5.5; Figures 7 and 4.21), one of the brightest objects in our sample. It shows a clear downward trend with a ∼ 0.75% hr⁻¹ in J band that is significant at >95% level. The same result is found when looking at narrower wavelength regions between 1.12 and 1.32 μm, indicating that the amplitude is quite uniform across this wavelength range. In the H band, a smaller trend (<0.5% hr⁻¹) remains, but it is not statistically significant or even tentative.

The second object showing clear variability is 2M0624−45 (L5, Figures 7 and 4.8). For this object, a linear model is a very poor fit. Both in J band H band (as well as the narrower flux peaks), the light curve appears to reach a minimum at ≈15 minutes. The variability is of similar shape and strength at all wavelengths within error bars. We therefore discuss only the white light curve integrated over the full spectrum from 1.12 to 1.66 μm in order to minimize the error bars. We divide the curve into two linear slopes and apply our Bayesian analysis only to the second, longer slope. We find a >95% significant slope with a > 1% hr⁻¹.

The other five objects do not show any indication of variability above the uncertainties (see Figures 4.5, 4.6, 4.9, 4.15, and 4.20). The 68% upper limits are below 0.5%–1% hr⁻¹ for most wavelengths for 2M0421−63, 2M0539−00, 2M0801+46, and 2M1515+48 and 2M1711+22 (Table 4). The small trends in the longest wavelength bins for 2M0421−63 were found to originate from small shifts in the wavelength at the spectrum edge (see Section 2.4) and are therefore not real.

We observed five objects with spectral types between L7.5 and L9.5 and find two objects with significant, and one object with tentative variability. An interesting significant detection is made for the binary 2M0310+16 (L9+L9; Figures 8 and 4.4), where we find a very strong slope a ∼ 2% hr⁻¹ in the 1.26–1.32 μm J-band flux peak with >95% confidence. A separate extraction of the two spectra was not possible. Since it is likely that the variability stems from only one component, the true amplitude for this object is likely even larger. This is reminiscent of the variable component in the very nearby binary WISEJ104915.57−531906.1AB (Gillon et al. 2013). For 2M0310, we could not derive the broadband J and H variability because of a large number of bad pixels shortward of 1.26 μm and longward of 1.65 μm.

Another significant variable is 2M0825+21 (L7.5; Figures 8 and 4.11). For this object we find trends with a ∼ 1% hr⁻¹ that are similar over the whole wavelength range. Finally, 2M1219+31 (L8; Figures 8 and 4.15), is a curious case. At 1.12–1.20 μm, we find a quick ∼2% drop and then a slower ∼2% rise. Formally, the Bayesian analysis on only the rising slope shows that it is significant at >95% confidence, with a slope of 3%–6% hr⁻¹ depending on the prior. However, for this faint source the analysis is dependent on only very few points. In the 1.22–1.32 μm region, we also find an upward slope at >68% confidence, but no downward slope. At longer wavelengths, there is no evidence of variability. We therefore classify this source only as a tentative, and not as a confident variable.

The remaining two objects are not found to be variable: for both 2M0908+50 and 2M1632+19 (Figures 4.12 and 4.19) the 68% upper limits are mostly <0.6% hr⁻¹. The tentative trend in the H band for 2M1632+19 is not real for the same reason as for 2M0421−63.

Our observed sample contains only three early-T dwarfs, and none show variability at >95% significance. We find tentative trends (>68% significance) for one of the three objects, but only in one particular narrow wavelength region: for 2M1039+32 (T1; Figures 9 and 4.14) we find a tentative trend at 1.55–1.61 μm, which corresponds to the H-band peak. For 2M1324+63 (T2.5pec; Figure 4.16), the upper limit on the variability at all wavelengths is <0.6% (68% confidence). For the faint 2M0909+65 (T1.5; Figure 4.13), there is a trend in the 1.61–1.65 μm CH₄ dip in the H band, but it does not fall into our classification of tentative variability.

With two confident and three tentative variables out of seven objects, the mid-T dwarfs show the greatest variety in variability (Figure 10). The most prominent variable is 2M0559−14 (T4.5; Figures 10 and 4.7). It shows significant variability in the H band with a > 1% hr⁻¹. The strongest change is seen in the CH₄ dip in the 1.61–1.65 μm region, but the H-band peak shows tentative variability as well. There is also a small tentative variability trend in the J-band peak, but in opposite direction. If real, this would indicate that the light curves are out of phase.

The second significant variable is 2M0817−61 (T6; Figures 10 and 4.10), where we find confident variability in the peak region of the J band at 1.22–1.30 μm with a slope of ∼0.7% hr⁻¹. There is also tentative variability in the H-band peak in the same direction, while wavelengths beyond 1.60 μm are missing for this source.

An interesting source is 2M1624+00 (T6; Figures 10 and 4.18), which shows a strong maximum likelihood slope of a ∼ 14% in the deep water absorption band at 1.35–1.44 μm. Because the individual error bars are about 2%, the likelihood function is broad. For our chosen priors such large slopes are very unlikely, therefore the posterior probability is strongly influenced by the prior for this case, making the variability formally tentative for the weaker prior and non-existent for the stronger case. However, because there is no prior data on the occurrence rate of variability in deep water absorption bands, the prior, which is based on broadband observations, may be too stringent here. We know that a slope of several percent is not impossible: for 2M2228−43, Buenzli et al. (2012) found a maximum slope of a = 7.6% hr⁻¹ at the same wavelengths. With a flat prior the variability would be confident. Because errors are large, the flux is very low and the prior important, we adopt this variable as tentative. We do not find variability at other wavelengths, but several regions are missing due to overlap with a background star.

We find two other tentative variables. One is 2M2339+13 (T5; Figures 10 and 4.23), which shows a tentative downward trend for the whole integrated spectrum, most of which is due to variability in the J-band peak. On the other hand, 2M1750+17 (T3.5; Figures 10 and 4.22) shows a tentative trend in the H band with a > 1%. Because of the sources' faintness and the very short duration of the visit, the prior influences the resulting probability distribution.

The remaining two sources, 2M0000+25 (T4.5; Figure 4.2) and 2M0243−24 (T6; Figure 4.3) are not variable above the uncertainty level. For 2M0423−24, the small trends in the H band are artificial and due to wavelength shifts.

Four of our targets have previously been monitored for variability in the same wavelength range and here we compare our findings to earlier results.

A target of particular interest is 2M0559−14. This bright mid-T dwarf is overluminous across the near- and mid-IR (e.g., Dahn et al. 2002; Dupuy & Liu 2012) and has therefore been suspected to be an unresolved flux-equal binary. Liu et al. (2008) speculate that the object may be marginally resolved in their HST/WFPC2 image. Our direct image with an FWHM of 015 shows no indication of binarity, but the PSF is undersampled and we cannot rule out a very tight binary. Alternatively, it may represent the peak of the cloud clearing, however this does not explain the mid-IR overluminosity. 2M0559−14 was monitored for variability by Clarke et al. (2008) in J band. They find it to be non-variable at a level of 0.5%. While we see a small trend in J band, the observed flux change is<0.5%and tentative only when integrated over the narrow flux peak. The confident variability is seen only in H band and most strongly in the absorption band at 1.61–1.65 μm, which has not been monitored before. It seems unlikely that the variability stems solely from remnant clouds, because a cloud impacting an absorption band at low pressure level would also affect the deeper layers probed by the J band. Alternatively, circulation patterns may result in temperature perturbations in some layers as suggested for 2M2228−43 by Buenzli et al. (2012).

For 2M1711+22, Khandrika et al. (2013) report very strong (∼20%) but only marginally significant variability in a short 0.7 hr sequence. Our observations cannot confirm this detection. We find no evidence of variability with an upper limit of ≈2%. They also monitored 2M0825+21 in J band, but for only 0.5 hr and did not find variability, whereas we find a change of about 0.8% in 40 minutes. This is likely below the precision level of their very short measurement.

Nakajima et al. (2000) attempted to measure spectral variability for 2M1624+00 over an 80 minute time-span. They find tentative evidence for variability in water absorption features between 1.5 and 1.6 μm. We find potential strong water variability between 1.35–1.44 μm for this source, a wavelength region that is not accessible at sufficient precision from the ground. For the same source, Koen et al. (2004) did not find signs of variability in J band, consistent with our findings, although we are missing part of the J band due to overlap with a background star.

We derive the occurrence rate of brown dwarf variability for mid-L to mid-T spectral types from our sample by determining the binomial confidence interval considering k detections out of a sample of n = 22 brown dwarfs. We found a total of six confident (>95%) and five tentative (>68%) detections of variability. Because our observations only span fractions of a rotation period and are still precision-limited, the derived variability fraction must be considered to be a lower limit. We derive the 1σ confidence interval following the Appendix in Burgasser et al. (2003b). With 6 confident detections for a sample of 22, the minimum variability fraction is $f_{{\rm min}} = 27^{+11}_{-7}\%$. Including the tentative detections, if they are all real, the minimum fraction raises to f_min = 50% ± 10%. The true number of the minimum fraction of variables from our sample is therefore likely to lie between about one third to half of brown dwarfs.

Furthermore, it must be noted that the fraction of brown dwarfs with intrinsic heterogeneities is expected to be larger than the measured variability fraction: we can only detect rotationally asymmetric components and the variability signal is reduced for inclinations toward pole-on. Additionally, the survey may miss variables with very long periods or amplitudes below the photometric precision. We therefore expect that heterogeneity is a very important property of the condensate clouds and one that should be accounted for in ultracool atmosphere models.

In order to derive limits on the possible periods and amplitudes for our sample, we determine how the measured slope in an observation of 37 minutes (the average length of our observations) duration relates to the actual amplitude and period of the object. We simulate sinusoidal light curves for a grid of amplitudes and periods at random phases and measure the distribution of different slopes. We find a strong peak in the frequency of slopes close to the maximum possible slope when assuming an equatorial view (Figure 11, top panel). The maximum slope for any period/amplitude combination for a sine curve is given analytically by the derivative at zero phase: a_max = πA/P, where A is the peak-to-peak amplitude and P the period. For any measured slope a, we can therefore estimate the amplitude to be at minimum A_min = aP/π when assuming a period P. For example, for 2M0559−14, with a = 1.2% hr⁻¹ in H band, we expect a peak-to-peak amplitude A > 1% for P > 2.7 hr. At 1.61–1.65 μm it is even A > 1.5%. Similarly, it is likely that A ≳ 1% for 2M0825+21 and 2M0624−45 in J and H bands if P > 3 hr and for 2M1750−00 in J band if P > 4 hr. However, these calculations assume variability with a sine curve shape, which does not necessarily have to be the case.

Zoom In Zoom Out Reset image size

If a brown dwarf in not seen from an equatorial view, the observed amplitude will be reduced. We therefore also account for the inclination of the brown dwarf that reduces the amplitude of the light curve by cos i, assuming random orientations in the sky. In this case, we find that for any combination of period and true amplitude A^⋆, the probability of measuring a particular slope is distributed uniformly between 0 and the maximum possible slope (Figure 11, bottom). Although unlikely, an inclination near pole-on will yield a non-detection of variability even if the brown dwarf atmosphere is highly patchy.

Figure 12 shows the maximum possible slope as a function of both period and amplitude. Overplotted are contours that correspond to typical measured slopes. We find slopes of at most ∼2% hr⁻¹, with the exception of 2M1624+00 in the deep water band, where the slope could be >10% hr⁻¹ if we do not adopt a prior that strongly favors smaller slopes. Period/amplitude combinations to the upper left of a contour line are not possible for that particular slope under the assumption that curves are sinusoidal. Combinations to the lower right, i.e., shorter period and larger amplitude, while not impossible, become increasingly unlikely further away from the contour. Extreme combinations with low periods and high amplitudes can generally be excluded for robust non-detections because they would exhibit significant curvature, even at maximum and minimum phases. This range is indicated by the dashed line.

Zoom In Zoom Out Reset image size

We calculate the number of detections one would expect in our survey for two cases: all brown dwarfs have heterogeneities that result in variability with A = 1.5% when seen equatorially, or brown dwarf amplitudes are uniformly distributed between 0%–3% when seen equatorially. We assume a uniform period distribution between 2–10 hr and calculate the distribution of slopes for a large sample of 10,000 objects with random orientation in the sky. We then calculate for 22 targets the number of objects expected with a slope larger than a given number (Figure 13). For both cases, the distribution is similar and ≈7–8 targets would have slopes >0.5%, while about 1 object would have a slope >1.5% (this is consistent with our detections). Our results therefore are consistent with all brown dwarfs hosting low-level variability, but are also consistent with a uniform distribution of amplitudes from zero to a few percent.

Zoom In Zoom Out Reset image size

Previous surveys have discovered a handful of confident high-amplitude near-IR variable brown dwarfs, nearly all of these with spectral types between T0 and T2.5 at the L/T transition. However, when including less confident detections of variability, the combination of previous ground-based surveys does not indicate an excess of variables in the L/T transition with respect to before and after (Khandrika et al. 2013). One notable confident variable is the T6.5 dwarf 2M2228−43 (Clarke et al. 2008; Buenzli et al. 2012) which lies beyond what is generally regarded as the L/T transition region. Ground-based surveys are typically limited by their precision to detect broadband variables with amplitudes A ≳ 1%–2%. With broadband filters, which cover a broad pressure range in the atmosphere, a strong signal originating from a thin atmospheric layer may be diluted. This is true in particular for 2M2228−43, where narrow wavelength regions through the J and H bands have significantly different phases and amplitudes. On the other hand, the two early-T dwarfs that have been studied with HST/WFC3 spectroscopy (Apai et al. 2013), SIMP0136+09 and 2M2139+02, both show variability with a surprisingly weak wavelength dependence. There, the amplitudes and phases in narrow wavelength regions do not differ much (except in the deep water absorption bands where the variability is lower), therefore broadband observations are not at a disadvantage for these objects (not surprisingly, as these were both identified from precision-limited broadband surveys).

In this survey, we reach sub-percent precision and are able to study narrow wavelength regions for the first time over a statistically significant sample of brown dwarfs. It is therefore instructive to determine whether the trend of finding the majority of confidently variable brown dwarfs in the L/T transition holds, or whether we can confirm a relatively uniform distribution of variability throughout the L and T spectral types. As shown in Figure 6, our survey supports the latter. Indeed, we do not find a convincing variable brown dwarf in our (arguably small) sample of early-T dwarfs, although one of the three objects shows a tentative trend in the 1.5–1.7 μm wavelength region. On the other hand, we find convincing variables in all other spectral bins: two confident mid-Ls (out of seven), two confident and one tentative late-Ls (out of five), and two confident and three tentative mid-Ts (out of seven). Counting only the confident variables we derive the (minimum) frequency of variability in each spectral bin: $f_{\mathrm{min,mid\hbox{-}L}} = 29^{+20}_{-10}\%$, $f_{\mathrm{min,late\hbox{-}L}} = 40^{+21}_{-16}\%$, f_{min, early-T} = 0⁺³⁷%, and $f_{\mathrm{min,mid{\rm -}T}} = 29^{+20}_{-10}\%$, with error bars giving the 1σ confidence interval. Figure 14 shows the number and fraction of variables per spectral bin. We do not find significant differences between the different spectral bins, but with three to seven objects per spectral bin the uncertainties are still large. Furthermore, our result is consistent with the statement that high-amplitude variability (several percent) is rare both inside and outside the L/T transition.

Zoom In Zoom Out Reset image size

A second notable result is the fact that we find diverse variability as a function of wavelength. For the L dwarfs 2M0825+21 and 2M0624−45 we find similar variability levels across the whole spectrum, this is similar to the known variable early-T dwarfs. For the other two confidently variable L dwarfs, 2M0310+16 and 2M1750−00, the variability originates in the J band while it is not evident in the H band. For the T dwarfs, we find variability sometimes in the J band peak, the H band peak, or methane or water absorption bands. It is not clear where these differences arise from. It may simply be that the amplitude of variability is lower in one or the other wavelength region, perhaps because the perturbations stem from a specific pressure level in the atmosphere. Or the variability is shifted in phase, and it is therefore not detected in our short observations for some wavelengths. In order to understand the color dependence of the variability, which can pose strong constraints on atmospheric models, longer observations that cover most of a rotation period are required.

Because of the low amplitudes and sometimes limited wavelength regions where variability occurs, broadband photometric surveys are likely to miss it for many brown dwarfs and may therefore underestimate the occurrence rate of variability. High-precision spectroscopy from space-based instrument such as HST/WFC3 makes it possible to use spectral mapping to study clouds and weather phenomena in three dimensions for many brown dwarfs across the L- and T spectral types. In the near future, with next generation adaptive optics systems and with the James Webb Space Telescope, the same technique can be applied to extrasolar planets (e.g., Kostov & Apai 2013).

What is the origin of this low-level variability that seems to be a common occurrence for many brown dwarfs? Within the L/T transition, models have shown that patchy cloud cover can explain the optical and near-IR spectra of brown dwarfs (e.g., Marley et al. 2010). Asymmetric distribution of patches then results in variability at the several percent level. However, we also clearly find variability for mid- to late-L type dwarfs where models generally predict a thick cloud cover, and for mid- to late-T dwarfs where silicate clouds are expected to be below the visible photosphere. For the L dwarfs, it seems reasonable to think that the variability could stem from heterogeneous cloud cover of thicker and thinner cloud patches, which Apai et al. (2013) also suggested for the early-T dwarfs 2M2139+02 and SIMP0136+09. For these two objects, the variability amplitudes are either constant as a function of wavelength or slightly larger in the J band than the H band, in agreement to what we find for our variable L dwarfs. While for early-L dwarfs magnetic spots may also contribute to variability (e.g., Lane et al. 2007), there is no indication that this may also be true for later type brown dwarfs. For example, Gelino et al. (2002) argue that these photospheres are largely neutral, and they do not find any correlation between I-band variability and Hα emission. For the later T dwarfs, cloud types other than silicates may appear that could also be patchy (Morley et al. 2012). For the known variable T6 dwarf 2M2228−43 models with sulfide and chromium clouds, as well as with clouds of a species with similar optical properties to iron, provided a good match to the average spectrum (Buenzli et al. 2012), while cloud free models did not. However, the light-curve phase shifts as a function of wavelength show that not only clouds may be responsible for variability: other possibilities include temperature fluctuations from circulation patterns or gas opacity perturbations. The diversity in variability that we find as a function of wavelength may point toward several mechanism where different ones dominate for different objects. Considering that even the solar system gas giants appear variable at some wavelengths (e.g., Gelino & Marley 2000; Karkoschka 2011), it may be reasonable to expect low level heterogeneities in all atmospheres where condensates form.