# TESS Observations of the Luhman 16 AB Brown Dwarf System: Rotational Periods, Lightcurve Evolution, and Zonal Circulation

Dániel Apai, Domenico Nardiello, Luigi R. Bedin

► To cite this version:

Dániel Apai, Domenico Nardiello, Luigi R. Bedin. TESS Observations of the Luhman 16 AB Brown Dwarf System: Rotational Periods, Lightcurve Evolution, and Zonal Circulation. The Astrophysical Journal, 2021, 906, 10.3847/1538-4357/abc97. insu-03690847

HAL Id: insu-03690847

https://insu.hal.science/insu-03690847v1

Submitted on 2 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# TESS Observations of the Luhman 16 AB Brown Dwarf System: Rotational Periods, Lightcurve Evolution, and Zonal Circulation*

Dániel Apai¹,², Domenico Nardiello³,⁴, and Luigi R. Bedin⁴

¹Steward Observatory, The University of Arizona, 933 N. Cherry Avenue, Tucson, AZ 85721, USA; apai@arizona.edu

²Lunar and Planetary Laboratory, The University of Arizona, 1629 E. University Blvd., Tucson, AZ 85721, USA

³Aix Marseille Univ., CNRS, CNES, LAM, Marseille, France

⁴INAF—Osservatorio Astronomico di Padova, Vicolo dell'Osservatorio 5, I-35122, Padova, Italy

Received 2020 April 29; revised 2020 October 26; accepted 2020 October 27; published 2021 January 7

## Abstract

Brown dwarfs were recently found to display rotational modulations, commonly attributed to cloud cover of varying thickness, possibly modulated by planetary-scale waves. However, the long-term, continuous, high-precision monitoring data to test this hypothesis for more objects is lacking. By applying our novel photometric approach to Transiting Exoplanet Survey Satellite data, we extract a high-precision lightcurve of the closest brown dwarfs, which form the binary system Luhman 16 AB. Our observations, which cover about 100 rotations of Luhman 16 B, display continuous lightcurve evolution. The periodogram analysis shows that the rotational period of the component that dominates the lightcurve is 5.28 hr. We also find evidence for periods of 2.5, 6.94, and 90.8 hr. We show that the 2.5 and 5.28 hr periods emerge from Luhman 16 B and that they consist of multiple, slightly shifted peaks, revealing the presence of high-speed jets and zonal circulation in this object. We find that the lightcurve evolution is well fit by the planetary-scale waves model, further supporting this interpretation. We argue that the 6.94 hr peak is likely the rotation period of Luhman 16 A. By comparing the rotational periods to observed $v \sin(i)$ measurements, we show that the two brown dwarfs are viewed at angles close to their equatorial planes. We also describe a long-period ($P \sim 91$ hr) evolution in the lightcurve, which we propose emerges from the vortex-dominated polar regions. Our study paves the way toward direct comparisons of the predictions of global circulation models to observations via periodogram analysis.

*Unified Astronomy Thesaurus concepts: Exoplanet atmospheres (487); Atmospheric circulation (112); T dwarfs (1679); Atmospheric variability (2119); Stellar rotation (1629); Photometry (1234)*

## 1. Introduction

Planetary and brown dwarf atmospheres are fundamentally impacted by the presence of condensate clouds (e.g., Burrows et al. 2001; Robinson & Marley 2014; Zhang 2020), and several studies have showed that the properties of clouds themselves are impacted by atmospheric dynamics (e.g., Showman & Kaspi 2012, 2013; Showman et al. 2020). In weakly irradiated atmospheres (such as directly imaged exoplanets and brown dwarfs) the interplay of rotation and heat transfer will set the nature of large-scale atmospheric circulation (e.g., Showman & Kaspi 2012; Zhang & Showman 2014). The presence or absence of clouds will impact the local pressure-temperature profiles which, in turn, can lead to the formation/evaporative dissipation of clouds (Tan & Showman 2017), triggering self-sustaining cycles of cloud formation and, possibly, driving planetary-scale waves. Atmospheric layers can also be heated by the breaking of atmospheric waves emerging from high-pressure turbulent layers, generating planetary-scale waves (Tan & Showman 2019).

The presence of planetary-scale waves has recently been invoked to explain data from long-term monitoring of L/T spectral type transition brown dwarfs (Apai et al. 2017). In three brown dwarfs (two of which are of planetary mass), Apai et al. (2017) observed continuous lightcurve evolution over timescales of 1–1000 rotational periods. They showed that bright or dark features alone (i.e., one to few Great Red Spot-like features) cannot explain the symmetric, wavelike lightcurves (as spots tend to introduce asymmetric variations).

Instead, Apai et al. (2017) demonstrated that planetary-scale waves can provide a simple explanation for the observed lightcurve evolution—a simple model of a sum of two to three sine waves (with unconstrained amplitude, phase, and period) successfully fitted all lightcurve segments. This study raised two questions: (A) Do brown dwarfs other than the three studied exhibit lightcurve evolution that is consistent with the planetary-scale model? (B) What processes drive the planetary-scale waves? Answering these two questions requires long-term, high-precision, photometric monitoring of a larger sample of brown dwarfs, but the end of the Spitzer Space Telescope mission led to the loss of continuous infrared monitoring capabilities. The study presented here demonstrates a novel approach to obtain the required data.

Although the studies described focus on brown dwarfs, it is paramount for the reader to understand that brown dwarfs—including planetary-mass objects—also provide a gateway to understanding directly imaged exoplanets. In fact, L/T spectral type transition brown dwarfs are excellent analogs to the directly imaged exoplanets in the HR 8799 system (Marois et al. 2010).

This study focuses on the L/T spectral transition brown dwarf binary WISEJ104915.57-531906.1 or, in the following, Luhman 16. Even though it is the closest known brown dwarf system to the Sun ($d = 1.9955 \pm 0.0004$ pc, Bedin et al. 2017), Luhman 16 was only discovered in 2013 (Luhman 2013). The binary consists of a $34.2_{-1.1}^{+1.3} M_{\text{Jup}}$ L-type primary and a $27.9_{-1.1}^{+1.0} M_{\text{Jup}}$ T-type secondary brown dwarf component (Ammons & Garcia 2019). The primary’s spectral type is L7.5 and the secondary’s spectral type is T0.5

<footnote>* Based on observations collected with the NASA TESS.</footnote>
(Burgasser et al. 2013; Kniazev et al. 2013). Right after its discovery Luhman 16 was identified as one of the most highly variable brown dwarfs (Biller et al. 2013; Gillon et al. 2013) and follow-up spectroscopic observations showed spectrophotometric modulations (Burgasser et al. 2014) that resembled those of the few other L/T transition brown dwarfs studied with similar methods (Apai et al. 2013). High-precision spectrophotometric observations in combination with high spatial resolution revealed that both A and B components are variable (Biller et al. 2013), although the integrated modulations are dominated by the modulations in the T0.5-type B component of the system and not by the L7.5-type A component (e.g., Burgasser et al. 2014; Crossfield et al. 2014; Buenzli et al. 2015b). Detailed modeling of the spectrophotometric modulations showed that spatially correlated cloud thickness–temperature modulations (invoked for other L/T transition brown dwarfs, see Radigan et al. 2012; Apai et al. 2013) can also explain the behavior of Luhman 16 (Buenzli et al. 2015a, 2015b). Luhman 16 is one of the relatively small number of brown dwarfs bright enough to allow time-resolved medium- and high-resolution spectroscopy that can probe different atomic and molecular tracers to further constrain the properties of the heterogeneous cloud cover. Buenzli et al. (2015a) reported that FeH absorption is very strongly correlated with the continuum brightness, possibly arguing for a lack of deep opacity holes that would allow sight lines into the hotter interior where FeH would be present in gas phase. Crossfield et al. (2014) used a Doppler-imaging CO inversion technique to present a possible map of the cloud pattern, but one which is not sensitive to potential longitudinal banded structures. Kellogg et al. (2017) found that the K I line is inversely correlated with the continuum brightness variations, suggesting that the cloud thickness variations significantly change the pressure–temperature profile in the upper atmosphere.

In a new study Millar-Blanchaer et al. (2020) published a novel H-band, spatially resolved polarimetric study of the Luhman 16 system. This study measured a 0.031% ± 0.004% polarization level for Luhman 16 A and 0.010% ± 0.004% polarization level for Luhman 16 B. Detailed modeling showed that while the lower-level polarization of Luhman 16 B can be explained by either oblateness or by banded cloud structures, Luhman 16 A's stronger polarization is only consistent with the presence of banded clouds and not with oblateness alone. The Millar-Blanchaer et al. (2020) study assumed a rotation period of 5 hr or 8 hr for Luhman 16 A and 5 hr for Luhman 16 B.

Beyond the spectrophotometric information that can shed light onto the cloud properties, the most perplexing property of Luhman 16 B is the prominent evolution of its lightcurve. Modeling of multiepoch lightcurve snapshots (typically covering 1 rotational period) identified a recurring lightcurve feature and was able to explain the lightcurves with few large spots (Karalidi et al. 2016). However, the study of Apai et al. (2017) showed that while single-rotation lightcurves can often be explained with spots or waves, lightcurve evolution observed in brown dwarfs over multiple rotations could not be explained by spots and required planetary-scale waves. As a nearby, bright, and highly variable L/T transition brown dwarf, Luhman 16 provides the ideal target to test the planetary-scale model.

In this paper we present a unique, long-term monitoring observation of Luhman 16 with NASA's Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014), covering about 100

rotational periods. In the following we review the observations and data reduction, we present our lightcurve analysis, and review and discuss the results. Finally, we summarize the conclusions of our study.

## 2. Observations and Data Reduction

### 2.1. TESS Photometry

Luhman 16 AB were observed by the TESS (Ricker et al. 2014) during Sector 10, between 2019 March 26 and 2019 April 22. The Luhman 16 system was only observed on the full-frame images (FFIs) mode, with a cadence of 30 minutes.

Because of the low resolution, the Luhman 16 AB system is not resolved on the TESS FFIs; as a consequence we extracted the combined lightcurve for the system.

For the extraction of the lightcurve, we used our point-spread function (PSF)-based approach described in detail in Nardiello et al. (2019) for the case of TESS time series, and already adopted for ground- and space-based time series (see Nardiello et al. 2015, 2016a, 2016b and Libralato et al. 2016a, 2016b; Benatti et al. 2019).

Briefly, by using FFIs, accurate PSF models, and a high-resolution input catalog (in this work the one from Gaia DR2, Gaia Collaboration et al. 2018), we extracted high-precision photometry of the target source after subtracting adjacent field stars using two different photometric methods: aperture photometry (with 1-, 2-, 3-, 4-TESS pixel radii) and PSF-fitting photometry.

Because Luhman 16 AB is a close-by high-proper-motion source, we calculated the position of the system at the epoch of observations (~2019.26) by using the astrometric positions calculated by Bedin et al. (2017). We adopted ($α_{2019.26}$, $δ_{2019.26}$) = (162.303282427, -53.317573814), that corresponds to a time-averaged (over the observing sequence) detector position with the raw pixel coordinates of (1767.4, 7.4) on the CCD 1 of Camera 2 (see Figure 1).

After the extraction, we corrected the systematic trends that affect the lightcurve by using cotrending basis vectors (CBVs) extracted from a set of 2000 stars across the whole CCD. We applied the CBVs as in Nardiello et al. (2019).

The corrected, normalized lightcurves extracted with all the photometric methods are shown in Figure 2. We computed the P2P rms, as calculated in Nardiello et al. (2019), for all the photometric methods to select the best lightcurve<fnref n="5" />. We found that the lower P2P rms are for the lightcurves obtained with 1 px aperture photometry (~9.8 mmag) and PSF-fitting photometry (~10.6 mmag), while the other photometric apertures show higher P2P rms (~14.1 mmag, ~28.7 mmag, ~42.8 mmag, for 2 px, 3 px, and 4 px aperture photometries, respectively). In this work we use the 1 px aperture photometry lightcurve (in magenta in Figure 2).

In much of the current study we explore the periods present in the TESS lightcurve of Luhman 16 AB. Given the relatively large pixels of TESS, false periods appearing due to spacecraft positional jitter may be a concern. In other words, a brighter star within or close to the aperture may contribute to the measured light and, if the telescope's position is oscillating, it may introduce an apparent periodicity. Given TESS's extremely stable pointing, such a contamination is very unlikely. Nevertheless, in order to conservatively test this possibility, we

<footnote n="5">5 P2P rms is not affected by the variability of the source</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (Left) Finding chart of Luhman 16 AB on TESS image (tess2019107105933). Axes are aligned north up, east to the left, and the field of view (FOV) is about 66" × 60". (Right) Same FOV from the Hubble Space Telescope (HST) (from Bedin et al. 2017). Green lines indicate the pixel grid on the TESS images. The blue lines are the astrometric solution for the barycenter of the Luhman 16 AB system. The red circles highlight the adopted position for the barycenter at the epoch 2019.26 during which TESS observations were acquired, the smallest of which has a radius of 220 mas (indicative of the smallest separation between the two components). The green squares are the sources in the Gaia DR2 catalog.</figcaption></figure>

calculated the Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) the same way as we do for the lightcurve analysis in subsequent sections (Section 4). The resulting periodogram is shown in Figure 3, along with the window function (top panel). We include the Lomb–Scargle periodograms of the x- and y-positional differences, as well as the difference from the aperture center. The Lomb–Scargle periodogram does display peaks, demonstrating a very slight oscillation with multiple periods. However, comparison of these periodograms to the peaks in the periodogram of the astrophysical source shows that none of the peaks identified in Luhman 16 AB’s lightcurve correspond to peaks in the periodogram of the aperture position. The lack of a match between the peaks in the Luhman 16 AB and the positional modulation periodograms demonstrates that the Luhman 16 AB lightcurve periodicity does not emerge from spacecraft oscillations.

The lightcurve will be released on the MAST archive as HLSP under the project PATHOS<fnref n="6" /> (see Nardiello et al. 2019 for details).

## 2.2. HST Lightcurves

Another potential, though unlikely, source of contamination in our TESS Luhman 16 AB lightcurve could be introduced by background stars within or close to the TESS aperture. Any star that would introduce such a contamination would need to be relatively bright (comparable to Luhman 16 AB) and would need to display strong variability (at a multipercent level). Although TESS does not allow for easy disambiguation of multiple sources in the aperture, we were able to utilize already published HST Wide Field Camera 3 ultraviolet/visible channel (WFC3/UVIS) image series, obtained in the broadband F814W filter. The image mosaic covers an area of about 160" × 160". These images have been collected between 2014 and 2016 in 12 epochs, with the goal to monitor the astrometric

evolution of the binary system. They have been analyzed and published in Bedin et al. (2017). These data sets included high-precision photometric measurements for 808 stars in the vicinity of Luhman 16 AB.

Our goal was to test the hypothesis that one or more relatively bright star exists within the vicinity of the Luhman 16 AB aperture. We thus searched the multiepoch photometry for bright, high-amplitude variable stars close to Luhman 16 AB. In Figure 4 we plot the measured brightness (F814W band) of stars against the standard deviation of their lightcurves (as a measure of their potential variability). For reference, we also plot Luhman 16 B’s observed brightness and semiamplitude from (Bedin et al. 2017). In the figure we highlighted the six stars that are of comparable brightness to Luhman 16 AB ($m_{814W} \sim 15.2$ mag, Bedin et al. 2017). For each of these stars we also include labels that show the projected distances of these stars from Luhman 16 AB. Based on this information, we conclude that there are no stars of comparable brightness to Luhman 16 AB within the TESS aperture; and that the closest, relatively bright star is 60" away, i.e., about $3 \times$ the radius of the TESS aperture used in our study. These facts make it unlikely that any of the stars contaminate the lightcurve. Furthermore, as explained in the following sections, the Luhman 16 AB TESS lightcurves display large-amplitude variations (up to ~10%); this level of variation is about $4 \times$ higher than the most variable stars seen in our data. We note a minor caveat to this analysis: the central pixels of the few brightest stars may be nonlinear in the HST images, thus, their brightness may be slightly underestimated (Bedin et al. 2017). Nevertheless, the rest of the pixels are well sampled and the low standard deviations further bolster confidence in our findings. As one final step and additional test, we also phase-folded the lightcurves of the 10 stars with the 90.8 hr period, a particular period identified later in this analysis. This analysis also turned up negative: none of the stars showed even low-amplitude variability with the 90.8 hr period.

<footnote n="6">6 https://archive.stsci.edu/hlsp/pathos, doi:10.17909/t9-es7m-vw14</footnote>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Lightcurves of the Luhman 16 AB system extracted from TESS data obtained with different photometric methods: PSF-fitting (black) and aperture photometry with radius 1 px (magenta), 2 px (blue), 3 px (green), and 4 px (orange).</figcaption></figure>

In summary, we conclude that due to the lack of bright stars in the vicinity of Luhman 16 AB and due to the lack of highly variable stars anywhere in the field, it is extremely unlikely that the periods observed in the Luhman 16 AB TESS lightcurves are significantly contaminated by a background star.

**3. Results: Modulations and Lightcurve Evolution**

During the entire duration of the data set, the lightcurve of Luhman 16 remained variable. The greatest relative brightness difference observed in the data set was about ~13% (between the maximum and minimum of Segment 1), but the typical peak-to-trough brightness changes over the timescale of a single rotation were typically about 4%. The lightcurve shows long-term modulations, as well as short-term (~1 rotational period) modulations. The Luhman 16 lightcurve displays a complex, yet familiar behavior (see Figure 8), reminiscent of those reported in short-duration data sets for Luhman 16 AB

(Gillon et al. 2013; Buenzli et al. 2015a, 2015b; Karalidi et al. 2016) and other L/T brown dwarfs (Apai et al. 2017). The lightcurve shows dozens of peaks over the nearly 50 day long baseline of the observations, but no obvious, simple periodicity; instead, continuous evolution is apparent in the lightcurve.

Based on the inspection of the lightcurve we conclude the following: (1) Luhman 16 remains variable during the entire data set, covering over 540 hr, suggesting that the atmosphere itself remained heterogeneous. (2) The lightcurve is not strictly periodic but continues to evolve, confirming previous reports of rotational modulations and lightcurve evolution (Gillon et al. 2013; Karalidi et al. 2016) in Luhman 16. (3) Both short-term (comparable to the rotational period of ~5.4 hr, see Section 6) and long-term (60–100 hr) variations are apparent.

In the following we will explore the lightcurve evolution via a Lomb–Scargle periodogram analysis (Section 4.4), followed by an exploratory lightcurve modeling with the planetary-scale wave model (Section 5).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Lomb–Scargle periodograms of the window function (top panel), and the X- and Y-positions of the TESS aperture, and for its distance from the median (X, Y) coordinate. The LS periodograms of the spacecraft position have multiple peaks in the explored period range. However, the short-period (*P* < 20 hr) range is clean of major peaks as is the vicinity of the long-period (*P* ~ 90.8 hr) peak identified in the lightcurve. Therefore, we conclude that slight spacecraft positional jitter does not explain or is expected to impact the periods identified as physical in the lightcurve.</figcaption></figure>

**4. Periodogram Analysis**

We now turn to a generalized Lomb–Scargle periodogram analysis (see Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009, and, for a review of applications and limitations, VanderPlas 2018) to explore the underlying periodic components within the Luhman 16 intensity variations. This section presents the periodogram analysis, along with the study of potential window function artifacts, and the results of the analysis.

**4.1. Lomb–Scargle Periodogram and Window Function Analysis**

Lomb–Scargle periodogram analysis can be performed on data sampled on an irregular grid, but such data sets will inevitably introduce a nonuniform power distribution in the periodogram. Although the TESS data set provides an excellent lightcurve with close to perfect sampling, the impact of the window function on the periodogram must always be carefully considered. To aid the interpretation of the Lomb–Scargle periodogram we calculated the window function’s generalized Lomb–Scargle periodogram (using the PyAstronomy package, Czesla et al. 2019, with Zechmeister–Kürster normalization, Zechmeister & Kürster 2009). The window function was

calculated over a finely and regularly sampled array where the function’s value was set to unity at the point when measurements were available and zero elsewhere. This window function was then analyzed through the Lomb–Scargle method to allow the identification of the sampling’s imprint on the periodogram.

In Figure 5 we show the resulting Lomb–Scargle periodograms for the Luhman 16 data set (in blue) and for the window function (red). The two upper panels show the 0–140 hr period range, while the two lower panels show the same for the 0–20 hr period range, allowing a finer inspection. As a comparison, the bottom panel shows the periodogram of a noise-added sine-wave function sampled at the same time as the TESS observations. The period of the sine wave is *P* = 5.28 hr and its amplitude is 0.05. The noise added to this sine wave resembles the photon noise estimated for our data (1σ = 0.013).

Based on Figure 5 we conclude that the window function’s periodogram is nearly featureless for periods between 1–20 hr, due to TESS’s efficient sampling of the lightcurve. The window function, however, displays prominent broad and prominent peaks at 86 and 120 hr, and additional multiple narrower and weaker peaks corresponding to periods between 20 and 64 hr.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Brightness and standard deviation of the stars in the Luhman 16 AB field from HST photometry (12 epochs, Bedin et al. 2017). All stars that are comparable in brightness to Luhman 16 AB and show some variability (blue symbols) are far in projection from Luhman 16 AB and very unlikely to contaminate the TESS aperture. In addition, none of the stars shows variability approaching that detected in the TESS lightcurve of Luhman 16 AB. We conclude that contamination from a field star is extremely unlikely.</figcaption></figure>

## 4.2. Interpretation of the Peaks

Our periodogram analysis leads to the following conclusions:

(1) A confident detection of periodicity corresponding to 5.28 hr, which likely corresponds to the rotational period of the Luhman 16 component that dominates the rotational modulations.

(2) The peak corresponding to the *P* = 5.28 hr period modulations is slightly broader in the data than that of the sine wave (bottom panel). We explore this finding in more detail in Section 4.3.

(3) Two additional, but much weaker, peaks are observed at *P* = 2.5 hr and *P* = 6.94 hr. These do not correspond to peaks that would emerge from the window function or from the *P* = 5.28 hr sinusoidal modulation. Peaks in the power spectra corresponding to half-periods have been reported in Spitzer lightcurves of brown dwarfs (Apai et al. 2017), and this may explain the *P* = 2.5 hr period peak. The *P* = 6.94 hr period peak, we speculate, may correspond to the rotational period of the Luhman 16 component that has a lower amplitude in the TESS band. We note that through injection of sine waves with these periods we verified that these two peaks are *present independently* of each other, i.e., neither the 2.5 hr or the 6.94 hr peak is an alias of the other peak, or emerges from a combination of the other peak and the 5.28 hr peak.

(4) We also identify two peaks in the power spectrum that correspond to long-term modulations: a peak at period 90.81 hr and a peak at period 126.6 hr. Both of these peaks lie very close to peaks prominent in the window function's periodogram, which raises the possibility that these are not physical features. To help their interpretation, we calculated the false alarm probabilities (FAPs) for these peaks. We found that the *P* = 5.28 hr and *P* = 90.81 hr peaks have very low FAPs (log(FAP) of −46 and −44, respectively). In contrast, the

*P* = 126.6 hr peak's log(FAP) is much higher (−13). Given these considerations, we conclude that the peak at period 90.8 hr is likely physical, while the peak at period 126.6 hr is likely not physical. Although the window function complicates the interpretation of such long-period peaks, supporting our interpretation is the fact that the TESS lightcurve shows very clear long-term evolution (see, for example, Figures 2 and 8). Although a period close to 90.8 hr is very likely real, given the contamination from the window function, we consider its exact period not well determined.

## 4.3. Fine Structures of the 5.28 hr and 2.50 hr Peaks

We will now explore the nature of the *P* = 5.28 hr peak in the periodogram. In Figure 6 we plot the periodogram data (in blue, top panel) for the *P* = 4.0–6.00 hr period range), along with simple model fits (panels below top panel, in green). The periodogram data (blue) displays multiple peaks: cursory inspection suggests possible peaks around 4.75, 5.00, 5.22, 5.28, and 5.45 hr. Of course, the periodogram peaks are not trivial to interpret: a single sinusoidal modulation will introduce multiple peaks in the periodogram (see, for example, the bottom panel of Figure 5 or Figure 6(b)). We attempt to interpret the observed complex periodogram signature by models of increasing complexity consisting of periodogram signatures of sine waves. Our models are a single sine wave (panel (b)), three sine waves (panel (c)), four sine waves (panel (d)), and six sine waves (panel (e)). The model fits were carried out by minimizing the sum of the squares of the data–model residuals through a Levenberg–Marquardt optimization. Our models are linear combinations of the periodograms of a sine wave, shifted in period space and scaled in intensity
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Lomb–Scargle periodograms of the TESS data (blue), its window function (red), and a noise-added sine wave (green) with the same temporal sampling as the data. The top four panels show the same data and window periodograms but for different period ranges, to allow closer inspection. The 1–20 hr period range is not affected by the window function. A clear physical periodicity is identified at *P* = 5.28 hr. Additional possible periods may be present at *P* = 2.5 hr and *P* = 6.94 hr.</figcaption></figure>

(corresponding to shifts in period and intensity). In order to ensure that the Levenberg–Marquardt optimization is not trapped in a local minima we repeated the optimization several thousand times, starting from randomly determined initial guesses, and accepted the best-fit model. The similarities of the solutions for different model complexities demonstrate that our modeling procedure is very robust.

We find that a single sine wave fails to match the complex periodogram signature, i.e., there is clear evidence for the presence of multiple periods. We find that models of increasing complexity succeed in reproducing well the multipeaked structure of the periodogram (Figures 6(c)–(e)). Our most

complex model has six periods (panel (e)), which fits the periodogram of the data very well, although not perfectly. Our model fits provide compelling evidence for the presence of multiple, similar periods in Luhman 16 B.

Our fits show the following periods (in order of decreasing power): 5.29, 5.22, 4.98, 5.46, and 5.02 hr. We stopped at a model with six sine waves and did not increase the model’s complexity further. We opted to do this because, while more complex sine-wave based fits are possible, the high-quality data on the fine structure of the periodogram calls for future models more deeply rooted in physical models, such as a comparison to predictions by global circulation models (Showman et al. 2019),
Luhman 16B Rotational Period Periodogram and Models

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Lomb–Scargle periodograms of the TESS data (blue), focusing on the rotational period of Luhman 16 B, and models consisting of sine waves modeled in the period space.</figcaption></figure>

the predictions of which may be converted to periodograms and compared to the observations.

As a first-order metric to describe the width of the substructure, we provide an estimate for the relative range of periods contained in it. Expressed as $\frac{(P_{\text{max}} - P_{\text{min}})}{P_{\text{min}}}$ and adopting $P_{\text{min}} = 4.75$ hr, $P_{\text{max}} = 5.77$ hr values, we find that the relative period range is ~22%.

The planetary-scale wave model by Apai et al. (2017) would predict that, in addition to $k = 1$ wavenumber waves, $k = 2$ waves would be present. It is thus instructional to explore whether the $P = 2.5$ hr peak also shows a multipeaked substructure as the $P = 5.28$ hr peak does. In Figure 7 we magnified the Lomb–Scargle periodogram, focusing on the period range 2.0–3.0 hr. Indeed, the periodogram displays a substructure here, too, with multiple peaks in the approximate period range 2.4–2.7 hr. We note that as there is much less power in this peak than in the 5.28 hr peak, the signal-to-noise ratio is more limited, and the substructure is less well defined. Therefore, we will not attempt the detailed modeling of the substructure here. Instead, to allow a first-order comparison to the 5.28 hr peak structure, we calculate the approximate

relative period range in which peaks appear to be above the background noise. Expressed as $\frac{(P_{\text{max}} - P_{\text{min}})}{P_{\text{min}}}$ and adopting $P_{\text{min}} = 2.4$ hr, $P_{\text{max}} = 2.7$ hr values, we find that the relative period range is ~12%.

In summary, we conclude that both the $P = 5.28$ hr and the $P = 2.5$ hr peaks display fine structures, with multiple peaks, and that the peaks in the 2.5 hr group are likely to be the $k = 2$ (half-period) waves corresponding to the $k = 1$ waves in the 5.28 hr group. The high-quality $P = 5.28$ hr peak group reveals the presence of multiple—at least six—distinct, slightly different rotational periods.

4.4. Long-period Modulation ($P = 90.8$ hr)

We note here, however, that there is a remote possibility that the long-period modulation ($P = 90.8$ hr) is introduced by a background star that falls within the same aperture. Although we have done extensive testing based on the available data and found no indication for such a scenario, this cannot be excluded with the present data at hand. In the future TESS campaigns, however, by revisiting the high-proper-motion Luhman 16,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Lomb–Scargle periodograms of the TESS data on Luhman 16 B, focusing on the peaks around period *P* = 2.5 hr.</figcaption></figure>

contamination by a background star (with negligible proper motion) will be possible to identify.

**5. Lightcurve Modeling**

Our Lomb–Scargle analysis showed that in Luhman 16 B's lightcurve multiple peaks are present with slightly different periods, in two groups: one centered on the rotational period of the object ~5.28 hr and one centered at the half-period ~2.52 hr. These findings are fully consistent with the lightcurve evolution model proposed by Apai et al. (2017). That model is based on the comprehensive Spitzer/IRAC monitoring obtained in Apai et al. (2017) and explains the lightcurve evolution with planetary-scale waves with similar amplitudes, but slightly different periods (likely due to differential rotation). Apai et al. (2017) also reported the presence of *k* = 1 and *k* = 2 waves (where *k* = 1 corresponds to a single peak during a rotational period). We will now explore how well this model can actually fit the TESS lightcurves of Luhman 16 AB presented here.

We applied the same model described in Apai et al. (2017) to the TESS lightcurve with a few minor modifications. These changes are required to run the algorithm efficiently on the TESS lightcurve, but are not expected to lead to any significant difference between the model published in Apai et al. (2017) and our implementation. For the details of the model and extensive test results we refer to Apai et al. (2017) and its supplementary online material, and here only summarize the fundamental setup of the model and the minor changes unique to our implementation for the TESS data.

The planetary-scale wave model describes the evolving lightcurve evolution as a sum of three sine waves. The periods, amplitudes, and phases of the three sine waves are mostly unconstrained. The three sine waves are fit through a two-stage fitting process: first a nonheuristic algorithm (Genetic Algorithm) is used to identify the likely global optimum, which is then passed on as initial guess to an efficient multidimensional gradient optimizer to find the exact best-fit solution.

We made two modifications to the original procedure: First, we constrained the range of the periods fitted to 2–6 hr, the full range of periods identified by the Lomb–Scargle analysis. This ensured that the parameter space searched is consistent with the periodogram results, but without prescribing or favoring actual periods within that range. Second, in our model we did not attempt to fit the very long-period (*P* = 90.8 hr) component of

the lightcurve, but focused on the evolution of the short-period component. Therefore, we fitted a copy of the lightcurves with the long-term evolution removed by subtracting a 10-sample (5 hr window) boxcar-smoothed version of itself (fits using a 7.5 hr boxcar window led to identical results).

Figure 8 shows the entire lightcurve split into two segments (Segment 1 and Segment 2, separated by a data downlink gap). Within this figure we marked five segments to which we applied the planetary-scale wave models (Segments 1abc and 2ab). The beginning and end times of the five segments were chosen somewhat arbitrarily, with the intent to identify overall representative segments of different lengths, to enable testing the planetary-scale wave model. Multiple other segments have been fitted, although not shown in the paper, with qualitatively identical results.

Figures 9–13 show our model fits to five lightcurve segments. In the upper panels of these figures, filled circles show the segment-normalized TESS photometric measurements (with colors reflecting the brightness of the object, i.e., the ordinate of the data points). The panel also includes the key parameters of the fits (periods and amplitudes of three sine waves). The lower panels show the residuals (data–model) in black filled circles, and the model itself (offset to 0.0 mean ordinate value).

We find that the lightcurve evolution over all five segments fitted is very well reproduced by the model. The standard deviation of the residuals remains very close to 0.01 for the fits, including the longest segment (100 hr, Segment 1c). It is instructive to see how this simple model succeeds in matching even segments with apparently complex lightcurve evolution. For example, in Segment 1b (Figure 10) in the course of 40 hr the lightcurve evolves from a very low-amplitude, seemingly aperiodic section into a high-amplitude, flat-topped section with deep, short minima. Or, for example, in Segment 1 (Figure 11), where in 100 hr the lightcurve goes from low-amplitude sine-like variations to high-amplitude sine-like variations, before returning. In these cases, as in all other segments, the model reproduces the lightcurves well.

A key question in the model fits is whether the fitted periods are, in fact, consistent with the Lomb–Scargle periodogram analysis’ finding that the periods should form two groups (2.5 and 5.28 hr), specifically: how often do the fitted periods fall within the period ranges identified in the Lomb–Scargle periodogram, and how often are periods are inconsistent with those ranges?
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Overview of the TESS lightcurve of Luhman 16 AB. Within the two segments four regions are marked (Segments 1abc, 2ab), which are magnified and modeled in Figures 9–13.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Upper panel: Segment 1a of the Luhman 16 lightcurve (minus the long-term trend) and the planetary-scale waves model (blue curve). Bottom panel: residuals (data–model) and the model. The planetary-scale waves model reproduces the evolution of the lightcurve well with fitted period values that are fully consistent with those found in the Lomb–Scargle analysis.</figcaption></figure>

We find, in fact, that all (15) periods fitted by the model are in the ranges identified by the Lomb–Scargle periodogram analysis. The most deviant value is period 3 (P₃ = 5.94) in the shortest-fitted lightcurve, Segment 1a (Figure 9), but even this is relatively close to the group 4.75–5.75 hr period group.

In short, our lightcurve fits demonstrate that: (a) even long and complex sections of the lightcurve are very well reproduced by the planetary-scale wave model; and, (b) the best-fit periods concentrate into the same two groups identified by the Lomb–Scargle analysis.

# 6. Discussion

## 6.1. Luhman 16 A or Luhman 16 B? Or both?

An important challenge in interpreting the TESS data presented here is the fact that the TESS aperture includes both A and B components of the Luhman 16 system. Given the very

small projected separation (<0.′′3, Bedin et al. 2017) of the two components, it is not possible to identify the relative contributions of the two components to the TESS lightcurve based solely on the TESS data.

Nevertheless, a comparison of our results to published spatially resolved observations allows us to argue that the lightcurve presented here is (1) very likely dominated by a single component, and (2) that that single component is Luhman 16 B.

First, our study of the Lomb–Scargle periodogram identifies a single very prominent periodicity (P = 5.28 hr) that is responsible for almost all of the short-term modulations. A much lower-amplitude periodicity is also seen at P = 6.9 hr, which we tentatively attribute to the other component in the system. These data show that only a single component dominates the lightcurve: the alternative, i.e., both components contributing significantly, would require their rotational periods
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Upper panel: Segment 1b of the Luhman 16 lightcurve (minus the long-term trend) and the planetary-scale waves model (blue curve). Bottom panel: residuals (data–model) and the model. The planetary-scale waves model reproduces the evolution of the lightcurve well with fitted period values that are fully consistent with those found in the Lomb–Scargle analysis.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Upper panel: Segment 1c (100 hr) of the Luhman 16 lightcurve (minus the long-term trend) and the planetary-scale waves model (blue curve). Bottom panel: residuals (data–model) and the model. Covering about 20 rotational periods, our longest segment modeled is still reproduced well by the planetary-scale waves model, with fitted period values fully consistent with those found in the Lomb–Scargle analysis.</figcaption></figure>

and amplitudes to match very closely (periods to within a few percent and amplitudes to a factor of ~2). Such a coincidence in rotational period and amplitudes is extremely unlikely, leaving us with the conclusion that a single component (with $P = 5.28$ hr rotation period) dominates.

Second, the comparison to the literature measurements strongly suggests that this component is Luhman 16 B. Although both components are variable (e.g., Biller et al. 2013; Buenzli et al. 2015b), the few spatially resolved studies all found B to display a 2–3 times higher amplitude in the near-infrared than A does. For example, spatially resolved HST WFC3 G102 and G141 spectrophotometry showed that, for wavelengths longer than 0.8 µm, component B has at least twice as great an amplitude (4.6%) than component A (2.2%) (Buenzli et al. 2015a, 2015b). Although continuous, spatially resolved, complete-rotation lightcurves are not available for wavelengths shorter than 0.8 µm, spatially resolved multiepoch HST photometry of the system (Bedin et al. 2017) showed that, while components A and B have similar amplitude

(~0.05 mag) in the F606W filter ($\lambda_c$ = 606 nm), in the redder F814W filter component B's amplitude (~0.12 mag) is more than twice of that of component A (~0.05 mag).

In comparison, TESS's single band covers wavelengths between approximately 600–1000 nm, with its central wavelength of $\lambda_c$ = 786 nm. Therefore, although no complete, spatially resolved lightcurves are available for components A and B in the TESS band, the data from multiple studies consistently and very strongly suggests that component B has at least a 2 times higher band-integrated amplitude in the TESS band than component A. In addition, the Luhman 16 system is a flux-reversal binary: the overall brighter primary Luhman 16 A is actually slightly fainter in the Y and J-bands than component B, further decreasing its relative contribution to the TESS band.

Finally, there is also evidence from the periodicity of the objects. While A's period remains poorly constrained in the 5–8 hr range (e.g., Millar-Blanchaer et al. 2020), B's amplitude has been established to be very close to 5.2 hr: Gillon et al. (2013)
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Upper panel: Segment 2a of the Luhman 16 lightcurve (after subtraction of the long-term trend) and the planetary-scale waves model (blue curve). Bottom panel: residuals (data–model) and the model.</figcaption></figure>

found a rotational period of $4.87 \pm 0.01$ hr. Shortly after, Burgasser et al. (2014) fitted the lightcurve with a period of $5.05 \pm 0.1$ hr, a longer period than estimated by Gillon et al. (2013). Mancini et al. (2015) presented additional data supporting a longer rotational period ($5.1 \pm 0.1$ hr). These rotational period measurements are completely consistent with the rotational period of 5.28 hr we found in the Lomb–Scargle periodogram analysis (see Section 4.4), and very close to the value ($P \sim 5.4$ hr) found in the fits of the lightcurve segments via the planetary-scale wave model (Section 5).

Therefore, we conclude that evidence strongly suggests that the lightcurve is dominated by modulations of a single component, and the facts that component has a $P = 5.28$ hr rotational period and that B has larger amplitude in the TESS band strongly suggest that the component dominating the lightcurve is, in fact, component B.

We identify this period with the rotational period of 5.28 hr derived from the Lomb–Scargle periodogram analysis (see Section 4.4), which we interpret as the rotational period of Luhman 16 B. Therefore, the $P = 5.28$ hr rotational period is fully consistent with recent literature values as well as the results of our periodogram analysis and the planetary-scale wave model fits.

## 6.2. The Rotational Periods and Inclinations of Luhman 16 A and B

The combination of our finding of one dominant peak ($P = 5.28$ hr) in the periodogram with literature-based identification of B's rotational period of 5.2 hr represents a step forward in understanding the Luhman 16 AB system. Based on our results B's rotational period is $P = 5.28$ hr (and it likely exhibits a small range of differential rotation). Furthermore, we tentatively attribute the 6.9 hr peak in our periodogram to the rotational period of component A.

By combining the rotational periods for components A and B with measured radial rotational velocities for each component, we can also explore possible spin axis orientations in the system, i.e., constrain the inclinations of the two brown dwarfs. Assuming brown dwarf radii between 0.90 and $1.10 R_{Jup}$ (e.g., Burrows et al. 2001), the equatorial rotational velocities expected from the rotational periods would be $v_{A,eq} = 16.3–19.9$ km s$^{-1}$

and $v_{B,eq} = 21.2–25.9$ km s$^{-1}$. In contrast, Crossfield et al. (2014) measured projected rotational velocities of $17.6 \pm 0.1$ km s$^{-1}$ for Luhman 16 A and $26.1 \pm 0.2$ km s$^{-1}$ for Luhman 16 B. Thus, the observed $v \sin(i)$ measurements are within $1\sigma$ of the range of the equatorial rotational velocities calculated. This comparison is informative and shows that the measured rotational periods are fully consistent with the observed $v \sin(i)$ values.

Furthermore, our analysis suggests a close to equatorial viewing angle for both brown dwarfs. Based on our simple model, Luhman 16 A is viewed from an angles within 28° from its equatorial plane (i.e., $i > 62°$); and Luhman 16 B is viewed almost exactly equatorially (within a few degrees, i.e., $i \simeq 90°$). Our analysis does not provide evidence for misaligned spin axes in the system; on the contrary, the two brown dwarfs may have very similarly aligned axes. As the Luhman 16 system's orbital motion has been closely monitored, it is possible to compare the spin axes of the brown dwarfs to that of their orbital spin axis: Bedin et al. (2017) found the orbital inclination to be $79.21 \pm 0.45$, which suggests that the three rotational and orbital spin axes in the system may be well aligned.

We note that the power of the $P = 6.94$ hr peak, which we attribute to Luhman 16 A, is significantly ($> 10\times$) lower than that of the $P = 5.28$ hr peak, corresponding to Luhman 16 B's rotation. While Luhman 16 B's amplitude is lower, as discussed above, the relatively large difference between the two peaks in the TESS data is surprising. We speculate that this may suggest that the TESS observations occurred during a period when Luhman 16 A's amplitude was lower than usual. Future revisits of the Luhman 16 system may help test this hypothesis and shed more light on the periodogram signal from Luhman 16 A.

Finally, we note that a close-to-equatorial viewing angle is also consistent with the fact that the Luhman 16 B is among the highest-amplitude variable brown dwarfs known. High amplitudes are much more likely in systems seen in near-equatorial orientations (e.g., Metchev et al. 2015; Vos et al. 2017), as found here.

## 6.3. A Novel Look Into Brown Dwarf Atmospheres

A stunning feature of the Luhman 16 lightcurve presented here is its complexity, i.e., its nonperiodic, evolving nature.
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Upper panel: Segment 2b of the Luhman 16 lightcurve (after subtraction of the long-term trend) and the planetary-scale waves model (blue curve). Bottom panel: residuals (data–model) and the model.</figcaption></figure>

Although complex lightcurve evolution has been reported for a number of brown dwarfs (e.g., Radigan et al. 2012; Apai et al. 2013; Gillon et al. 2013; Metchev et al. 2015; Karalidi et al. 2016), most of these lightcurves have been too short and/or had too low signal-to-noise ratios to allow testing hypothetical models for lightcurve evolution. In contrast, the lightcurve shown here provides a relatively high signal-to-noise ratio and long temporal baseline. The only somewhat comparable lightcurve data is from the Spitzer Space Telescope Exploration Science program Extrasolar Storms (see Yang et al. 2016; Apai et al. 2017). The Extrasolar Storms data provided higher cadence and higher signal-to-noise ratio infrared (3.6 and 4.5 µm) photometry on six brown dwarfs, with overlapping HST spectrophotometric snapshots. While the Extrasolar Storms data covered up to ~1000 rotation periods for each target in eight visits, no visits covered more than four rotation periods. In contrast, the TESS data presented here provides a unique, quasi-continuous coverage of a single target over ~100 rotation periods, thereby providing data on the previous unexplored 4–40 rotation timescale. Furthermore, the current data set is the only high-precision, sustained monitoring carried out in the visible wavelength regime for a brown dwarf with continuously evolving lightcurves. Thus, the data set provides a powerful comparison to the shorter-coverage infrared Spitzer monitoring data that are available for a larger set of objects.

## 6.4. Lightcurve Evolution Model

Based on the qualitative observation of the lightcurve (Figure 8 and Section 3) and the quantitative model fits described in Section 5, we conclude that the visual lightcurve of Luhman 16 is very similar in nature to the lightcurves observed in the infrared for this and other brown dwarfs (Buenzli et al. 2015a, 2015b; Karalidi et al. 2016; Apai et al. 2017). Specifically, we identify four properties that are shared between the visual lightcurve of this object and the infrared lightcurves of other objects: (1) The lightcurves remain variable over long periods (years). (2) The lightcurve shape evolves, yet it displays characteristic period, which is likely the rotational period of the object (as found in Apai et al. 2017). (3) In spite of the rapid evolution of the lightcurve, the amplitudes over rotational timescales remain similar and characteristic to the

object. (4) The lightcurves tend to be symmetric in the sense of a similar number of positive-negative features, in contrast to, for example, a situation in which a single positive feature appears periodically on an otherwise flat lightcurve, which would indicate a single bright spot in the atmosphere. Therefore, the qualitative similarity of the Luhman 16 lightcurve to the data obtained for other L/T brown dwarfs in the Extrasolar Storms program motivates the test of the planetary-scale waves model developed for those objects (Apai et al. 2017). In Figures 9–13 we applied the planetary-scale wave model to lightcurve segments of various lengths, ranging from 25 hr (Segment 1a in Figure 9), 40 hr (Segment 1b, Figure 10), 50 hr (Segments 2a and 2b, Figures 12 and 13), and 100 hr (Segment 1c, Figure 11). Each of these figures includes the key fit parameters (wave periods and amplitudes), as well as the data–model residuals.

We found that not only does the simple planetary-scale wave model provide an excellent fit to all lightcurve segments fitted, but that best-fit periods in the lightcurve segments correspond closely the periods identified in the periodogram analysis. Thus, the combination of these findings strongly supports the applicability of the planetary-scale waves to Luhman 16 B, i.e., strongly argues for the presence of multiple zones in the atmosphere with slightly different rotational periods—in other words, it argues for zonal circulation.

## 6.5. Zonal Circulation in Luhman 16 B

As revealed by the periodogram and by the lightcurve segment fits, a multitude of similar, but slightly different effective periods exist in the modulations of Luhman 16 B. A very similar periodogram has been compiled from K2 (Kepler extended mission, Howell et al. 2014) photometry of Neptune (Simon et al. 2016), which was found to be similar to the periodogram derived from Spitzer lightcurves of two L/T transition brown dwarfs (Apai et al. 2017). However, the data presented here provides a much higher-quality periodogram in which much finer details are visible. The periodogram presented here reveals not just a simple split in the peak corresponding to the rotational period, but a multitude of peaks (Figure 6). Furthermore, for the first time, it also provides similar insights into the weaker peak that corresponds to *k* = 2
(half-period) waves (Figure 7). This “peak,” too, is resolved to multiple similar peaks.

The multipeaked nature of the periodogram peaks demonstrates that Luhman 16 B has zonal circulation, i.e., it sports latitudinal zones of eastward and westward jets. These high-speed jets rotating prograde and retrograde will modulate the rotational periods of any structures embedded in these zones: structures embedded in prograde (eastward jets) will display shorter rotational periods and structure embedded in westward jets will display longer rotational periods. Simulations (e.g., Showman & Kaspi 2012; Showman et al. 2019; Imamura et al. 2020) of giant planets and brown dwarfs, as well as observations of Jupiter (e.g., Limaye 1986; Porco et al. 2003) suggest that the equatorial jet is likely to be the fastest and also the widest jet, although its direction may oscillate between eastward and westward (Showman et al. 2019).

With these insights in mind, we interpret Figure 6 in the following way. The strongest period detected ($P=5.29$ hr) likely corresponds to the period of the equatorial jet, i.e., it is slightly different from the true rotational period (interior) of Luhman 16 B. Peaks that are shorter than the 5.25 hr are likely emerging from structures embedded in prograde (eastward) jets, while structures embedded in retrograde (westward) jets are responsible for the peaks at periods longer than ~5.29 hr. Whether the equatorial jet itself is prograde or retrograde cannot, at this point, be determined.

We note that the number of peaks fitted to the periodogram do not correspond directly to the number of jets presents in the atmosphere, as any number of jets with similar speeds will appear as a single peak. Nevertheless, the number of peaks may place a lower limit on the number of zones, under three reasonable assumptions: (1) wind speeds in a given zone are constants, (2) wind speeds are approximately symmetric to the equator, and (3) there are no transition regions between the zones. Such a simple model would suggest that Luhman 16 B harbors at least six distinct zones in velocity space—or, given equatorial symmetry—a total of 12 zones (or six band/belt pairs).

Although the model described above is useful as a first approximation, we note that the latitudinal wind velocity distribution in real objects is unlikely to be this simple. In fact, Cassini observations of wind speeds in Jupiter (Limaye 1986; Porco et al. 2003), as well as circulation models of brown dwarfs, show no sudden breaks in the latitudinal wind speed distribution, but rather smooth transitions. Thus, models based on a number of belts with distinct velocities (or effective periods) have only limited validity. We further explore these considerations with a simple model in Section 6.6 and note that future studies may attempt more direct comparisons between circulation models and periodograms.

Nevertheless, our data does allow the exploration of the relative wind velocities presented in Luhman 16 B. Although the exact period range in which peaks are present is not well determined, we find evidence for at least a 22% relative period range (Figure 6). Given the rotational period and the typical diameter of a brown dwarf, the maximum wind speed difference corresponds to about 4.5 km s⁻¹. Interestingly, this value is about an order of magnitude greater than that measured in Jupiter or predicted by baseline brown dwarf models (Showman et al. 2019). We note that, perhaps, this discrepancy may be resolved if the structures seen are not just passive tracers embedded in the jet systems, but truly perturbations that propagate within the jet

systems. For example, radiative-convective feedbacks have been proposed to introduce periodic cloud formation/dissolution (Tan & Showman 2017). If such perturbations do exist, then the measured velocities would correspond to the sum of the jet velocities and the propagation speed of the perturbations, thus significantly narrowing the gap between predicted wind speeds and observed periods. It is also useful to place the observed wind speed range of 4.5 km s⁻¹ in the context of the recent study by Allers et al. (2020), who derived the wind speed in a brown dwarf relative to its interior by comparing its near-infrared rotational period to its radio-based rotational period. This study found an eastward wind speed of 650 ± 310 m s⁻¹ for a rapidly rotating ($P=1.77 \pm 0.04$ hr) T6 spectral type brown dwarf. Considering the uncertainties, this is about a factor of three lower than the wind speed range we derived for the slower-rotating and hotter Luhman 16 B. This difference, however, should not be surprising as the radio-near-infrared period comparison is likely measuring the relative speed between the dominant jet and the brown dwarf’s interior, while our measurements capture the velocity differences between the fastest eastward and the fastest westward jets, and thus are expected to be higher.

## 6.6. Comparison to Jupiter Model

To place the Luhman 16 B periodogram results in context, we will now explore how Jupiter—if observed similarly to Luhman 16 B—may compare to them. Jupiter is, of course, not a perfect analog to Luhman 16 B, as they differ in effective temperature, surface gravity, and—very likely—in composition. Yet, their general atmospheric circulation may share similarities and a comparison may provide some insights. There is no identical data set available, and we attempt here to build an admittedly simple model for Jupiter, aiming at a qualitative, exploratory comparison.

Using wind speed measurements (as a function of latitude) for Jupiter, we will calculate a histogram of the effective period across the atmosphere of the planet. The fundamental, simplifying assumption we make here is that every visible atmospheric element contributes to the signal that emerges from Jupiter equally, adding its own periodicity. In other words, the larger fraction of the atmosphere has a given wind speed, the stronger the effective period that correspond to that wind speed value will appear in the histogram.

Figure 14(a) shows the wind speed values as a function of latitude, as measured by the Cassini mission (from Porco et al. 2003). The highest wind speeds (~130 and ~70 m s⁻¹) were observed at low latitudes, close to the equator, while more moderate wind speeds (typically between −10 and +50 m s⁻¹) were observed at mid- to high latitudes. We note that a very similar wind speed pattern was also measured in Voyager 1 and Voyager 2 image pair differences (Limaye 1986) about two decades earlier, therefore the wind speed pattern shown is characteristic to Jupiter’s current atmospheric circulation state. As our goal is to translate the wind speed pattern into a histogram of effective periods, we map the wind speeds ($v_{wind}$) observed onto a simulated disk of Jupiter (approximating its shape with a sphere), as shown in Figure 14(b). We then calculate an effective period map for our synthetic Jupiter, by adopting an equatorial rotational period of $P_{eq} = 9$ hr 56 min, a radius of 1 Jupiter radius, and calculating, for each atmospheric element, the rotational period that emerges from the combination of the local wind speed and the bulk rotation
<figure><img src="image_16.png" /><figcaption>Figure 14. A simple model for Jupiter to support comparison to the Luhman 16 B periodogram. (a) Wind speeds measured as a function of latitude in Jupiter by the Cassini team (Porco et al. 2003). (b) The observed wind speeds mapped to a spherical disk. (c) Effective periods calculated from the sum of the Jupiter’s bulk rotation and the latitudinally varying wind speeds.</figcaption></figure>

of the planet:

$$P_{\text{eff}}(\lambda) = \frac{2\pi R_{\text{Jup}} * \cos(\lambda)}{v_{\text{rot}} + v_{\text{wind}}(\lambda)},$$

where λ denotes the latitude, R$_{Jup}$ is 1 Jupiter radius, v$_{rot}$ = 2πR$_{Jup}$/P$_{eq}$ is the rotational velocity corresponding to the bulk rotation of Jupiter. The resulting map is shown in panel (c). In this synthetic map, rotational periods range from 9 hr 49 min to 9 hr 59 min.

We will now contrast the predictions of this very simple Jupiter model to the observations of Luhman 16 B. In the top panel of Figure 15 we plot the histogram of the effective periods of the photospheric elements in our Jupiter model. In the lower panel of the same figure we plot, for comparison, the Lomb–Scargle periodogram results for the primary rotational peak (P = 5.28 hr) of Luhman 16 B (from Section 4.3). We find that, in spite of the simplicity of our Jupiter model and the limited amount of information on Luhman 16 B, there are some interesting qualitative similarities between the results. First, in both cases there is a range of effective periods present in the vicinity of the objects’ rotational periods. Second, the distributions are asymmetric: there are longer tails toward shorter periods than toward longer periods, suggesting that larger parts of the atmosphere are characterized by high-speed prograde jets than by retrograde ones, and that the prograde jets have higher velocity. Third, in both distributions there appears to be evidence for distinct peaks, i.e., there are likely some typical, characteristic jet speeds rather than a smooth

distribution of wind speeds without peaks. The comparison also reveals two differences. First, the relative range of periods predicted for Jupiter is (as mentioned in Section 6.5) about an order of magnitude lower than that seen in Luhman 16 B. Second, the periodogram of Luhman 16 B shows more distinct period peaks than those seen in our synthetic Jupiter model. The latter point may simply be a result of our simple model, where we assumed that each atmospheric element contributes equally to the periodogram. In reality, Jupiter’s thermal emission is very unevenly distributed and opacity holes— infrared bright spots—are strongly correlated with the atmospheric circulation (e.g., Showman & Dowling 2000). Due to this correlation between wind speed and local thermal infrared brightness, it is very likely that less simplistic models of Jupiter’s periodogram will display more prominent peaks.

The former difference, however, between the relative range of periods in Jupiter and the Luhman 16 B lightcurve is very unlikely to be the result of our simplistic model, as this aspect of the prediction is the direct consequence of the measured wind speeds. The difference between Luhman 16 B and Jupiter, in terms of their relative period range, does suggest that an additional, different process is occurring in the brown dwarf than in Jupiter. The propagation of planetary-scale waves (Apai et al. 2017), and/or waves caused by radiative-convection feedbacks and cloud formation (Tan & Showman 2017) may be candidates to explain the different relative period ranges. Our model assumes a perfectly equatorial viewing angle, which is —as explained in Section 6.2—likely very close to Luhman 16 B’s orientation. Nevertheless, we will briefly address here how
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Top panel: histogram of the effective periods in our Jupiter model (Figure 14(c)). Bottom panel: the multipeaked structure in the periodogram of the Luhman 16 B rotational peak (*P* = 5.28 hr). In spite of the important differences between Jupiter and Luhman 16 B and the simplicity of our Jupiter model, there is an overall qualitative similarity between the distributions of effective periods, supporting the conclusion that Luhman 16 B has zonal circulation, as Jupiter does.</figcaption></figure>

the predicted period distribution may change for more inclined (10°–30°) viewing angles. Unlike for a Doppler-shift-based measurement, the change in the viewing angle for rotational period-based measurements will not result directly in the compression of the deduced velocity/period differences; this is because the measurements probe the rotational periods of different latitudes, which are independent of the viewing angle. The viewing angle will, however, lead to a second-order effect. Specifically, it will impact which latitude range is visible in the atmosphere. This, in turn, will weaken the power of the periodogram peaks that correspond to the less-visible latitudes, while leaving the power of other peaks essentially unchanged. Therefore, we conclude that our model’s fundamental predictions are mostly insensitive to the exact viewing angle (as long as they are generally close to equatorial, as expected for Luhman 16 B), and that slight inclinations would impact the periodogram power of the peaks but not their periods or the range of periods observed.

Thus, given the above, we conclude that, in spite of an important quantitative difference between the simulated Jupiter periodogram and the observed Luhman 16 B periodogram, there are multiple qualitative similarities that support the interpretation that Jupiter and Luhman 16 B share similarities between their atmospheric circulation. Therefore, we take the qualitative similarities as further support for the conclusion that Luhman 16 B displays zonal circulation. Figure 16 shows a qualitative sketch that illustrates the possible appearance of Luhman 16 B.

## 6.7. Long-timescale Evolution

Given our analysis, we have now viable models and explanations for the natures of the short-period peaks (*P* = 2.5 hr, *P* = 5.28 hr, *P* = 6.94 hr) in the Luhman 16 AB system. However, the TESS lightcurve also includes a long-period component, which appears in the Lomb–Scargle

periodogram as a *P* = 90.8 hr peak, although its exact period is not well determined. Our analysis of spacecraft positional jitter (Section 2.1) and background star variability (Section 2.2) strongly argue against such contamination and, therefore, for the genuine nature of the observed intensity variations. We will now briefly speculate about the possible nature of these variations.

The overall amplitude of the observed long-term modulations is about 10%, almost twice as large as the amplitude of the short-term modulations (mainly *P* = 5.28 hr). Given that we identified the rotational periods of both components A and B, and that these are 6.94 hr and 5.28 hr, respectively, we can establish that the long-period variability—regardless of its precise period—is more than an order of magnitude longer than the rotational periods. This observation means that the part of the atmosphere that is introducing the changes is visible throughout a rotational cycle.

Given these observations, we speculate that the long-period changes observed here emerge from one (or both) of the following sources: (a) global atmospheric evolution affecting both hemispheres in a coordinated manner; or (b) permanently visible polar regions.

With the information at hand we cannot distinguish between these two scenarios nor determine whether the long-period modulations emerge from components A or B (or both). Nevertheless, we briefly discuss both scenarios and the possible source. Scenario I—global evolution of atmospheric properties—is interesting as such a global evolution (e.g., cloud cover, atmospheric chemistry, or oscillations in the pressure–temperature profile) would seemingly naturally explain the long-period modulations. However, the ~90 hr timescale is likely too short to be explained by such a global evolution. Disturbances that propagate with the sound speed (~1 km s⁻¹) may circle a brown dwarf in ~100 hr. However, global changes are likely to propagate at significantly subsonic velocities, especially in the latitudinal direction, where—presumably—circulation zones must
## Sketch of the Possible Brightness Distribution and Zonal Circulation in Luhman 16B

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Sketch of the possible appearance of Luhman 16 B, based on the emerging evidence. Zonal circulation models and comparison to Jupiter suggests that low-latitude regions are dominated by the fastest jets, and that wind speeds at midlatitude are significantly lower. Circulation at the polar regions is likely to be vortex- and not jet-dominated. Cloud cover is likely to be correlated with the atmospheric circulation.</figcaption></figure>

be crossed. An example for long-period oscillations in brown dwarfs, as predicted by Showman et al. (2019), are quasi-biennial oscillations. Simulations predict that such oscillations would have a timescale of 1000–10,000 Earth days, i.e., about three orders of magnitude longer than the *P* ~ 90 hr period observed in our lightcurve. Therefore, we consider that global atmospheric evolution is an unlikely cause of the *P* ~ 90 hr lightcurve evolution.

In scenario II the permanently visible polar regions of the brown dwarf evolve. We note that no particular inclination is required for the polar regions to be permanently visible (as even an equatorial view will probe the polar regions), but brown dwarfs viewed more pole-on will also allow circumpolar regions to be observed. There are reasons to think that the polar regions will be morphologically distinct from the low- and mid-latitude regions: such a difference exists for both Saturn and Jupiter. Furthermore, it is predicted by three-dimensional circulation models for brown dwarfs, too: while the low- to mid-latitude regions are typically dominated by jets, the polar regions are dominated by wave dynamics and characterized, in appearance by vortices. The latitude of transition between the jet- and the vortex-dominated circulation depends primarily on the dynamical friction, wind speeds, and rotational speed. The polar regions are smaller in size than the two hemisphere areas scenario I required to alter; and the polar regions are not delineated by high-speed atmospheric jet boundaries. Therefore, naively, these regions appear to be more likely candidates for the source of the long-period evolution. If so, then the changes observed may correspond to rearrangements of vortices, i.e., possibly the evolution of large polar storms. We further speculate that, if indeed we are observing polar storm systems, these may more likely be observed on Luhman 16 A, which appears to be observed at a slightly inclined viewing angle (Section 6.2).

Although our study cannot determine the nature and source of the long-term evolution of the lightcurve, we note that the hypothesis that the long-term evolution probes permanently or quasi-permanently visible polar regions could be testable with future monitoring of Luhman 16 AB. In addition, a likely prediction of this hypothesis is an anticorrelation between short-term and long-term variability: sources seen closer to equatorial viewing angles may display greater rotational modulations, while sources viewed closer to the polar angles would be dominated by the slower-evolving polar storm

systems. If this hypothesis is confirmed, it will become possible to elevate our understanding of atmospheric circulation of brown dwarfs by studying both their jet-dominated zonal circulation and their vortex-dominated polar regions.

### 6.8. State of the Art and Outlook

The TESS lightcurve presented here—and potential future, even higher-quality data—will allow for detailed tests of brown dwarf circulation models. Most importantly, the quasi-continuous nature of the TESS lightcurve allows the study of the time evolution of the wave components present in a way no other existing data set does. The exploratory Spitzer monitoring by Apai et al. (2017) provided lightcurve segments that were too far apart to sample the lightcurve evolution between segments, revealing that the timescale of wave component evolution is shorter than the 20–100 rotations that separated the Spitzer lightcurve segments.

The data presented here strongly supports the general model put forward by Apai et al. (2017). The current TESS data set provides a relatively well-sampled, quasi-continuous lightcurve that enabled the identification of different circulation zones within Luhman 16 B. This data set sets the stage for direct tests of large-scale circulation models (e.g., Zhang & Showman 2014; Tan & Showman 2017, 2019) through the comparisons not of the lightcurves, but the periodogram generated from the simulated data and from the observations.

In Section 6.7 we also proposed the hypothesis that the long-term evolution of the lightcurve probes the evolution of the vortex-dominated polar regions. This is a testable hypothesis. If confirmed, the combination of the short-term and long-term variability would open a window to both the jet-dominated and the vortex-dominated regions of brown dwarf atmospheres.

In fact, future TESS observations have the potential to further refine the periods present in the lightcurve and also to constrain the long-term (~1000 rotations) evolution of the periodogram. Given that Luhman 16 is a uniquely bright variable brown dwarf, it is likely that it will remain the best target for future TESS observations. We note, furthermore, that to maximize the science return of the unique TESS data, simultaneous multi-wavelength lightcurves are required. Such multiwavelength data sets have been used to constrain, for example, the likely latitudinal width of the zonal circulation-driven banded structure (Apai et al. 2017) or, in Buenzli et al. (2012), simultaneous HST
and Spitzer observations were used to uncover large-scale latitudinal–vertical organization in a late T-type brown dwarf. A possible revisit of Luhman 16 by TESS will be a unique opportunity to obtain simultaneous multiwavelength observations (even if much shorter in coverage). Finally, we note that detailed and high-precision observations of the effective periods of the jet-dominated regions, such as those presented here, may be complemented by radio-wavelength period measurements that probe the rotation rate of the interior of the brown dwarfs. In a recent study, Allers et al. (2020) demonstrated the potential of comparing radio periods to near-infrared rotation periods and, from their difference, deduced an eastward wind speed of $650 \pm 310 \text{ m s}^{-1}$. The TESS observations presented here do not probe the rotation rate of the interior and, thus, lack an absolute reference frame for wind speed measurements. Complementing them in the future with radio observations will enhance their diagnostic power for studies of atmospheric dynamics in Luhman 16 B. This technique can of course also be extended for other, sufficiently bright, brown dwarfs.

## 7. Conclusions

Our study presents TESS lightcurves of Luhman 16 that cover about 22 days. The key findings of our study are as follows:

(1) The observations cover approximately 100 rotations of Luhman 16, during which the lightcurve remained non-constant, i.e., exhibiting significant rotational modulations.

(2) The lightcurve of Luhman 16 displays clear evolution during the observations: the lightcurve evolves rapidly and continuously, even over timescales of one to two rotational periods.

(3) We find peak-to-peak variability amplitudes over time-scale of a single assumed rotational period to be ~4% or less, but detect brightness differences as large as ~10% over longer baselines (100 hr).

(4) Based on Lomb–Scargle periodogram analysis we identify the rotational period of the component that dominates the short-term modulations as 5.28 hr. We also identify peaks at 2.5 and 6.94 hr, as well as a long-period peak with a period around 91 hr.

(5) Based on the periodogram analysis, we conclude that the lightcurve is very likely dominated by one of the binary components, but not both.

(6) We find strong evidence that the 2.5 and 5.28 hr peaks emerge from Luhman 16 B, while the 6.94 hr peak likely corresponds to the rotation period of Luhman 16 A. The 2.5 hr peaks likely correspond to $k = 2$ (half-period) waves, as seen in other brown dwarf lightcurves (Apai et al. 2017).

(7) Our analysis demonstrates that, in the periodogram, both the 2.5 hr and the 5.28 hr peaks consist of multiple adjacent peaks, a robust signature of the combined effects of differential rotation and planetary-scale waves. The fitted periods demonstrate an at least 20% rotational period range in Luhman 16 B, likely emerging from the combination of strong differential rotation and high-speed winds.

(8) We find that the derived rotational periods are fully consistent with observed $v \sin(i)$ values and that their combination suggests closely equatorial viewing angles (within 25° for A and within a few degrees for B).

(9) The lightcurve evolution observed in Luhman 16 is very well fit by the planetary-scale model proposed by Apai et al. (2017).

(10) We present a simple model for Jupiter, based on Cassini wind speed measurements, and show that the predicted periodogram is qualitatively similar to that of Luhman 16 B. This finding further supports the conclusion that Luhman 16 B's atmosphere is shaped by zonal circulation.

(11) We cannot determine whether the long-period (~90.8 hr) variability originates in Luhman 16 A or B, but propose that these modulations emerge from the (quasi-)permanently visible vortex-dominated polar regions. This is a testable hypothesis that may pave the way toward a more complete understanding of atmospheric circulation in brown dwarfs. A TESS revisit of the Luhman 16 system —currently scheduled for 2021 March–April—will provide a valuable data set to contrast to that presented here, allowing the comparison of the nature of the lightcurve evolution over a baseline of two years. Given our study and the recent results by Millar-Blanchaer et al. (2020), we can now conclude that both brown dwarfs in the Luhman 16 system display zonal circulation. Our study opens up the possibility, for the first time, to directly compare predictions of circulation models to observations via periodogram analysis of long-term monitoring data sets of brown dwarfs.

We are grateful to the referee whose suggestions have fundamentally changed and enhanced this study. L.R.B. acknowledges support by MIUR under PRIN program #2017Z2HSMF. D.N. acknowledges support from the French Centre National d'Etudes Spatiales (CNES). This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA Explorer Program. This study made use of data obtained in the Hubble Space Telescope programs GO-13748 and GO-14330. This work makes also use of results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia.

**Facilities:** TESS, Cassini, Hubble Space Telescope.

**Software:** astropy (Astropy Collaboration et al. 2013). IDL implementation of Genetic Algorithm (https://www.ncnr.nist.gov/staff/dimeo/idl_programs.html). Planetary-scale wave model (Apai et al. 2017). SWarp (Bertin et al. 2002). VARTOOLS (Hartman & Bakos 2016). ASTROMETRY. NET (Lang et al. 2010). AstroPy (Czesla et al. 2019).

## ORCID iDs

Dániel Apai https://orcid.org/0000-0003-3714-5855

Domenico Nardiello https://orcid.org/0000-0003-1149-3659

## References

Allers, K. N., Vos, J. M., Biller, B. A., & Williams, P. K. G. 2020, *Sci*, **368**, 169

Ammons, S. M., & Garcia, V. 2019, *AAS Meeting*, **233**, 114.07

Apai, D., Karalidi, T., Marley, M. S., et al. 2017, *Sci*, **357**, 683
Apai, D., Radigan, J., Buenzli, E., et al. 2013, ApJ, 768, 121

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Bedin, L. R., Pourbaix, D., Apai, D., et al. 2017, MNRAS, 470, 1140

Benatti, S., Nardiello, D., Malavolta, L., et al. 2019, A&A, 630, A81

Bertin, E., Mellier, Y., Radovich, M., et al. 2002, in ASP Conf. Ser. 281, The TERAPIX Pipeline, ed. D. A. Bohlender, D. Durand, & T. H. Handley (San Francisco, CA: ASP), 228

Biller, B. A., Crossfield, I. J. M., Mancini, L., et al. 2013, ApJL, 778, L10

Buenzli, E., Apai, D., Morley, C. V., et al. 2012, ApJL, 760, L31

Buenzli, E., Marley, M. S., Apai, D., et al. 2015a, ApJ, 812, 163

Buenzli, E., Saumon, D., Marley, M. S., et al. 2015b, ApJ, 798, 127

Burgasser, A. J., Gillon, M., Faherty, J. K., et al. 2014, ApJ, 785, 48

Burgasser, A. J., Sheppard, S. S., & Luhman, K. L. 2013, ApJ, 772, 129

Burrows, A., Hubbard, W. B., Lunine, J. I., & Liebert, J. 2001, RvMP, 73, 719

Crossfield, I. J. M., Biller, B., Schlieder, J. E., et al. 2014, Natur, 505, 654

Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python astronomy-related packages, Astrophysics Source Code Library, ascl:1906.010

Gaia Collaboration, Mignard, F., Klioner, S. A., et al. 2018, A&A, 616, A14

Gillon, M., Triaud, A. H. M. J., Jehin, E., et al. 2013, A&A, 555, L5

Hartman, J. D., & Bakos, G. Á. 2016, A&C, 17, 1

Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398

Imamura, T., Mitchell, J., Lebonnois, S., et al. 2020, SSRv, 216, 87

Karalidi, T., Apai, D., Marley, M. S., & Buenzli, E. 2016, ApJ, 825, 90

Kellogg, K., Metchev, S., Heinze, A., Gagné, J., & Kurtev, R. 2017, ApJ, 849, 72

Kniazev, A. Y., Vaisanen, P., Mužić, K., et al. 2013, ApJ, 770, 124

Lang, D., Hogg, D. W., Mierle, K., Blanton, M., & Roweis, S. 2010, AJ, 137, 1782

Libralato, M., Bedin, L. R., Nardiello, D., & Piotto, G. 2016a, MNRAS, 456, 1137

Libralato, M., Nardiello, D., Bedin, L. R., et al. 2016b, MNRAS, 463, 1780

Limaye, S. S. 1986, Icar, 65, 335

Lomb, N. R. 1976, Ap&SS, 39, 447

Luhman, K. L. 2013, ApJL, 767, L1

Mancini, L., Giacobbe, P., Littlefair, S. P., et al. 2015, A&A, 584, A104

Marois, C., Zuckerman, B., Konopacky, Q. M., Macintosh, B., & Barman, T. 2010, Natur, 468, 1080

Metchev, S. A., Heinze, A., Apai, D., et al. 2015, ApJ, 799, 154

Millar-Blanchaer, M. A., Girard, J. H., Karalidi, T., et al. 2020, ApJ, 894, 42

Nardiello, D., Bedin, L. R., Nascimbeni, V., et al. 2015, MNRAS, 447, 3536

Nardiello, D., Borsato, L., Piotto, G., et al. 2019, MNRAS, 490, 3806

Nardiello, D., Libralato, M., Bedin, L. R., et al. 2016a, MNRAS, 463, 1831

Nardiello, D., Libralato, M., Bedin, L. R., et al. 2016b, MNRAS, 455, 2337

Porco, C. C., West, R. A., McEwen, A., et al. 2003, Sci, 299, 1541

Radigan, J., Jayawardhana, R., Lafrenière, D., et al. 2012, ApJ, 750, 105

Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Proc. SPIE, 9143, 914320

Robinson, T. D., & Marley, M. S. 2014, ApJ, 785, 158

Scargle, J. D. 1982, ApJ, 263, 835

Showman, A. P., & Dowling, T. E. 2000, Sci, 289, 1737

Showman, A. P., & Kaspi, Y. 2012, arXiv:1210.7573

Showman, A. P., & Kaspi, Y. 2013, ApJ, 776, 85

Showman, A. P., Tan, X., & Parmentier, V. 2020, arXiv:2007.15363

Showman, A. P., Tan, X., & Zhang, X. 2019, ApJ, 883, 4

Simon, A. A., Rowe, J. F., Gaulme, P., et al. 2016, ApJ, 817, 162

Tan, X., & Showman, A. P. 2017, ApJ, 835, 186

Tan, X., & Showman, A. P. 2019, ApJ, 874, 111

VanderPlas, J. T. 2018, ApJS, 236, 16

Vos, J. M., Allers, K. N., & Biller, B. A. 2017, ApJ, 842, 78

Yang, H., Apai, D., Marley, M. S., et al. 2016, ApJ, 826, 8

Zechmeister, M., & Kürster, M. 2009, A&A, 496, 577

Zhang, X. 2020, RAA, 20, 099

Zhang, X., & Showman, A. P. 2014, ApJL, 788, L6