Latitude-dependent Atmospheric Waves and Long-period Modulations in Luhman 16 B from the Longest Light Curve of an Extrasolar World

In this work, we present the longest photometric monitoring of up to 1200 hr of the strongly variable brown dwarf binaries Luhman 16 AB and provide evidence of ±5% variability on a timescale of several to hundreds of hours for this object. We show that short-period rotational modulation around 5 hr (k = 1 wavenumber) and 2.5 hr (k = 2 wavenumber) dominate the variability under 10 hr, where the planetary-scale wave model composed of k = 1 and k = 2 waves provides good fits to both the periodograms and light curve. In particular, models consisting of three to four sine waves could explain the variability of the light-curve durations up to 100 hr. We show that the relative range of the k = 2 periods is narrower compared to the k = 1 periods. Using simple models of zonal banding in solar system giants, we suggest that the difference in period range arises from the difference in wind-speed distribution at low and mid-to-high latitudes in the atmosphere. Last, we show that Luhman 16 AB also exhibits long-period ±5% variability, with periods ranging from 15 hr up to 100 hr over the longest monitoring of this object. Our results for the k = 1 and k = 2 waves and long-period evolution are consistent with previous 3D atmosphere simulations, demonstrating that both latitude-dependent waves and slow-varying atmospheric features are potentially present in Luhman 16 AB atmospheres and are a significant contribution to the light-curve modulation over hundreds of rotations.

Prior to 2017, no long-term infrared monitoring data was available.Short-term (few hours-long) photometric light curves could be fitted well via elliptical spots (e.g., Apai et al. 2013;Karalidi et al. 2015).However, the long-term Spitzer monitoring -which became available in 2017 through a dedicated large Spitzer program -revealed complex and continuously evolving light curves (Apai et al. 2017a).The nature of the evolving light curves was inconsistent with and could not be fitted by elliptical spots: Rather, Apai et al. (2017a) showed that the light curves are very well reproduced by a new model: planetary-scale waves.In this model, planetaryscale waves trapped in zonal circulation modulate cloud thickness which, in turn, produces rotational modulation in the rotating atmospheres.The initial study, however, was limited to continuous observations of only four rotations, and only on three L/T dwarfs (although each object was visited eight times).
In 2021, Apai et al. (2021) presented a dataset that covered 20× more rotational periods continuously than any previous dataset.These TESS (Ricker et al. 2015) light curves probed rotational modulations in Luhman 16 AB -an L7.5+T0.5 binary brown dwarfs system that is also the closest to Earth.The Luhman 16AB light curve showed complex evolution over 100 rotations.These changes were successfully modeled with the planetary-scale wave model (Apai et al. 2021), just like previously the Spitzer light curves.However, the longer light curves allowed more comprehensive analysis and revealed the presence of not one, but multiple similar rotational periods.The authors attributed the range of periods observed to zonal circulation and differential rotation, i.e., wave-modulated cloud structures trapped in zones at different latitudes will have slightly different periods due to differential rotation.
Studies using other techniques also found evidence for zonal circulation and planetary-scale waves: Millar-Blanchaer et al. (2020) used VLT time-resolved polarization measurements to constrain the surface brightness distribution in the atmosphere of Luhman 16A and B. They found evidence for asymmetry (net polarization signal) which is consistent with the presence of bands and zones.Mukherjee et al. (2021) also reproduced the measured polarization of Luhman 16AB using polarization radiative transfer modeling that is coupled with general circulation model (GCM) outputs, and those GCM exhibits enhanced cloud abundances near the equatorial zones.Targeting another brown dwarf with high-amplitude rotational modulations, Zhou et al. (2022) also found that its light curve is evolving rapidly but not irregularly, and showed that the modulations can be well fit with the planetary-scale waves model.The presence of zonal circulation and planetary-scale waves should come as no surprise, as they are present in most, if not all, planetary atmospheres in the Solar System, including Earth and Jupiter (e.g., Fletcher et al. 2020).Zonal circulation and jets are predicted in rotation-dominated brown dwarf atmospheres (Zhang & Showman 2014;Showman et al. 2019;Tan 2022;Hammond et al. 2023) and planetary-scale waves are also seen in the most comprehensive general circulation models (e.g., Tan & Showman 2021) We note that one prominent measurement seemingly contradicts the presence of bands and zones in Luhman 16AB: Crossfield et al. (2014), which inverted time-resolved CO line profile modulations via a method developed for starspots, to identify the most likely surface brightness distribution at the CO-probed low-pressure region.The surface brightness model from this study does not show evidence for bands and zones -but due to its nature, the method is insensitive to such structures (see supplementary online material of Crossfield 2014).Karalidi et al. (2016) showed that the same data, when considering the uncertainties does, in fact, not contradict surface brightness models derived from other methods.
Arguably, important advances over the past ten years have been due to the improving temporal coverage and precision of time-resolved observations of rotational modulations.Partial-or single-rotation modulations do not provide enough constraints, due to the information loss inherent to hemisphere-integrated signal and the inverse problem of mapping exoplanets (e.g., Cowan et al. 2013).Modulations over four back-to-back rotations (in 32 epochs) provided strong evidence against elliptical spots playing dominant roles in shaping light curves (Apai et al. 2017a).A hundred-rotation coverage provided by TESS Apai et al. (2021) allowed initial characterization of the zones, wind speeds, and differential rotation in a brown dwarf atmosphere (Apai et al. 2021), and the tentative identification of k = 1 as well as k = 2 waves.
However, changes over very long timescales (thousands of rotations) and periods shorter than the rotational period could not be studied due to the lack of data.Our paper presents a new, longest monitoring dataset of brown dwarf atmospheres thus far that provides higher cadence data and allows detailed characterization of the k = 2 waves for the first time.Furthermore, we also offer very long baseline data, which enables the exploration of changes on timescales well exceed the rotational modulations -and therefore give access to a yet unexplored region of new atmospheric processes in brown dwarf atmospheres.
This paper extends the methodology of Apai et al. (2021) and uses higher-cadence and longer-baseline data to provide improved analysis of short-and long-period modulations.The paper is organized as follows: In Section 2, we discuss the TESS monitoring of Luhman 16 AB in Sectors 36-37 and strategies implemented to assess photometric contamination.In Section 3 we discuss our generalized Lomb-Scargle periodogram to explore the temporal structure and period components in the variability.In Section 4, we describe a multi-sine, planetary-scale waves model to explain the behavior of the light curve in short periods for periods under 10 hours.In Section 5, we will discuss the evolution of the long-period light curve from 15-100 hours.In Section 6, we discuss the interpretation of short-period components, provide a tentative comparison with the Solar System gas giants toy model, and discuss the long-period results in the context of 3D atmosphere studies.Finally, we summarize our key findings in Section 7.

TESS PHOTOMETRY
In this paper, we are presenting TESS (Ricker et al. 2015) photometry data (600-1,000 nm bandpass) of Luhman 16AB in Sectors 36 and 37. Photometry data is collected from March 07 2021 to April 28 2021 in two consecutive TESS orbits, covering a baseline of more than 50 days -about 48 of which are continuous science data not interrupted by data downlink gap.The downlink period totaled 2.12 days during two TESS orbits for sectors 36 and 37 1 .
We used the PATHOS pipeline from Nardiello et al. (2020) for the extraction and correction of the light curves from TESS full-frame images.PATHOS is a pipeline designed to extract high-precision photometric products for objects in crowded star fields using empirical PSFs and subtraction of neighbors.This helps to minimize potential flux contamination for faint objects.
The photometry presented in this work has a 10minute cadence compared to a cadence of 30 minutes in Sector 10 data from Apai et al. (2021).Figure 1 shows the full light curve of Luhman 16AB extracted with PATHOS with different photometric extraction apertures: the point-spread-function (PSF), and circular aperture with radii of 1, 2, and 3 pixels.Photometric noise becomes more significant as the aperture radius increases.In order to minimize photometric noise, we primarily use the PSF-aperture-extracted light curve throughout this study.
The average photometric error is 4.5%.The Earth and the Moon enter the field of view at the start of every orbit and introduce a significant amount of scattering photons.Background sources' scattered light is filtered out by removing photometric points with local sky noise factor σ SKY > 140 e s −1 and bad quality flag DQUALITY̸ = 0.This is in order to remove noisy photons from background sources, particularly from scattered lights of the Earth and Moon.
The temporal data itself is contaminated by various instrumental and astrophysical artifacts.We assessed possible sources of instrumental and astrophys-1 TESS Cycle 3 Data Release Notes ical contamination: A) the spacecraft positional fluctuation -how the spacecraft pointing fluctuates on a pixelby-pixel scale during observations; and B) background sources variations.The window function contamination, e.g., how the gaps in data sampling bias acquisition of certain periods of variability, will be discussed more in detail in Section 3.
The centroid of the PSF extraction on Luhman 16AB fluctuates on a pixel-by-pixel scale due to spacecraft jitters that introduce sources of periodic contaminants into data.The gradual, slight motion of the photometric aperture across the starfield might introduce a trend between pointing position and measure intensity.
Thus, we analyzed the time variation of these centroid positions and their distance to the actual source with periodograms, as shown in Figure 2. In order to quantify the relative strength of periodic components in the time series data, we use the Lomb-Scargle (GLS) periodogram method to transform time-series data into the frequency domain and explore their period distribution.
Slight spacecraft pointing drifts lead to variations in the celestial (X/Y) pointing positions and distance to the source.These variations are small compared to short-period brightness changes as shown in Figure 2.This is true for powers with periods below 25 hours but is not valid for powers with longer periods.The degree to which long-period powers in the data are affected is not uniform: gaps in the positional variation powers sometimes coincide with peaks in the data (i.e., around 72 and 90 hours).For periods larger than 125 hours, particularly from 150-200 hours, the periodogram powers of the XY positional variations are at least twice as large compared to powers at 100 hours.To avoid potential contamination, we take a 100-hour period to be the upper limit when assessing the long-period variation in the light curve in Section 5.
Another potential source of long-period contamination is the temporal variation of background sources.TESS has a pixel edge of 21 arcseconds -thus neighboring sources that fall into a region of ±1 pixel or 20-40 arcseconds (1-2 pixels) will potentially contaminate the signal.Past work of Apai et al. (2021) has shown that no bright stars are in the vicinity of Luhman 16 And that the strongest signal with the highest amplitude comes from the aperture centered around Luhman 16.They used deep HST multi-epoch photometry to assess the brightnesses, variability amplitudes, and distances of background stars around Luhman 16AB.By combining 12 HST epochs data, they showed that none introduced amplitude similar to that of Luhman 16AB to the TESS data.The GLS (Generalized Lomb-Scargle) periodograms of the TESS data (blue, first panel) and the instrumental pointing accuracy: X & Y position through time (green, second and third panel), and drift distance from the median center pixel through time (purple, fourth panel).The data periodogram shows strong power under 10 hours, which does not appear in the instrumental pointing over time.Peaks in the data curve which correspond to troughs in the spacecraft positional variation curves suggest that this period range is likely not contaminated, and vice versa.The green dashed line indicated the false-alarmprobability (FAP) level with a 10% power threshold, where peaks under the FAP level are potentially spurious.

PERIODOGRAM ANALYSIS
The generalized Lomb-Scargle (GLS) periodogram is a method to efficiently calculate the Fourier power spectrum estimator using unevenly sampled data, providing a way to understand the underlying periods of wavelike signals.Here, we used the GLS periodograms to explore the contribution of each periodic component to the overall power spectrum of the data.In the following sections, we compare the power spectrum of all timevarying components, including the reduced photometric data, the window function (which describes data collection windows), and Luhman 16AB coordinates given in the TESS detector pixel coordinates.We then generated synthetic fits for the periodogram of the data to understand the makeup of its variable components.

Brown dwarfs data
Rotational Period of Luhman 16A? Rotational Period of Luhman 16B k=2 wave?
Luhman 16AB: 0-20 hr Figure 3.The GLS periodograms of the window function describing gaps in data collection (red) and the data (blue).Two upper panels are for the period range from 0-200 hours.Two lower panels are for the period range from 0-20 hours.Little power contamination to the data exists for periods shorter than 20 hours while contamination is significant for longer periods.Inset texts indicate the literature value for the rotational periods of Luhman 16 A (≈ 7-hour) and B (≈ 5-hour), as well as the potential k = 2 wavenumber periods around 2.5 hours.Peaks in the data curve which correspond to troughs in the window function curves suggest that the corresponding period range is not contaminated, and vice versa.The green dashed line indicated the false-alarm-probability level with a 10% power threshold.

Comparing time-varying components with the GLS periodograms
Due to the nature of our long baseline data extending to about 50 days (∼ 1200 hours), variability may be seen from short (under 20 hours) to longer timescales (up to 200 hours).We first created GLS periodograms to compare the long-period variability of the light curve with the positional variation of TESS, which is shown in Figure 2. To distinguish between real and spurious peaks, we added the false alarm probability (FAP) with a power threshold of 10%.The same threshold of FAP is applied through the data periodograms in Figure 2  Periodogram fits for Luhman 16 B rotational period around 5 hours using different multi-sinusoidal models.For frequency-domain analysis, phase information is not needed, only the periods and amplitude.For each model plot, the best-fit periods are notated in the legend.One-sine, three-sine, and six-sine models (green curves) were used to fit the periodogram of the data (blue curve) following Apai et al. (2021).This periodogram fits show that a three-sine model is the simplest model that captures the highest amplitude structures (at periods of 4.7-hour and 5.25-hour) around the Luhman 16 B rotational period of 5 hours.and 4. It can be seen that peaks around the rotational period of , as well as long-period peaks from 70-125 hours, are well above the FAP level.
The positional variations are an important metric to assess pointing accuracy and to rule out time-dependent photometric contamination coming from spacecraft jitters or position changes.Positional variations are analyzed using the X, and Y pixel coordinates of Luhman 16AB derived from the TESS full-frame image.The distance drift is calculated using the median coordinates of all frames.
The window function is an important metric to assess the reliability of the periodogram signal against potential sampling biases.To assess potential contamination from sampling, we created GLS periodograms to compare at the same time the periodicity of the light curve and the window function.In Figure 3, the GLS periodogram of the data and the window function are shown.Upper panels are for long periods of up to 200 hours and lower panels are for short periods of up to 20 hours.
To create the window function, we generated an evenly-spaced array with a resolution equal to the 10minute cadence of the TESS data.The array is defined to be 1 where data exists and 0 elsewhere.Generally, the window function shows four significant gaps of no data collection, corresponding to a total of four data down-link gaps in the two TESS orbits.The ability to capture robust periodogram power will be inhibited for periods similar to or larger than the timescale of these gaps in the window function.
3.2.Synthetic sine fits for the periodogram Apai et al. (2021) showed that the periodogram calculated from shorter TESS data segments (∼500 hours) can be well fit by the planetary-scale wave model (Apai et al. 2017a).We also test this model on our longer and higher-cadence data, following the synthetic periodogram fit outlined in Apai et al. (2021).Our goal was to find a linear combination of sine waves whose periodogram matches the observed pattern.We generated synthetic light curves composed of one-, three-and six-period multi-sine models in the form Σ i α i sin(2πω i t) and added random, uniformly distributed noise matching the amplitude of the noise in the TESS data.For each of these models, we searched the parameter space for a best-fit solution with Sequential Least-Squares Programming gradient fitting2 .For frequency-domain analysis, phase information is not needed and therefore not included.The fit results are shown in the second, third, and fourth panels of Figure 5.  ple timescales of periodicity corresponding to short and long-period variations.The top panel of Figure 2 shows that the largest contribution to Luhman 16 AB variability comes from periods under 10 hours and periods above 50 hours.There is little power between periods of 125 and 200 hours in the Luhman 16 AB data.The strongest powers are found at around 5 hoursa well-supported rotational period value of Luhman 16 B from numerous past observations (Apai et al. 2021;Buenzli et al. 2015;Gillon et al. 2013;Burgasser et al. 2013).In the HST analysis by Buenzli et al. (2015) and TESS analysis by Apai et al. (2021), a strong power peak around 7.5 hours has been attributed to Luhman 16 A rotational period.

Results for periodogram analysis
However, the rotational period of Luhman 16 A around 7.5 hours in our TESS data is significantly smaller in power relative to the ∼5-hour power of Luhman 16 B and is well below the FAP level, as shown in the bottom panel of Figure 3, and in Figure 4.The power around Luhman 16 A rotational period is also significantly smaller in our TESS Sector 36 & 37 data compared to Sector 10 data from Apai et al. (2021).Figure 4 shows four periodograms for every consecutive 13-day data collection segment that make up the overall 52-day dataset and compares the power amplitude of Luhman 16 A and B variability.The power in the periodogram close to the Luhman 16 A rotation period remains consistently very low -nearly flat -through the duration of the TESS observations.
The cause for the lower rotational modulation amplitude of Luhman 16 A in our 2021 dataset (compared to the 2019 dataset in Apai et al. 2021) is not known.Long-term changes in Luhman 16 A atmosphere could be responsible for the difference in the photometry taken two years apart, as shown by Bedin et al. (2017).Due to the very low power around the rotation period of component A, in this study, we attribute most of the binaries' variability to the more dominant component B.
Strong power in window function: The window function periodogram, similar to the positional variation, also shows significant power for periods from 50 hours to 200 hours, and very little power for periods under 25 hours, as shown in the first and second panels of Figure 3.The short-period variation is thus significantly less contaminated compared to the longer-period data and could be a safe domain for analysis.
Significant long-period power in spacecraft pointing variation: The spacecraft pointing variation is evaluated using the X position, Y position, and me-dian variation of the position of the target on the detector.
The second, third, and fourth panels of Figure 2, show that the power of these positional variations only exists for periods longer than 25 hours.Some of these power peaks coincide with power peaks from the data periodograms (i.e., 63-hour peak), while others create 'gaps' with no significant power (i.e., 72 and 90-hour peaks).Lastly, there exist very large window function power peaks from 125 to 200 hours not seen in the data.
Periodograms analysis shown in Figures 2 and 3 show a dataset rich in periodicity ranging from periods under 10 hours to periods as long as 200 hours.At the same time, the spacecraft pointing variations and the window function also show a strong, non-uniform presence across a range of periods, which are potential contaminants in specific parts of the data but not others.
Here, the long-period analysis (above 50 hours) would require more careful biases assessment than the shortperiod (under 10 hours) variability.
Here, we discuss in detail the analysis of short-period variability in the periodogram of Luhman 16 B.
Periodogram fit using synthetic sine waves could explain the primary 5-hour power peak: Using the simple multi-sine model outlined in Section 3.2, we generated periodogram fit using synthetic light curves with the single-sine, three-sine, and six-sine models and generated GLS periodograms for each.The second, third, and fourth panels of Figure 5 are the periodograms of the corresponding synthetic light curves.
The three-sine and six-sine models capture the power peaks better than the single-sine model, and more sine waves brought the fits closer to the 1200-hour duration periodogram.The number of peaks between periods of 4-6 hours will decrease with shorter data duration from which the periodogram is generated, showing that the period distribution in the light curve changes and becomes more complex with time.Following Apai et al. (2021), we started with the three-sine model as a simple approach to fitting the light curve, which works well for segment up to 30-hour in duration.In section 4 we will discuss generating light curve fit in the time domain using a three-sine light curve model with phase offset information.
The primary structure with the strongest power around 5 hours seems to contain multiple sub-structures.The first panel of Figure 5 shows that the data periodogram contains many peaks between 4 and 6 hours.Two maxima are located near 4.7 hours and 5.25 hours respectively.A comparison with the 2019 data from Apai et al. (2021) shows that there are significantly more period peaks in the period range from 4 to 6 hours, due to the three times higher cadence and longer baseline of the new 2021 data set.
Power in k = 1 and k = 2 wavenumbers: We identify prominent peaks in the periodogram power around 2.5 hours, which is half of the rotational period (5-hour).The strongest periodogram peak is attributed to Luhman 16 B rotational period of around 5 hours (Apai et al. 2021;Buenzli et al. 2015).Following (Apai et al. 2021), we interpret the 2.5-hour peak group as half-periods, corresponding to k=2 wavenumbers, as in standing waves: λ ∝ 2L n where n is a natural number.We examined in detail the primary 5-hour power peak and the half-period 2.5-hour peak in Figure 6a.Here onward we refer to the primary peaks around the 4-6 hours range as the k = 1 periods and the secondary peaks around the 2-3 hours range as the k = 2 periods.
The period distributions of k = 1 and k = 2 periods both contain dominant double-peak features and multiple smaller peaks.These periods are centered at 2.5 hours and 5 hours respectively, suggesting that they are multiples of each other.
Figure 6a shows that the k = 2 periods have a relative range smaller than that of the k = 1 periods: 10.5% versus 22.2%.The relative period range is defined to be (P max − P min )/P min , where P max and P min is 4.5-5.5 hours and 2.375-2.625hours for the k = 1 and k = 2 periods, respectively.Upon closer inspection, the locations of the largest k = 2 peaks are not multiples of 2 of the largest k = 1 peaks.For example, multiplying the peaks at 2.43 and 2.575 of k = 2 periods does not reproduce the location of the peaks at 4.70 and 5.25 hours of k = 1 periods.
It is possible that the half-period k = 2 peaks represent the higher-order wavenumbers propagating zonally in the atmospheres.The k = 2 wavenumbers have been previously observed in Luhman 16 B (Apai et al. 2021) Multi-peaked structures both appeared clearly around the half-period and primary period of Neptune's rota- tion, i.e. 9 and 18 hours.However, compared to Luhman 16 B, the k = 1 and k = 2 periods in Neptune's power spectra are much more evenly matched -their relative range remained the same.The k = 1 and k = 2 periods of the Neptune data are almost perfectly in a one-totwo relationship, in contrast to Luhman 16 B.The exact mechanism to explain this is not well-known, but it is possible that the Neptune variability is mainly due to high-atmosphere cloud-top ices, intensified by solar UV radiation during solar maximum (Simon et al. 2016) that contributed to the variability observed with Kepler/K2.Zonal circulation is highlighted by ice particles as tracers of wind and results in the variability structure over the 49-day observation of Neptune.If deeper zonal circulation lower in the atmosphere had been visible, it would likely introduce more varied windspeed distribution which would not result in perfect one-to-two relationship between k = 1 and k = 2 periods as observed in the Neptune data.
Further analysis of this correspondence will be presented in Section 6.3, where we will examine Solar System gas giants windspeed distribution as a function of latitude.

SHORT-PERIOD LIGHT CURVE FITS
In the following analysis, we refer to "short periods" when speaking about variability with periods shorter than 10 hours.The previous periodogram analysis shows that 1) injecting multiple sine waves with periods close to the rotational period of the object could explain the periodogram, and 2) there exists multiple short and long-period variations in the data.The intensity variations on timescales shorter and longer than rotational periods are likely to have different physical origins; this is a motivation to examine them separately.In the following section, we will explore the data for short periods on par with the rotational period and attempt to fit the short periods with a light curve model composed of planetary-scale waves.

Fitting the Short-period light curve with a multi-sine model
To capture only the short-period peaks, a two-step process is employed: filtering the peaks with periods longer than 20 hours with a box-car average smoothing algorithm, and then subtracting the long-period light curve from the original data.
We selected three segments in the light curve for light curve fitting.We identified three different segments about 20 to 30 hours in duration.These segments correspond to 30-50, 183-203, and 409-433 hours, counting from the beginning of data collection on BTJ day 2280.30 in Sector 36.We also selected a longer segment at 410-490 hours with a duration of 80 hours.
For the 20 to 30 hours segments, we fitted the normalized light curve with a model comprising three sine waves in the form: with a total of 9 parameters for fitting, following Apai et al. (2021).
The three-sine model is the simplest initial model to explore the validity of planetary-scale waves.For the The fitting process is comprised of two steps: 1) finding an initial guess and 2) MCMC (Markov-Chain Monte Carlo) fitting.Firstly, we find initial guesses for the MCMC model from a fast decision tree optimizer hyperOPT (Bergstra et al. 2015).We found that fitting methods based on gradient descent do not work optimally with the sinusoidal model.Due to the oscillatory nature of sinusoids, these fits tend to converge at a local minimum.The decision tree optimizer mitigates this problem by exploring the parameter space simultaneously to find a better parameter set that will minimize the difference between the model and data.This step ensures the best initial guesses are obtained, but also takes less time to find best-fit parameters.
Secondly, we used MCMC to sample and explore the parameter space with emcee (Foreman- Mackey et al. 2013).From the initial guess given in the previous step, 'random walks' or random variations in parameters are taken to explore the parameter space and find the best fit.Weights can be given to each parameter to decide their 'random walk' ranges -well-constrained parameters need not deviate too much while less-constrained parameters could benefit from more varied guesses.

Results from Light-curve Fitting
Three-sine model fits the light curve well for the duration of ∼30 hours: We show in Section 3.3 that the light curve could be fitted using combinations of period components as a multi-sine model in the frequency space.To demonstrate this in the temporal space, we created test segments with a duration of 20 to 30 hours and fitted them with the multi-sine model.The resulting light-curve fits for the 30-50 hour, 183-203 hour and 409-433 hour segments are shown in Figures 7, 8, and 9, respectively.All three segments are fitted with three-sine wave models.The fit residuals are also shown within the bottom panel of each plot.
Results of the light curve fitting for the first segment within 30-50 hours are shown in Figure 7.The threesines model shows very good correspondence with the light curve data, with two periods around 5 hours and one period around 2.5 hours that best match the data.This period distribution reflects the periodogram analysis in Section 3.3, where the highest-power peaks for periods under 20 hours both correspond to the k = 1 and k = 2 waves, respectively.
Our analysis of the other two segments showed similarly consistent results.Figure 8 and Figure 9 show the light curve fits for the 180-210 hours segment and the 409-433 hours segment, respectively.Particularly, the best-fit periods correspond well to the peaks in the periodograms that are identified as k = 1 and k = 2 waves in Section 3.3.
We calculate the standard deviation of the fit residual (data minus fit) to estimate the goodness of the model.The respective standard deviation of fit residual for each model in Figure 7, 8, 9 are 1.02%, 1.17%, 1.14% respectively.These small residual shows that the sine-wave models well-captured all the components present in our light curve data on a timescale of 20-30 hours.The fit residuals for the light curve fit with our sine wave model are consistent with past results of Apai et al. (2021), which also found period components about 2.5 and 5 hours via the multi-sine model to fit the light curve best.Nominal photometric noises in the data are similar (≈ 4%), and standard deviations of the fit residual are also similar (≈ 1%).
Parameters in fit show small variation from one test segment to another: Throughout each segment shown in our analysis, the amplitudes in the multi-sine model are all non-constant as we move from one segment to another.This suggests that the amplitudes of period components experience a time evolution, pointing to a dynamic picture in the data.Moreover, the phase offset parameters also experience changes throughout the segment fits although these variations seem to be within a given bound.The phase parameters ϕ i (in radians) for each segment in Figure 7 Segment: 183 -203 hours | Amplitude 1, 2, 3 = 2.7, 1.6, 1.9 % | Phase 1, 2, 3 = 2.82, 0.21, 2.80 radian | Period P 1, 2, 3 = 4.9, 4.9, 2.5 hours | 2 = 0.2675 1.49), and (−2.81, −0.21, 2.80), and (0.39, 1.90, 2.72).
Four-sine model fits the light curve well for the duration of 80 hours: From Figure 2 and 3, the periodogram shows that the period distribution becomes more complex as one considers a longer duration of light curve evolution, which points to a non-static, actively evolving picture of the atmosphere.In the 410-490 hours segment -the longest segment we have in this analysis -we used a four-sine model to try to fit the light curve as having one more sine wave helps produce better fits.
Moving from a 20-hour to an 80-hour window, the amount of photometric data quadrupled for this test segment compared to the rest.We binned the light curve down from a 10-minute cadence to a 20-minute cadence after having applied a box-car median filter.The result for the 80-hour fit is shown in Figure 10.
The four-sine model fitted the light curve well and resulted in a 1.14% standard deviation of the residual, with a reduced chi-square χ 2 ν ∼ 0.74.Notably, the solution converged on a range of period components very similar to the primary and secondary power peaks in the periodogram: 2.4, 2.6 hours (k = 2 waves), and 4.8 and 5.4 hours (k = 1 waves).This correspondence is consistent with our analysis in Section 3, and also with the periods found for shorter segments with the three-sine model.The residual of the 80-hour fit is comparatively more varied than the shorter segment fits, hinting that remnants of periodic structures seem to persist within the fit.
We find that the four-sine model matched well the lightcurve behavior from 410-490 hours, better than the three-sine model.In Appendix B, we discuss the quality difference between the three-sine model and the four-sine model fit for this 80-hour data segment.We also used a Bayesian Information Criterion score (BIC) Bauldry (2015) to quantify the goodness of fit as well as over-fitting risk.Through that analysis, we found that the four-sine model produced a better fit compared to the three-sine model while remaining a simple enough model to explain the light curve behavior for this particular data segment.

LONG-PERIOD LIGHT CURVE EVOLUTION
Luhman 16 AB has been shown to exhibit variable behavior on timescales longer than 20 hours (Apai et al. 2021).In this section, we present the analysis of longperiod variability components on Luhman 16 AB tens of times the rotational period.Long-period variations in the light curve are responsible for brightness variability of about ±5% in the light curve, similar to the shortperiod variation.In the following analysis, we refer to "long periods" when discussing variability with periods longer than 15 up to 100 hours.

Frequency filtering of long-period variations
Before we proceeded with the analysis, we assessed to understand the potential contamination of background sources to the long-period data.We obtained a list of sources within the vicinity of 2 arc-minutes of Luhman 16 along with the periodogram of their PATHOSextracted light curves.We find that the sources bright and variable enough to present significant contamination comparable to Luhman 16 AB are all separated by more than 200 arcseconds (about 10 pixels).Sources under 200 arcseconds separation otherwise have featureless periodograms and no significant variability (see Appendix).Looking at the periodograms of the background sources' light curves, no variability is found in the long periods above 40 hours to under 100 hours, and the periodogram appears featureless.In agreement with the analysis conducted with background sources using HST from Apai et al. (2021), this indicates that background source contamination in the vicinity of Luhman 16AB is of minor concern for long-period data under 100 hours.
The periodogram analysis in the first panel of Figure 3 shows that the long periods from 50 to 125 hours contain the highest amplitude power, comparable to the amplitude of the short periods around 5 hours.However, the long period range coincides with strong contamination from the window function, as evident in the second panel of Figure 3.
Considering the window function contamination, there is likely a significant spurious signal in the periodogram above 150 hours and relatively significant power contamination for periods between 50 to 125 hours.However, there are 'gaps' of periods retrievable in the light curve data, and this range is also where the data shows the strongest power in long-period variability.
For potential lightcurve signals from telescope positional drift, the same conclusion can be drawn for periods above 150 hours in Figure 2: There is a strong possibility for periods contamination at this range.
Thus, we considered only the 15-95 hours range for our long-period analysis since this range contains the highest power with relatively less contamination coming from the windowing function and potential pointing errors.
In order to filter out periods within 15-95 hours, we employed a Fourier bandpass filter process.First, we interpolated the processed light curve data to create a uniformly-space light curve.Then, we applied a Fourier transform to the light curve to the frequency space and applied a bandpass such that values outside the target period range are zero.Lastly, we applied an inverse-Fourier transform to obtain a time-domain, long-period light curve.The resulting bandpass-filtered light curve is shown in Figure 11, where the long-period light curve is plotted over the original data, demonstrating that this long-period extraction captures the mean of the data very well.

Analysis Results of the Long-period Light Curve
The evolution of the Luhman 16 AB light curve is markedly different over long periods from its behavior over short periods.We tried fitting the multi-sine wave model and found that the model could not explain the long-period light curve in any segment longer than 100 hours, hinting that the long-period light curve is not time-stationary (i.e. a time-stationary series has statistical properties (e.g., mean and variance) that do not vary in time: for example, a sine wave with con-  First, we compare the flux values below and above the equilibrium level of 1.By assuming a baseline level of flux and evaluating flux increase and decrease around the baseline, we can identify excess flux that indicates potential atmospheric features, for example, storm spot that produces excess flux on a timescale close to selfrotation.
Figure 12 is the comparison of flux value above and below the equilibrium value of 1, in which flux below 1 is vertically mirrored.Numerically integrating the flux, we find the values 5.4 and 4.8 for the flux below 1 and flux above 1, respectively.Figure 12 shows that a dip is usually followed by a peak in flux such that the light curve appears relatively even across the baseline level.
Figure 13 shows a number of shorter segments demonstrating that the long-period filter captures the mean of the data very well.We note a strong correspondence in the long-period evolution with GCM (general circulation model) results from Tan & Showman (2021) where different inclined viewing angles for a 3D atmosphere produce different short and long-period evolution.
These timescales of variation are markedly very different from the short-period modulation characterized in the previous sections, hinting at their different physical origin.With this data set, we are demonstrating for the first time that Luhman 16 AB is continuously variable in periods up to 95 hours over a baseline of 1200 hours.In Section 6.7 and Section 6.8 we will discuss the potential origins of long-period atmospheric evolution and relevant results from GCM studies.In this section, we discuss the results obtained from the periodogram analysis showing strong power around the rotational period of Luhman 16 B (Section 3.3) and characteristics of the narrower range of the k = 2 waves with respect to the k = 1 waves.We show via a toy model of Jupiter and Saturn that the narrowing of the period range might be related to the different windspeed distribution between the equatorial region and the midto-high latitudes region.Then, we discuss and review possible mechanisms in the literature for the long-period light curve observed in Luhman 16 AB (Section 5.2).

Periodograms
Our periodogram analysis of Luhman 16 AB in Sectors 36 & 37 of TESS monitoring over a time span of 1,200 hours -covering more than two hundred rotations of the targets -identified strong peaks around the rotational period of Luhman 16B.In the range between 4.7 to 5.3 hours in our periodogram, there is strong evidence for a rich, multi-peaked distribution of period components (see Figure 3 and 6).Via our analysis in Sections 2 & 3.3 of spacecraft pointing drift and temporal sampling bias (modeled via a window function), we conclude that periods under 10 hours are unlikely to be contaminated, and emerge from Luhman 16B.This result agrees with the findings of Apai et al. (2021), in which variability at the same periods was also attributed to Luhman 16 B. Moreover, as in Apai et al. (2021), our study also supports the applicability of using the planetary-scale wave model as fits to the periodogram.
Sectors 36 and 37 TESS data have a cadence of 10 minutes, which is 3 times higher than the 30-minute cadence data used in Apai et al. (2021).The higher cadence, in combination with the longer baseline, reveals more sub-structures in the periodogram of Luhman 16 B.A multi-peaked distribution of periods with peaks ranging from 4.50 to 5.50 hours is identified (Figure 6).This period distribution is an interesting correspondence to the Solar System's giant planets, for example, Jupiter (Figure 15 of Apai et al. 2021) and Neptune (Figure 3B of Apai et al. 2017a), which both have multi-peaked power spectra around the rotation periods.

Persistent Period Distribution
Our new Sector 36 and 37 data show a periodogram for Luhman 16 B that is consistent in all of its properties with the earlier datasets from 2019, presented in Apai et al. (2021).This is an important finding as it demonstrates that the periodicities present in the targets are sustained over more than three thousand of rotations.Although past data are taken prior to TESS (e.g., Gillon et al. 2013;Buenzli et al. 2015;Karalidi et al. 2016) provide snapshots in time rather than multi-period light curves, we note that the qualitative behavior of Luhman 16 B described in those studies is consistent with those found here.Thus, there is no evidence for a change in the nature of the rotational modulation and the evidence available is consistent with Luhman 16 B displaying rotational modulations of the same nature since its discovery in 2013 (Luhman 2013), over the past 10 years.
Therefore, we conclude that the atmospheric processes that modulate the brightness distribution in the atmosphere of Luhman 16 B remained active and dominant over the timescales of at least thousands of rotations (over the course of the available TESS observations between 2019 and 2022), and likely well over 10,000 rotations (since 2013).This persistent nature of the modulations demonstrates that the nature of atmospheric brightness distribution is not due to a transient state, such as a chance cloud alignment, but tied to a fundamental, powerful mechanism within the atmosphere that does not change significantly even over long timelines (>10 years).The most obvious source of such a mechanism is atmospheric circulation, which we will explore in the following section.

k=2 Wavenumber
In this work, we confirmed the existence of group of periods half of the rotational period of Luhman 16 B in the periodograms, which we identify as k = 2 wavenumber waves.We find strong similarity between the periodogram power structure of the primary k = 1 period group and the k = 2 period group: They are similar in their multi-peaked nature and approximately bimodal distributions of periodogram power (see Figure 6).The presence of k = 2 periods in our data is consistent with results from two past studies: In Apai et al. (2017a), a periodicity analysis of Spitzer light curves of two L/T transition brown dwarfs showed evidence for peaks at half the rotational period.The TESS light curves of Luhman 16 B, presented in Apai et al. (2021), displayed strong peaks consistent with k = 2 waves in the generalized Lomb-Scargle periodogram for Sector 10 TESS data.
The light curves of Solar System gas giants, such as Jupiter and Neptune, all display strong power around the rotational period.For example, the effective period distribution of winds in Jupiter reveals differential rotation (Porco et al. 2003, Fletcher et al. 2020).Apai et al. 2017b demonstrated that periodogram peaks at half the rotational period also exist in the power spectrum of the Neptune (Simon et al. 2016) in ways similar to their brown dwarf counterparts.The general circulation models (GCMs) of Jupiter and Saturn also display prominent banded structures (Schneider & Liu 2009;Lian & Showman 2010;Young et al. 2019;Spiga et al. 2020), similar to banded structures seen in brown dwarf GCMs (Showman et al. 2019;Tan 2022).Apai et al. (2017a) showed that the peaks around the half rotational period are well fit via the planetary-scale waves model: k = 2 waves with periods of half of that of the primary k = 1 waves.The planetary-scale waves modulate the brightness of zones and belts which, in turn, are rotating at different rotational periods due to a combination of zonal circulation and differential rotation.The planetary-scale waves could modulate the thickness and/or structure of condensate clouds (e.g., Apai et al. 2013;Vos et al. 2022).
In the general circulation models (GCMs) of Tan & Showman (2021) with periods of 2.5, 5, and 20 hours, the synthetic periodogram displays large band structures in the equator as a result of strong circulation.It is notable that these models also contained k = 2 periods.Weak drag seemed to be an important precursor for the formation of large-scale jets: in Tan ( 2022), the weak-drag atmosphere model (characteristic drag timescale τ drag = 10 7 s) contains much larger k = 2 period power compared to strong-drag atmospheres (τ drag = 10 5 s).Tan (2022) postulated that two planetary-scale waves with the same zonal frequencies but different zonal phases could result in an apparent feature of zonal wavenumber 2 within the power spectrum.
However, via a synthetic light curve under the assumption of planetary-scale waves, we show that a k = 2 period component must correspond to an actual wave in the atmosphere and is not an effect of aliasing.Figure 14 shows a simulated light curve composed of two closely-spaced period components (4.7 and 5.25 hours) with a phase shift of 1.36 radians, with uniformed noise added.The figure shows that -unlike as it was assumed in some past works -the two waves do not lead to k = 2 period component.We repeated this analysis for twoperiod components with the same value for around 5 hours and cycled through a number of phase shift values.We found that all the power at 2.5 hours was at least 2 orders of magnitude smaller than the power at the primary period of 5 hours.The k = 2 periods suggest that Luhman 16 AB would thus have waves at half the period of k = 1 waves.

Planetary-scale Waves vs. Fourier Series
With k=1 and k=2 waves identified in our analysis, it may be relevant to contrast the planetary-scale waves model with the Fourier series.Superficially, it may appear that the success of the planetary-scale wave model is due to its mathematical similarity to the Fourier series.However, upon closer inspection, the two are substantially different, as explained below.
Despite the similarities, the planetary-scale wave model is mathematically different from a Fourier series expansion in two important ways.First, in the Fourier series, the function is expanded as a sum of waves with decreasing periods that form a series (P , P 2 , P 3 ,..., P n ).In essence, as the number of series elements increases, the frequency space (periodogram) is increasingly covered.However, the planetary-scale wave model is different: Periods are not fixed absolute or relative values but are mostly unconstrained.They converge to two groups of periods (group of periods around k=1 and k=2), rather than to a Fourier series.This is significant because these waves do not cover a broad frequency space, but the opposite: They concentrate into 1 or 2 groups.
The second key difference we will highlight is that Fourier expansion is applicable to periodic functions and the fundamental mode of the series (the longest period) is set to equal the entire length of the data.Therefore, the Fourier series is not used or should be expected to reproduce aperiodic functions beyond the extent of the  fundamental mode (longest wavelength).This, in our case, for example, a Fourier expansion's periods would be 1, 200 h, 600 h, 400 h, 240 h, 200 h, etc.
In contrast, our planetary-scale wave models converge on a sum of sine waves with periods grouping around 2.5 and 5.2 h, i.e., about 200 times shorter periods than the fundamental modes of the Fourier series would be.
In short, the planetary-scale wave model fits the lightcurve with a sum of sine waves with periods as free parameters, which is similar but mathematically and functionally different from a series of sine waves with periods that decrease as P n .

Models of Gas Giant Windspeed Distribution
Returning to Solar System gas giants, the comparison of Luhman 16 AB and Neptune's k = 1 and k = 2 highlighted an interesting difference: The period multiples of the strongest peak in the period distribution of Luhman 16 AB do not match, whereas these period multiples match relatively well in the Neptune data (Si-mon et al. 2016).In the next section, we will explore this difference.Both Neptune's and Uranus' zonal jets are not strongly visible in wavelengths available to Voyager (see comprehensive review by Sánchez-Lavega et al. 2019, Fletcher et al. 2020), although they show more structures and zonal banding in the NIR wavelength at epochs where cloud-top variability is greatest, i.e.Neptune (Irwin et al. 2023); Uranus with Keck H-band & Voyager 2 joint reanalysis (Sromovsky et al. 2015).Most of Neptune's deeper tropospheric zonal circulation is not visible in Kepler photometry.The variability in Simon et al. (2016) is dominated by the top-of-atmospheres cloud bands and not from waves coming from deeper zonal bands.
In the case of Luhman 16 B, the difference in period multiples of k = 1 and k = 2 waves hint at two different ranges of period distribution which may contain spatial information.In order to explore why the period distribution may be narrower for the k = 2 period group, we will use a simple model of Jupiter's period distribution as a function of latitude.
In our model of Jupiter, following Apai et al. (2021), in order to calculate the effective period of winds as a function of latitude, we used published windspeed measurements of Jupiter from the Cassini mission taken from Porco et al. (2003).In a fly-by in 2003, Cassini took multiple images of Jupiter over one rotational period and tracked the clouds to measure windspeed.The revolution time for each latitude is the effective period at that latitude.
From the wind speed measurements as a function of latitudes, we generated a synthetic model composed of the period-band structure.The effective period P eff is related to the latitudinal windspeed by: where L = 2πr cos(θ) is the circumference at each latitude θ, the rotational speed u rot = L/P where P is Jupiter's self-rotation period, and u wind the measured windspeed at each latitude θ.Thus we obtain an equirectangular projection of the period-band, which could be mapped back onto a sphere.This is achieved via the 3D mesh method in the visualization package mayavi 3 (Ramachandran & Varoquaux 2011) using a spherical mesh function.The final image used for analysis is an actual 3D sphere that is projected back onto a plane.This implementation provides more flexibility 3 mayavi Spherical Projection.
to explore the impact of different inclinations compared to the 2D model in Apai et al. (2021).
Luhman 16 B is inclined at about 20 • (Apai et al. 2021).We included a similar inclination in the 3D model of Jupiter shown in Figure 15.A different inclination changes the fraction of area visible between the polar and the equatorial regions, thus changing the period distribution in the visible disk.
We hypothesize that the k=2 period group may be emerging from spatially different distribution than the k=1 period group does.To test this, we masked out different regions of latitude and calculated the periodogram that emerges from the non-blocked regions.The results of this analysis and period histograms are shown in Figure 15.
The highest and lowest wind speeds correspond to the shortest and longest effective periods, forming the tail ends of the full-disk (unmasked disk) period distribution.These shortest and longest period values are all located at the equatorial region within 35 • S to 35 • N. Hence, if we consider only the windspeed from the midlatitudes to the polar region, excluding the equator, the period distribution is narrower than that of the full (unmasked) disk.
To broaden our windspeed comparison beyond Jupiter, we created a second period-band model of Saturn, assuming a zero-inclination, circular disk -not taking into account the oblateness of Saturn.We utilized the Cassini windspeed data as a function of latitude on Saturn (averaged from 2004to 2009, García-Melendo et al. 2011) as shown in Figure 16.The Saturn measurements are averaged over a longer 5-year timescale.The resulting synthetic toy model at zero-inclination is shown in Figure 16.
The resulting period distribution of Saturn differs from Jupiter (Figure 15) in that Saturn's highest windspeeds (shortest periods) are all concentrated in the equatorial region, in contrast to the more varied windspeed values in the equatorial region of Jupiter.
In the case of Saturn, almost the entire shorter-period half of the distribution is cut off in the masked disk, in contrast with Jupiter where only the tail ends are cut off.Nonetheless, this shows that our simple models are consistent in one way: Both their period distribution shrinks in the period range once the equatorial region is discarded, showing that the mid-latitude-to-polar region has a smaller range of period distribution.
The narrowing of the period distributions at different latitudes in our Jupiter and Saturn model mirrors the difference in the period distributions of the k = 1 and k = 2 on Luhman 16 B. Thus, a natural explanation for the Luhman 16 B pattern is that the k = 1 and k = 2 waves arise from different latitudes, such that the k = 2 wavenumber arises from only the mid-latitude region.

Short-period Light Curve Fits
The fact that a simple model of multi-sine waves analogous to planetary-scale waves could fit the light curve over multiple timescales is a strong support of the rotation and 3D circulation driving zonal circulation on Luhman 16 AB.This result is in strong agreement with previous works (Apai et al. 2021(Apai et al. , 2017b) and shows the applicability of planetary-scale waves for certain brown dwarf atmospheres including Luhman 16 AB.
We showed that the planetary-scale wave model is a simple model that could explain rotational modulation in the data over tens of time the rotational period, i.e. up to 100 hours for Luhman 16 AB.

Long-Period Light Curve Evolution
Via applying a frequency bandpass to isolate longperiod varying components, we showed that the Luhman 16 AB light curve contains components that vary on timescales much longer than the rotational period, with periods up to 100 hours.In this section, we will discuss some potential mechanisms for long-period variability in brown dwarf atmospheres.
A review of theoretical models and observations from Showman et al. (2020) suggests that long-term atmospheric variation is a trademark feature of gas giants and brown dwarfs.Quasi-periodic oscillation in the equatorial jets on a multi-year timescale exists on both Jupiter and Saturn, where eastward and westward jets migrate in the latitudinal direction and change their positions over time.This quasi-periodic oscillation has a period of about 4 years for Jupiter and 15 years for Saturn.The 3D GCMs for gas giants' atmospheres (Showman et al. 2019) also captured these oscillations with periods ranging from a few years up to 12 years.
Motivated by the unexpected long-term observed variability of Jupiter by Orton et al. (2023), Hori et al. (2023) demonstrated theoretically that the torsional oscillations arising from the rapid rotation of ionized, conducting cylindrical layers in the deep interior of Jupiter can explain the timescale of cloud-level variability observed with Juno.This torsional force between the interior layers bends the magnetic fields of Jupiter and creates waves that propagate to the cloud top.This is a potential model to explain the long-period, multi-year variability seen in Jupiter's atmosphere.
There is a clear timescale difference between the longperiod oscillations in our Luhman 16 AB data, which go up to hundreds of hours, versus the multi-year timescale discussed in the literature.It could be speculated that the difference in mass, magnetic field strength, and the convective interior might play a role in timescale differences of long-term atmospheric variation between brown dwarfs and gas giants.Still, the nature of long-period variability remains inconclusive because there has not been a time-resolved, year-long monitoring dataset of brown dwarf atmosphere yet.Indication of long-term changes on brown dwarfs exists via comparing photometry taken several years apart (Luhman 16 B, Bedin et al. 2017), but it is unclear how these multi-year oscillations would arise on brown dwarf atmospheres and to what degree they are connected to shorter-period variation.In the near future, time-resolved extended multi-band monitoring will be the key to understanding the physics behind long-period atmospheric variation on brown dwarf atmospheres.6.8.Insights from GCMs for the explanation of short and long period light curve GCMs of brown dwarf atmospheres in Tan & Showman (2021) demonstrated that cloud radiative feedback generates east-west traveling large-scale waves and stochastic vortices.
With a weak basal drag which crudely represents interactions between the active weather layer and the deep quiescent interior, strong zonal jets can form and the jet speed depends on the drag strength.The inherent inhomogeneous cloud coverage naturally leads to rotational light curves.In addition, traveling waves and the statistical evolution of storms produce irregular features in the light curve which can broaden the power spectrum of the light curve around the rotation period; the presence of zonal k = 2 waves imprints strong k = 2 signatures in the light-curve power spectrum.
The modeled light curve and its power spectrum from a GCM in Tan & Showman (2021) are presented in Fig- ure 17.This particular model assumes a rotation period of 5 hours, similar to that of Luhman 16B, and a moderate basal drag with a drag timescale of τ drag = 10 6 s such that the model developed a broad westward equatorial jet and high-latitude eastward jets with a peak speed of ∼ 500 ms −2 (see Figure 9 in Tan & Showman 2021).The modeled light curve shows both high-frequency variations that are caused by the rotational modulation of the zonal k = 1 and k = 2 atmospheric features and lowfrequency variations that vary over timescales on the order of 100 hours.Powers near 5 hours are significantly broadened and discretized into multiple peaks that can extend to ±0.4 hours away from the major peak, quantitatively in good agreement with the observed power spectrum shown in Figure 6.The k = 2 power in the modeled light curve shows a narrower broadening than the k = 1 components, similar to the observed one; however, this model did not produce double k = 2 groups which is robustly seen in the observed one.The longer-timescale variability is from the statistical evolution of the ensemble storms, despite that individual storms could have a much shorter evolution timescale.This statistical fluctuation has been shown in box models that exclude equatorial waves and rotational modulations in Tan & Showman (2021).The modeled long-term variability is remarkably similar to the observed one for Luhman 16B despite that the model was not tuned for Luhman 16B.
Note that zonal jets in the model are not the main reason for the broadened power spectrum.The equatorial westward jet provides a redshift for the bulk power spectrum including the k = 1 and k = 2 components.High-latitude eastward jets contribute little to the power spectrum because low-latitude variability dominates the light curve.A similar model with a strong drag generates a similarly broadened power spectrum but no redshift; another similar model with an even weaker drag (therefore with a stronger westward equatorial jet) produces a similar power spectrum that is merely redshifted more.
Comparisons between the observed and modeled lightcurve power spectra indicate a somewhat alternative view for some features of the observed power spectrum of Luhman 16B: Luhman 16B may potentially develop a global storm system in the atmosphere that host selforganized equatorial waves and vortices driven primarily by the cloud radiative feedback; k = 1 light-curve components are broadened by the traveling waves and stochastically evolving storms, while longer statistically fluctuation of the storm system drive long-term lightcurve variability.
Critically, the presence of two k = 2 groups in the observed power spectrum (but lacks in the models) suggests the development of alternating jets on Luhman 16B that are capable of affecting the observed light curve.In the model framework of Tan & Showman (2021), this indicates that the alternating jets should form at lower latitudes where most variability is from, necessitating further modeling efforts to understand jet formation on Luhman 16B.

Potential Physical Interpretation
At this point, despite its success in explaining the available rotational modulations over many rotational periods, the planetary-scale wave model is a mathematical construct providing a simple phenomenological model.It is not a self-consistent physical model, although it is physically motivated and it aligns with the general predictions of detailed dynamical models (see Section 6.8).In the following, we will provide a potential physical interpretation without claiming this to be a complete (or even unique) interpretation.
Our interpretation emerges from a combination of four established facts: (1) Atmospheric waves are present in all atmospheres and their potential wavelengths range from molecular scales to planetary scales.(2) Rotationally-dominated, internally heated gaseous atmospheres will develop zonal circulation, differential rotation, and latitudinally varying wind speeds.(3) Condensation products (dust particles) form in brown dwarf atmospheres at locations where the pressuretemperature conditions lead to super-saturation.(4) The L/T spectral type transition in brown dwarfs is broadly interpreted to occur as the top of the sili- By combining the four facts above and following them to their conclusion, we provide the following potential physical interpretation for the planetary-scale waves in Luhman 16B: In this object, we observe the top of the cloud deck close to the observable optical depth.Smallscale turbulence and mixing drive waves of increasing wavelengths, all the way to the longest-possible wavelength standing waves with wavelength equaling the circumference of the object.As the circumference changes with latitude and the rotating gaseous also inhibits differential rotation, the longest wavelength and the rotational period corresponding to the longest waves will vary somewhat with latitude.Waves with vertical components will periodically move gas parcels up and down in the atmosphere.Such motion close to the saturation curve of the key condensates (silicates and iron) will result in periodic condensation/evaporation of these, similar to gravity-wave-induced spatially periodic water vapor clouds that are common in Earth's atmosphere.
Put together, the planetary-scale waves will introduce a spatially periodic modulation in cloud thickness, that will be observed to show a temporal periodicity with a range of periods similar to the overall average rotation period of the atmosphere, but showing variations due to differential rotation and wind speeds.Our data also shows very strong evidence for standing waves that correspond to half the circumference of the atmosphere, also showing multiple peaks such as the k=1 waves.
The above qualitative interpretation appears to match the body of evidence we are aware of, and aligns well with the general findings of state-of-the-art GCMs, without requiring an assumption of any individual process or phenomenon that is not already established.
Perhaps the only seemingly surprising aspect of this interpretation may be the assumption that the largescale waves coincide in their vertical location with the cloud top (i.e., pressure-temperature conditions that are required for super-saturation of silicates).Generally and in most atmospheres, there may not be an obvious physical reason for these two to spatially coincide.However, this apparent contradiction is readily resolved when one considers the target selection biases: Most rotational variability studies focus on L/T transition brown dwarfs because that is the spectral type range where the largestamplitude modulations are observed (e.g., Radigan 2014 -which observation then naturally aligns with the explanation above.(Luhman 16B is specifically targeted because it is one of the most variable brown dwarfs.)In other words, in atmospheres where the large-scale waves are well above or well beneath the condensate cloud deck's top, the waves will not modulate the cloud properties; not break the axisymmetry of the atmosphere, nor leading to large-amplitude rotational modulations.
While the above interpretation provides a compelling and qualitatively consistent picture, we caution that it is not yet a quantitative and self-consistent model.While individual model components exist and support the above picture, more work is needed to build a multicomponent, self-consistent model that can provide a quantitative model for the planetary-scale waves' impact on cloud thickness in L/T brown dwarfs.

CONCLUSIONS
Luhman 16 AB is the brightest brown dwarf binary observable from the Solar System and an excellent target for probing ultra-cool atmospheres.This observation helps to understand the broad parameter space of atmospheres to which ultra-cool atmospheres like those of Jovian planets and sub-Neptunian planets also belong.Below we summarized our findings of Luhman 16 AB photometric monitoring.

Longest photometric monitoring of Luhman 16
AB : With light curve data covering 1,200 hours (50 days) baseline with a 10-minute cadence, we found evidence for strong, persistent photometric variability with an amplitude of about ±5%.We found strong periodic modulations around the rotational period of Luhman 16 B centered at ∼5 hours and half the rotational period at ∼2.5 hours.We ruled out the likelihood of photometric contamination by spacecraft pointing drifts or background sources.We also ruled out contribution from Luhman 16 A around 7.5 hours, due to very low variability amplitude in our dataset around this period.
2. Complex structure in short-period (P<10 h) periodogram: Our generalized Lomb-Scargle periodogram analysis shows a range of amplitudes and periods, consistent with previous work (Apai et al. 2021).Short-period variability (P<6 h) originates from rotational modulations in the atmosphere of Luhman 16 B.We showed that: (a) Strong power exists at periods near the rotational period of Luhman 16 B: Both in the primary k = 1 periods in the range 4.25-5.75hours with peaks at 4.7 and 5.25 hours, and the k = 2 periods in the range 2.3-2.65 hours with peaks at 2.43 and 2.575 hours.
(b) A multi-sine model analogous to planetaryscale waves could fit the periodogram of Luhman 16 B centered at 5 hours; (c) There is a multi-peak patterns coinciding with the k = 2 (half rotational period) range.
The pattern is similar to the k = 1 periods.Via analysis of Jupiter and Saturn's effective period distribution and Neptune's power spectra, we show that the observed distribution is observed if the k=2 waves are mostly present in mid-to-high latitudes, while k=1 waves are also present in the equatorial regions, 3. The planetary-scale wave model explains the light curves: The sum of 3-4 sine waves provides a very good match to the light curves (1σ residuals ∼1%).
We found that the period distributions in our fits agreed well with the dominant components in the periodogram of Luhman 16 B.The periods converge around 2.5 hours and 5 hours for all light curve segments -similar to the k = 1 and k = 2 periods.Our fits also showed that the periods and amplitude of sine waves experience time evolution as we move from one segment to another.
4. Sustained Long-period variability: We found longperiod variability (different from rotational modulations) with an amplitude of about 5% (much longer than the ∼5.2 hour rotational period of Luhman 16 B).The long-period components likely originate from distinct regions in the atmosphere, independent of short-period zonal circulation.
5. Qualitative consistency with GCM results: We found that certain properties in the period dis-tribution of our observations are consistent with specific simulations of weak-drag, rapidly rotating atmospheres from Tan (2022) and Tan & Showman (2021).Comparing the data and simulations, both the k = 1 and k = 2 period components are found in the periodogram, as well as the narrowed period distribution between k = 1 and k = 2 periods.Long-period variation in the light curve of up to hundreds of hours is also found.
Our study demonstrates the power of very long-baseline photometric monitoring of brown dwarf and exoplanet atmospheres and opens a window on brown dwarf brightness evolution over tens and hundreds of rotations.Further high-cadence data will allow the k=2 waves described here to be studied in greater detail.(Apai et al. 2021).In the follow-ing section, we again show that it is unlikely Luhman 16 AB contains photometric contamination from background sources in our data -specifically when using the PSF-extracted light curve where the photometricextraction radii are less than 1 TESS-pixel (about 21 arcseconds).
In Figure 18, we plot background sources of comparable or brighter TESS-magnitude (magnitude for TESS bandpass from 600-1000nm) than Luhman 16 AB within 2 arc-minutes.In Figure 19, we plot the light curve and the resulting Lomb-Scargle periodograms for each of these objects.It can be seen that most objects shown do not contain strong period power under 5 hours where Luhman 16 AB variability is strongest.It is noteworthy that the periodogram power shown via the Lomb-Scargles periodogram does not represent absolute power but rather normalized power, which is important when comparing periodograms from two different objects.Thus, an inspection of the light curve is always necessary to assess whether actual variability exists in the data, not just solely from whether certain period power exists in the periodogram.

Figure 1 .
Figure 1.Significant variability exists in TESS 50-day light curve of Luhman 16AB in Sector 36 and 37.The cadence is 10 minutes.Gaps in data are due to downlink periods.Different curves show the result of aperture size adjustment from the PATHOS light curve extraction pipeline.The PSF-extracted light curve has the smallest photometric error.

Figure 5 .
Figure5.Periodogram fits for Luhman 16 B rotational period around 5 hours using different multi-sinusoidal models.For frequency-domain analysis, phase information is not needed, only the periods and amplitude.For each model plot, the best-fit periods are notated in the legend.One-sine, three-sine, and six-sine models (green curves) were used to fit the periodogram of the data (blue curve) followingApai et al. (2021).This periodogram fits show that a three-sine model is the simplest model that captures the highest amplitude structures (at periods of 4.7-hour and 5.25-hour) around the Luhman 16 B rotational period of 5 hours.

Figure 6 .
Figure 6.Comparing k = 1 and k = 2 wavenumbers on Luhman 16 AB periodograms and Neptune power spectra.The vertical lines represent the boundary used to calculate the relative range of period structures.Black, dashed lines represent multiple-of-two of the range; i.e.[4.5, 5.5]  hours are 2 times the [2.25, 2.75] hours range.The colored, dashed lines in each plot of Panel (a) represent the multiple-of-two locations of the peaks from the other plot; i.e. the k=1 plot contains location of k=2 peaks multiplied by 2, the k=2 plot contains location of k=1 peaks divided by 2. Panel (a) -Luhman 16 B data periodograms both contain multi-peaked period components at 4 to 6 hours (k = 1 wavenumber, blue) and 2 to 3 hours (k = 2, red).The period multiples do not match: the k = 1 largest peaks do not correspond to the k = 2 largest peaks as shown by the dotted lines.Inset plots show the relative powers and positions of two structures.The relative period range is defined as (Pmax − Pmin)/Pmin.Panel (b) -Neptune power spectra fromSimon et al. (2016) showing the power of each period component arising to rotational modulation of high-altitude cloud bands.There exists matching similarity in the largest peaks at 16-19 hours (k = 1 waves) and at 8.0-9.5 hours (k = 2 waves).
, in two other L/T transition brown dwarfs Apai et al. 2017a and also in the power spectrum of Neptune. Figure 6b shows the Neptune power spectra taken from the Kepler/K2 (Howell et al. 2014) observation of Neptune with a cadence of 30 minutes over 50 days (data from Simon et al. 2016, power spectra generated by Apai et al. 2017a).Chavez et al. (2023) compiled infrared observations of Neptune from 1994 through 2022 and determined the cloud activity is variable in an 11-year cycle -with minima in 2002 and maxima in 2015 when the Kepler/K2 data is taken -with strong correlation to solar activity and implication of cloud-top photochemistry by UV radiation from the Sun.
Figure7.Upper : The light curve fit of the 30-50 hours segment using the three-sine model (red curve) corresponding to the data (blue points).Thin vertical lines show the photometric error.Amplitude, period, phase information, and the reduced chi-squared (χ 2 ν = 0.2718) are displayed in the title.Lower : the fit residual (data minus fit).The standard deviation (STD) of the fit residual is displayed in the legend.

FluxFigure 8 .Figure 9 .
Figure8.Upper : The light curve fit of the 183-203 hours segment using the three-sine model (red curve) corresponding to the data (blue points).Thin vertical lines show the photometric error.Amplitude, period, phase information, and the reduced chi-squared (χ 2 ν = 0.2675) are displayed in the title.Lower : the fit residual (data minus fit).The standard deviation (STD) of the fit residual is displayed in the legend.

Figure 10 .
Figure10.Light curve fit for the 410-490 hours segment.The light curve fit uses a four-sine model (red curve) on the binned data with a 20-minute cadence (blue points), based on the 10-minute cadence original flux (thin blue line).Thin vertical lines show the photometric error.Amplitude, period, phase information, the reduced chi-squared (χν = 0.7412), and the Bayesian Information Criterion score(BIC=−1954.26)are displayed in the title (for a discussion of BIC interpretation, see Appendix B).The fit residuals (data minus fit) and the standard deviation (STD) of the fit residuals are also shown.For the result of a three-sine fit for the same data segment, see Figure20.stantamplitude, phase and period).As shown in, the patterns of flux variability are in the shape of sharp peaks that evolve sporadically in a seemingly non-stationary manner.First, we compare the flux values below and above the equilibrium level of 1.By assuming a baseline level of flux and evaluating flux increase and decrease around the baseline, we can identify excess flux that indicates potential atmospheric features, for example, storm spot that produces excess flux on a timescale close to selfrotation.Figure12is the comparison of flux value above and below the equilibrium value of 1, in which flux below 1 is vertically mirrored.Numerically integrating the flux, we find the values 5.4 and 4.8 for the flux below 1 and flux above 1, respectively.Figure12shows that a dip is usually followed by a peak in flux such that the light curve appears relatively even across the baseline level.

Figure 11 .Figure 12 .
Figure11.Resulting light curve from filtering only the periods within 15-95 hours using a Fourier bandpass (bold points).The background data points display the original light curve data.The inset plot shows that the long-period filter captures the mean evolution of the data.

Figure 13 .
Figure 13.Long-period light curve (red) as a result from frequency filtering of periods component from 15-95 hours.Top-tobottom panels show different 100-hour segments and the comparison between the original and long-period light curves.

Figure 14 .
Figure 14.Left: Synthetic light curve composed of 2 sine waves with roughly the same periods (5 hours) with added uniform noise.Right: The periodogram of this light curve.The synthetic light curve does not show the presence of k = 2 waves despite having zonally off-phased waves, suggesting k = 2 period components must come from actual waves in the atmosphere.
Synthetic Jupiter periods maps and the resulting histograms for two case: full-disks including equatorial and mid-latitudes, and only mid-latitudes to polar region.

Figure 15 .
Figure 15.Panel (a) -The synthetic period distribution of Jupiter based on Cassini windspeed data as a function of latitude from Porco et al. (2003).Panel (b) -The superimposed k = 1 and k = 2 periodograms of Luhman 16 AB, where periods of the k = 2 waves are scaled up two times such that the 2-3 hour range matches the 4-6 hour range.The k = 2 waves show a narrower period distribution compared to k = 1.Panel (c) -Left: Synthetic image of 20 • -inclined Jupiter showing effective periods distribution of all latitudes; Right: Periods distribution only from latitudes above 35 • N and below 35 • S; Center : Histograms of Jupiter effective period for the two cases.The period distribution significantly narrows in range once the equatorial windspeeds are removed since these latitudes contain the slowest and fastest winds.

Figure 16 .
Figure 16.The synthetic period distribution of Saturn based on Cassini windspeed data (García-Melendo et al. 2011) as a function of latitude.Left: Zero-inclination, synthetic image showing periods distribution of all latitudes.Right: Zeroinclination synthetic image showing periods distribution only from latitudes above 35 • N and below 35 • S. Center : Histograms of effective period for two cases.

Figure 17 .
Figure 17.Modelled light curve and periodogram power spectrum from a GCM presented in Tan & Showman (2021).The model behind this light curve assumes a rotation period of 5 hours and a moderate basal drag with a drag timescale τ drag = 10 6 s, such that meridionally broad westward and eastward jets with a peak speed of ∼ 500 ms −2 formed in the model.Upper panel : modeled light curve spanning over 1300 hours with a sampling interval of 10 minutes which is the same as that of the TESS dataset in this work.The black line is the original light curve and the red line is the filtered light curve in which signals with periods less than 15 hours are filtered out.Lower left: Periodogram power spectrum of the modeled light curve, and the dashed line highlights the underlying rotation period.Lower right: Periodogram power spectrum zooming into the vicinity of the rotation period (blue line), and the red line is the power spectrum with twice of the period, highlighting the k = 2 component.Dashed lines from left to right represent signals Doppler-shifted by a jet with eastward velocities of 2, 1, 0, -1, and -2 km s −1 .Both the k = 1 and k = 2 components show redshifts compared to the rotation period, and this is because of the Doppler shift by the broad, westward equatorial jet formed in the model.cate cloud deck just aligns/approaches the optical depth probed by observers.By combining the four facts above and following them to their conclusion, we provide the following potential physical interpretation for the planetary-scale waves in Luhman 16B: In this object, we observe the top of the cloud deck close to the observable optical depth.Smallscale turbulence and mixing drive waves of increasing wavelengths, all the way to the longest-possible wavelength standing waves with wavelength equaling the circumference of the object.As the circumference changes with latitude and the rotating gaseous also inhibits differential rotation, the longest wavelength and the rotational period corresponding to the longest waves will vary somewhat with latitude.Waves with vertical components will periodically move gas parcels up and down

Figure 18 .
Figure 18.TESS sources within 2 arc-minutes of Luhman 16 AB with comparable or brighter magnitude, assessed for potential contaminants.TESS Identification Number and TESS-magnitude of each source are included in the legend.