Optimal Frequency-domain Analysis for Spacecraft Time Series: Introducing the Missing-data Multitaper Power Spectrum Estimator

While the Lomb–Scargle periodogram is foundational to astronomy, it has a significant shortcoming: the variance in the estimated power spectrum does not decrease as more data are acquired. Statisticians have a 60 yr history of developing variance-suppressing power spectrum estimators, but most are not used in astronomy because they are formulated for time series with uniform observing cadence and without seasonal or daily gaps. Here we demonstrate how to apply the missing-data multitaper power spectrum estimator to spacecraft data with uniform time intervals between observations but missing data during thruster fires or momentum dumps. The F-test for harmonic components may be applied to multitaper power spectrum estimates to identify statistically significant oscillations that would not rise above a white noise–based false alarm probability. Multitapering improves the dynamic range of the power spectrum estimate and suppresses spectral window artifacts. We show that the multitaper–F-test combination applied to Kepler observations of KIC 6102338 detects differential rotation without requiring iterative sinusoid fitting and subtraction. Significant signals reside at harmonics of both fundamental rotation frequencies and suggest an antisolar rotation profile. Next we use the missing-data multitaper power spectrum estimator to identify the oscillation modes responsible for the complex “scallop-shell” shape of the K2 light curve of EPIC 203354381. We argue that multitaper power spectrum estimators should be used for all time series with regular observing cadence.


INTRODUCTION
Starting with Schuster (1898Schuster ( , 1906)), astronomers investigating periodic phenomena have relied on statistical methods for estimating power spectra. 1 The standard Schuster periodogram is where x n is the observed time series, n is the time index, N is the number of observations, f is frequency, t n is time, and • denotes a statistical estimator.Although there exists no unbiased, nonparametric power spectrum estimator Ŝ(f ), the Schuster periodogram is asymptoti-cally unbiased, meaning it approaches the true physical power spectrum S(f ) as the number of observations approaches infinity (Percival & Walden 2020): (where E{} is expected value).But the number of observations required for the Schuster periodogram to yield an acceptably accurate power spectrum estimate is often too large to be practical: Thomson & Haley (2014) show a barometric pressure time series with N = 1000 for which the Schuster periodogram is biased at high frequencies by seven orders of magnitude.
Ideally our power spectrum estimator would be consistent, such that its variance decreases with increasing sample size.Unfortunately, the variance of the Schuster periodogram follows Var{ Ŝp (f )} = E{| Ŝp (f ) − S(f )|2 } ∝ S 2 (f ) (4) (Bloomfield 2000, Ch. 9).Since the right-hand side of Equation 4 does not depend on N , the Schuster periodogram is inconsistent: its variance does not decrease as more observations are obtained.Any periodogram computed from a finite-length time series is therefore susceptible to variance-induced false positives.
In the mid-to-late twentieth century, statisticians developed consistent, low-bias power spectrum estimators (Bartlett 1948;Blackman & Tukey 1958;Welch 1967;Thomson 1982).But the fact that ground-based observations always have uneven observing cadence has kept such power spectrum estimators out of astronomy, since they are formulated for time series with uniform observing cadence (constant t n − t n−1 ).Instead, astronomers typically use the Lomb-Scargle periodogram developed by Lomb (1976) and Scargle (1982) and generalized by Zechmeister & Kürster (2009).The Lomb-Scargle periodogram has facilitated decades of groundbreaking discoveries, including double degenerate binaries (Bragaglia et al. 1990;Marsh et al. 1995), sunlike magnetic activity cycles in nearby stars (Baliunas et al. 1995), and planets orbiting Proxima Centauri (Anglada-Escudé et al. 2016;Damasso et al. 2020;Faria et al. 2022).But the Lomb-Scargle periodogram shares the inconsistency of the Schuster periodogram and also inherits bias from the uneven observing cadence (Dawson & Fabrycky 2010;VanderPlas 2018).We seek an astronomical power spectrum estimator that is both minimally biased and consistent.
The gold-standard method for power spectral analysis of time series with uniform observing cadence, the multitaper (Thomson 1982;Bronez 1992;Stoica & Sundin 1999;Percival & Walden 2020, Ch. 8), is a nonparametric, Frequentist method that optimally solves the bias and variance problems inherent to periodograms.Multitaper suppresses variance by averaging together multiple independent, low-bias estimators of the power spectrum.Each independent power spectrum estimator S (k) xx (f ) comes from premultiplying x n with a set of weights or taper w (k) n before computing the Fourier transform.With K tapers, the variance of the final power spectrum estimator Ŝxx (f ) is lower than the periodogram variance (Equation 4) by a factor of almost 1/K.For example, using K = 8 and a time-bandwidth product (to be defined in §2.1) of 5, one achieves a fivefold reduction of false positive peaks at all levels of statistical significance (Thomson & Haley 2014).
The classical multitaper method of Thomson (1982) can only be applied to time series with uniform timesteps, a condition that is almost never true in astronomy (though multitaper has been used to investigate the angular power spectrum of cosmic microwave background fluctuations; Fowler et al. 2010).However, Bronez (1988) and Chave (2019) extended the method to time series sampled at regular intervals that are missing some observations.The Bronez-Chave technique is well suited to time series from missions such as Kepler and K2, in which regularly sampled time series are gapped when the reaction wheels dump momentum (García et al. 2011;Vanderspek et al. 2018), the thrusters fire unpredictably (Saunders et al. 2019), or the observations suffer argabrightenings that the data pipeline can't remove (Jenkins et al. 2010;Vanderspek et al. 2018).
Here we explain basic multitaper principles and demonstrate how to apply the missing-data multitaper power spectrum estimator to spacecraft time series. 2We introduce the mathematics of the classical multitaper method in §2, then demonstrate its extensions-tapers for time series with missing data, confidence intervals, and the F-test for identifying statistically significant oscillations-in §3.In §4, we showcase both the accuracy and the precision of the missing-data multitaper method by using it to diagnose antisolar rotation in Kepler observations of the young sunlike star KIC 6102338.In §5, we identify the set of modes that produce one of the complex "scallop-shell" light curves discovered by Stauffer et al. (2017) and Stauffer et al. (2018).We present our conclusions in §6.A software package that performs classical and missing-data multitaper calculations, Multitaper.jl, is discussed on page 23.

MULTITAPER BASICS: THE CLASSICAL MULTITAPER METHOD
The multitaper power spectrum estimator mitigates what Percival (1994) described as the "tragedy of the periodogram". 3As discussed in §1, the tragedy is twofold: periodograms are both biased and inconsistent.Here we introduce basic principles of multitaper power spectrum estimation.We begin with the physical motivation for tapering ( §2.1), then proceed to taper calculation ( §2.2).In §2.3 we show a sample set of tapers, then in §2.4 we describe how to estimate the power spectrum.We discuss bias and accuracy in §2.5 and provide guidance for bandwidth choice in §2.7.We finish the section with an example multitaper power spectrum estimate in §2.8.

What is a taper?
The classical multitaper method introduced by Thomson (1982) employs the discrete prolate spheroidal sequences (DPSS) discovered by Slepian (1978) as tapers w (k) n for n = 0, . . ., N − 1 and k = 0, . . ., K − 1. Tapers are sequences that produce optimally shaped spectral windows for Fourier analysis.The expected value of the k th tapered power spectrum estimator is where * denotes convolution and is the k th spectral window.The goal of tapering is to minimize biasing leakage, that is, spurious power at a given frequency that is actually attributable to power from a far-away frequency.The ideal spectral window is W (f ) = δ(f ), for which the expected value of the power spectrum estimator would equal the true power spectrum.However, there is no finite-length taper w n that Fourier transforms to δ(f )-it's only possible to achieve W (f ) = δ(f ) for time series and tapers of infinite length.
The boxcar function w n = √ N , which is the "taper" that results from having a finite-length dataset that is otherwise untapered, has a sinc 2 function as a spectral window, which leaks power into its high sidelobes.For an introduction to tapering principles, see Harris (1978).See Percival & Walden (2020, Chapter 8) and Mudelsee (2014, Chapter 8) for textbook descriptions of the DPSS tapers.
Since the DPSS are orthogonal, such that w xx (f ) therefore produces a multitaper power spectrum estimator with variance reduced by a factor of ≃ 1/K from any single-taper power spectrum estimator. 4We can quantify the quality of our power spectrum estimator using the mean-squared error (MSE): (Percival & Walden 2020, Ch. 8).Multitaper has the lowest MSE of any nonparametric power spectrum estimator, including the Schuster periodogram, the Lomb-Scargle periodogram, Bartlett's method (Bartlett 1948), and Welch's method (Welch 1967;Dodson-Robinson et al. 2022).

Calculating the multitapers
The DPSS tapers w (k) n are defined by the eigenvalue equation To solve Equation 8, all observation timestamps must be divided by the observing cadence ∆t obs in order to construct a time series with integer timestamps and unit sampling, ∆t = t n+1 − t n = 1.When there are no missing values in the time series, the matrix sin(2πϖ ) is Toeplitz and commutes with a tridiagonal matrix, enabling the use of computationally efficient eigenvalue routines to solve for the eigenvectors w n (Grünbaum 1981).In Equation 8, ϖ is the bandwidth, or half-width of the main lobe of the spectral window.5 Oscillations in the time series will manifest as peaks in the power spectrum estimate that have full width 2ϖ.We will discuss bandwidth selection in more detail in §2.7.Throughout this paper we will use ϖ to denote bandwidth in rescaled units where ∆t = 1 and π = ϖ/∆t obs to refer to bandwidth in physical units.
The DPSS tapers have the desirable property of being maximally bandlimited, which means they provide the highest spectral concentrations inside −ϖ < f < ϖ that can be obtained for any finite-length, discrete sequence.In units where ∆t = 1, spectral concentrations are defined as where 1/2 cycle per time unit is the maximum frequency probed by the time series (Nyquist frequency) and the λ (k) in Equation 9are equal to the eigenvalues in Equation 8. Equation 9 measures the fraction of power in W (k) (f ) that is confined to the frequency interval −ϖ < f < ϖ, so that 0 < λ (k) < 1.The higher the value of λ (k) , the lower the bias contribution to the mean-squared error (Equation 7; see also §2.5).

A sample set of tapers
Figure 1 shows the Kepler Q0 long-cadence observations of the differential rotator KIC 6102338.Out of N = 1639 flux measurements in the PDCSAP6 light curve downloaded from the MAST archive,7 only 16 had been rejected by the Kepler pipeline; these were denoted NaN in the light curve table.To construct a time series with uniform observing cadence and no missing data on which to demonstrate the classical multitaper technique, we used a simple interpolation technique where we replaced each NaN with the immediately preceding PDCSAP flux.
Figure 2 depicts Equation 8 applied to the timestamps from the Kepler Q0 KIC 6102338 observations.The top panel of Figure 2 shows the six DPSS's w (k) t , k = 0, . . ., 5 that yield the highest matrix eigenvalues λ (k) .The bottom panel of Figure 2 shows the positive halves of the K spectral windows W (k) (f ) (spectral windows are symmetric about f = 0).Spectral concentrations are shown in the legend.The characteristics of the spectral window are very important, as we will see in §3.2.Note that the frequency axis in Figure 2 is logarithmic, so the main lobes of the spectral windows are quite narrow (ϖ = 0.002441, which covers only 0.48% of the frequency domain).

Estimating the power spectrum
Recall that the orthogonal DPSS's of Equation 8 provide K independent power spectrum estimators The multitaper power spectrum estimate is a weighted average of the eigenspectra where the averaging limits the variance of the estimator.At frequencies where the power spectrum is varying rapidly (i.e.|d Ŝxx (f )/df | is high), the weighting adjusts to favor minimizing bias over variance by downweighting eigenspectra with lower spectral concentrations, as leakage could otherwise render the estimated power spectrum slope artificially shallow.We will demonstrate how to calculate the bias of Ŝxx (f ) in §2.5 and estimate the variance of Ŝxx (f ) in §3.4.The frequency-dependent weights, d(f ), are chosen so that each Ŝ(k) xx (f ) satisfies the equation = 0 (12) (Thomson 1982, §V), where s 2 x is the sample variance of the time series: The coupled equations 11 and 12 are solved iteratively so as to produce estimates Ŝxx (f ) that match within a user-defined tolerance.All functionality for calculating Ŝxx (f ) according to Equations 11 and 12 is included in the Multitaper.jlsoftware package (Haley & Geoga 2020a) used for this paper, as well as in the R multitaper package and the python package by Prieto (2022).Multitaper outperforms the Welch (1967) power spectrum estimator adapted for astronomical use by Dodson-Robinson et al. (2022).If one holds two of bandwidth, bias, and variance at fixed values, the third quantity is smaller for the multitaper estimator than for any other power spectrum estimator (Bronez 1992).

Bias and accuracy
Multitapering optimally minimizes the broadband bias defined in Equation 3.An upper limit to the broadband bias is bias Ŝxx (Thomson 1982), where λ (K−1) = min{λ (k) } is the spectral concentration of the highest-order DPSS and min{•} denotes minimum.In Figure 2, our choice of K = 6 limits the bias to 0.0075 s 2 x (λ (5) = 0.9925).Equation 14 guarantees that the multitaper power spectrum estimate cannot be excessively far from the true power spectrum; no such guarantee exists for the Schuster and Lomb-Scargle periodograms.In practical application, multitaper estimators are likely to have orders of magnitude lower MSE (Equation 7) than any other nonparametric estimator when the dynamic range of the true power spectrum is large (see example in Thomson & Haley 2014, Figure 1).
Although astronomers often focus on estimating oscillation frequencies, there are applications for which accurately measuring the power in either the oscillations or the underlying broadband continuum is critical.For example, Alfvén wave excitation in the solar photosphere is traced by the motion of magnetic "bright points," which anchor flux tubes that reach up into the corona.From radiative magnetohydrodynamic simulations of granulation in the solar photosphere, Van Kooten & Cranmer (2017) predicted a bright-point motion power spectrum of S(f ) ∝ f −1 , in agreement with observations from the Swedish 1m Solar Telescope and the Solar Optical Telescope (Chitta et al. 2012).Without an accurate E{S(f )} from observations, it would be much more difficult to assess the validity of the Van Kooten & Cranmer (2017) model at describing both the coronal heating rate and the granule size spectrum.Asteroseismic extraction of stellar parameters also relies on accurately representing both the total power and the shape of the Gaussian oscillation envelope (e.g.Kallinger et al. 2010).In §4 we show an example where accurately estimating the power in each detected oscillation is crucial to the science goal.

Why multitaper does not use error bars
It may seem unusual that Equations 10-12 do not include uncertainties.Unfortunately, weighting the data based on nonuniform error bars, as in the generalized Lomb-Scargle periodogram and its Bayesian extension (Zechmeister & Kürster 2009;Mortier et al. 2015), produces a "taper" for which the spectral window is shaped by the error bars (Equation 6).The expected value of the power spectrum estimator then departs from true power spectrum in unquantifiable ways (Equation 5; see also the discussion of accuracy in §2.5).Furthermore, the Bayesian and generalized Lomb-Scargle periodograms become dimensionless when computed with uncertainty-derived weights; instead of power density, the quantity estimated is either an approximate signalto-noise ratio or a probability.Since the data are scaled by uncertainty-derived weights before the computation of trigonometric sums (Zechmeister & Kürster 2009;Mortier et al. 2015), there is no straightforward conversion back into physical units.
If the observer's goal is to detect a small number of statistically significant oscillations, it is less critical that the power spectrum estimate be minimally biased across all frequencies.For example, radial velocity planet hunters will often model a time series as a sum of 1-5 Keplerian orbits plus a Gaussian process to treat stellar activity.For such applications, it is sufficient that the periodogram peak corresponding to each oscillation sticks out sufficiently far above the continuum; accurate, dimensional estimates of each signal's power are not required.(Note, however, that spectral leakage from strong oscillations can hide weaker oscillations, so some control over bias is important even for signal detection.)But as we learned in §2.5, studying phenomena with broadband power spectra such as granulation and turbulence requires quantifying the contribution of sinusoids of all frequencies to the time series (not just the strongest oscillations).
Furthermore, we note that the error bars on most modern spacecraft photometry are both miniscule when compared to the quantity being measured and nearly uniform.The left panel of Figure 3 shows a histogram of 1σ errors on the Kepler Q0 observations of KIC 6102338 (Figure 4), while the right panel contains a histogram of fractional errors from the same dataset.The distribution of absolute flux error spans < 0.1 e − s −1 , while the fractional flux error does not exceed 4.65 × 10 −4 , in keeping with Kepler's mandate to deliver photometry with parts-per-10,000 precision.Twinkle, a forthcoming Earth-orbiting satellite that will deliver high-cadence, low-resolution spectroscopy spanning 0.5-4.5µm,will have ∼ 1 part-per-thousand precision (Stotesbury et al. 2022).It is acceptable to leave out error bars when estimating the power spectra of such high-quality data, a practice followed by the Kepler Asteroseismology Program (e.g.Kallinger et al. 2010).Finally, when the error bars are uniform, the weighted generalized Lomb

Choice of bandwidth
The choice of N ϖ allows the observer to trade bias for variance.Equation 8 yields 2N ϖ maximally bandlimited tapers w (k) n with λ ≃ 1 (Thomson 1982).Highorder DPSS with k > 2N ϖ have λ ≪ 1, so do not make useful tapers; the maximum allowable value of K is therefore 2N ϖ.To improve variance suppression, the observer can choose a larger bandwidth so as to increase K. To prioritize bias suppression, which is important if the physical process under study produces a power spectrum with high dynamic range, the observer should use a narrower bandwidth and fewer tapers.The control of the bias-variance tradeoff through the adjustment of bandwidth is widely recognized as the main performance advantage of the multitaper method (Bronez 1992;Haley & Anitescu 2017).For more details about how ad- justment of the bandwidth can be used to trade bias for variance, see Percival & Walden (2020).
To construct a multitaper power spectrum estimate, the observer sets ϖ such that the time-bandwidth product N ϖ is a small integer or half-integer, which creates a full bandwidth 2ϖ that is an integer number of Rayleigh resolution units.As in optics, the Rayleigh resolution R gives the minimum frequency separation required for two oscillations to be statistically distinguishable: is λ (7) = 0.6988, which means more than 30% of the power in an oscillation with frequency f 0 leaks out of the interval f 0 − ϖ < f 0 < f 0 + ϖ.Since λ (6) = 0.9367, which allows more than 6% spectral leakage outside of the bandwidth, an especially bias-averse observer might set K = 6, k = 0 . . .5, as we have done.A geophysics example that prioritizes low bias is the measurement of the elastic thickness of the central Australian lithosphere, for which Swain & Kirby (2003) To maintain the consistency of the multitaper estimator, it is important to increase K (and accordingly N ϖ) along with the sample size.Haley & Anitescu (2017) present a method for finding the bandwidth that minimizes the MSE using the second derivative of a spline fit to Ŝxx (f ) on a frequency grid with spacing 2ϖ.Heuristically, their experiments with a Lorentzian power spectrum yield optimal time-bandwidth products N ϖ = 5 for N = 500, N ϖ = 8 for N = 1000, and N ϖ = 11 for N = 1500.However, the best choice of N ϖ depends on the process power spectrum and the science goal.Since all applications in this manuscript involve detecting periodic signals, we use smaller timebandwidth products that allow more precise oscillation frequency estimation.Seismologists often use narrow bandwidths (small N ϖ) in order to detect oscillations: Chave et al. (2019) selected N ϖ = 8, K = 15 to identify closely spaced seismic modes in a time series of N = 597, 504 deep-ocean pressure measurements, while PéRez-Gussinyé et al. (2007) chose N ϖ = 3, K = 5 in their calculation of lithospheric elastic thickness from seismic tomography.For astronomical time series with 1000 ⪅ N ⪅ 5000, we recommend values of N ϖ between 4 and 8.For N ⪅ 1000, N ϖ =3-4 is a good starting point.Thomson (1982) recommends that the observer use higher N ϖ to characterize the continuum and lower N ϖ to study the oscillations, so that a detailed multitaper analysis of a single time series may involve several power spectrum estimates with different N ϖ.

A classical multitaper power spectrum estimate
The top panel of Figure 4 shows the multitaper power spectrum estimate Ŝxx (f ) based on the K = 6 tapers shown in Figure 2.For all time series shown in this paper, we subtract the sample mean x = (1/N ) N −1 t=0 x n -which is a zero-frequency term that can cause bias-before estimating the power spectrum.The light blue shading in Figure 4 shows the 95% confidence interval, the calculation of which we will discuss in §3.4.The vertical bar of the blue cross indicates the expected value of the jackknife variance estimate (see §3.4), while the horizontal cross bar indicates the bandwidth π. Green dots identify signals that exceed a significance threshold of 1 − 1/N = 0.99939 according to the F-test for harmonic components ( §3.5).The middle panel shows the p-value from the F-test conducted at each frequency; green dots indicate p-values below 1/N = 0.00061, for which we reject the null hypothesis that the dataset contains no oscillation component at frequency f .Figure 4 demonstrates the variance reduction provided by multitapering: the 95% confidence interval on Ŝxx (f ) does not come close to encompassing the variability in Ŝp (f ).For a time series x n consisting only of white noise, Ŝxx (f ) forms a χ 2 2α -distributed stochastic process with α = d 2 k (f ) ⪅ K, which has stronger correlations across frequency and thus fewer large excursions than does the χ 2 2 distributed periodogram Ŝp (f ).When equal weighting at all frequencies is used, α = K.The distribution of the power spectrum estimator is derived in the Brillinger (1975) textbook.

MULTITAPER EXTENSIONS: MISSING DATA, CONFIDENCE INTERVALS, AND THE F-TEST
Here we show how to compute multitaper power spectrum estimates from spacecraft data and identify statistically significant signals.In §3.1 we demonstrate tapering for time series with gaps using the Chave (2019) method.We compare the spectral windows produced by the Chave multitaper method and the generalized Lomb-Scargle periodogram in §3.2.Finally, we demonstrate two methods for assessing the statistical significance of oscillatory signals: jackknife confidence intervals ( §3.4) and the F-test for harmonic components ( §3.5).

Missing-data multitaper estimator
Almost all time-series observers have to contend with missing data due to instrument downtime or poor environmental conditions.In space-based astronomy, pieces of dust may reflect sunlight into the telescope barrel, yielding sudden background flux increases that cannot always be modeled out of the light curves.Furthermore, Kepler and K2 observations were interrupted once per month during data downlink (see García et al. (2011) for a description of all types of Kepler observing interruptions).The TESS reaction-wheel momentum dumps that occur every 2.5 days increase the pointing jitter to ∼ 1' for a few minutes at a time, creating low-quality light curve segments that some observers clip out (e.g.Plachy et al. 2021).
Chave (2019) extended the notion of maximally bandlimited orthogonal sequences introduced by Slepian (1978) to those sampled on a grid where there is missing data.Such sequences are defined only at timestamps for which data are available.As in §2, we set the bandwidth ϖ for the desired spectral window, and let the sequences be defined on the length N time grid {t n } N −1 n=0 with unit sampling (∆t = 1) except where there are gaps-for example, N = 7 and {t n } = {0, 1, 2, 3, 7, 8, 9}.The missing-data Slepian sequences (MDSS) solve an eigenvalue problem similar to Equation 8: except on the missing-data time grid {t n } N −1 n=0 instead of the consecutive integers 0, . . ., N − 1.The MDSS form an orthonormal set on {t n } N −1 n=0 and can be used as tapers to produce a missing-data multitaper (MDMT) power spectrum estimate: where F denotes the nonuniform fast Fourier transform (NFFT).8Although it is possible to use frequencydependent weighting (Equations 11 and 12) for missingdata multitaper, it often yields power spectrum estimates with too few degrees of freedom for the F-test ( §3.5), so we do not use it here.We note that Springford et al. (2020) presented a multitaper method for Kepler data that relies on interpolating the DPSS tapers.However, interpolated DPSS are in general not orthogonal on the missing-data time grid, so the independence of the Ŝ(k) xx is lost and jackknife confidence intervals cannot be computed.
Figure 5 shows the PDCSAP fluxes from the Kepler Q16 observations of KIC 6102338 (top) and a set of MDSS tapers w (k) n made for those observations (bottom), for which N = 3534, N ϖ = 4 ( π = 0.05539 day −1 ), and K = 6.Gaps are highlighted by gray shading.Note that all w (k) n approach zero at both edges of the large gap near the beginning of the time series: this behavior suppresses spectral leakage by removing the false high-frequency signal created by the sudden cutoffs in x n at the gap edges.(In fact, given the near-zero weighting for the earliest data, it would be reasonable to limit the analysis to BJD> 2456322.)Regardless of N and N ϖ, the two DPSS/MDSS tapers with the highest spectral concentrations, w  edges of major gaps.We also see w n , w n , w n , and w (4) n approaching zero at the edges of the smaller gap that begins at BJD 2456357.96.
The top panel of Figure 7 shows the missing-data multitaper power spectrum estimate Ŝxx (f ) from the KIC 6102338 Q16 data and tapers shown in  The k th MDSS spectral window W (k) (f ) can be calculated using Equation 6, replacing the standard fast Fourier transform with the nonuniform fast Fourier transform (Barnett et al. 2019;Barnett 2020).Figure 6 shows the positive half of each W (k) (f ) calculated from the w (k) n pictured in Figure 5 (as in the DPSS, the MDSS spectral windows are symmetric about zero).All W (k) (f ) are orthonormal over the frequency range −1/2 < f < 1/2 and orthogonal on −ϖ < f < ϖ.Importantly, the spectral window of the missing-data multitaper estimator can be written analytically in terms of the individual W (k) (f )'s.If one replaces the adaptive weights d 2 (f ) (Equation 12) with 1/K, one obtains the following expression for the missing-data spectral window: The MDSS spectral window delivers a lower dynamic range than the DPSS spectral window, as we can see by comparing the two versions of W (0) (f ) in Figures 2 and  7: the DPSS W (0) (f ) drops by over 10 orders of magnitude over the frequency range 0 < f < ϖ, while the MDSS W (0) (f ) drops by less than five orders of magnitude over 0 < f < ϖ.In general, missing-data tapers have lower spectral concentrations than classical tapers, leading to more spectral leakage.If the dataset is only missing a few scattered data points, as in the KIC 6102338 Q0 data, a bias-averse observer may wish to fill in the missing data and use the classical multitaper method instead of the Chave (2019) missing-data method.However, despite their different bias reduction capabilities, the missing-data and standard multitaper power spectrum estimators share the formal property that all power spectrum peaks have the same rectangular shape as the spectral window, which is in general not true for the Lomb-Scargle periodogram (see Scargle (1982, Appendix D) and §3.2 on pseudowindowing).

Spectral window comparison: Lomb-Scargle vs. missing-data multitaper
From Equation 5we know that the spectral window contributes bias by convolution with the true power spectrum, the characteristics of which require some special consideration when there are missing data.Figure 3 of Scargle (1982) shows dramatic comparisons between spectral windows from time series with unequal observing cadence vs. unit observing cadence.Scargle (1982) recommends the following heuristic to approximate the spectral window, or pseudowindow, at each frequency: compute the Lomb-Scargle periodogram of a sinusoid with frequency f (where f is far from zero) on the set of observation times {t n } N −1 n=0 .Plotting this pseudowindow, which is convolved with the true process power spectrum in the neighborhood of frequency f , gives an idea of the spectral leakage.
As an example of pseudowindowing, we take time stamps from CHEOPS observations of 55 Cnc obtained on 22-23 February 2021 , PI: B. Demory).We construct two synthetic signals of the form x n = sin(2πf 0 t) with f 0 = 25 day −1 and f 0 = 80 day −1 on the grid of observation times.The top panel of Figure 8 shows the synthetic signal with f 0 = 25 day −1 .From each signal, we compute a Lomb-Scargle periodogram and MDMT power spectrum estimate.The bottom panel shows the MDSS computed for N ϖ = 5, K = 7, π = 7.61 day −1 , with the legend showing the spectral concentrations λ (k) .We selected N ϖ = 5, rather than N ϖ = 4 as for the Kepler Q16 data, in order to widen the bandwidth and increase the spectral concentrations: with N ϖ = 4, the CHEOPS timestamps give λ (0) = 0.643 for W (0) (f ), compared with λ (0) = 0.769 for N ϖ = 5.For a given N ϖ, time series with more gaps (e.g.CHEOPS) will have lower MDSS spectral concentrations than time series with fewer gaps (e.g.Kepler Q16).
The top panel of Figure 9 compares the spectral window of the MDMT, W (f ) (purple), with a proxy for the zero-frequency spectral window of the Lomb-Scargle periodogram W LS (f ) (green dashed).We calculate W LS (f ) using the non-uniform fast Fourier transform: where the indicator function is I N (t) = 1 for t = t 0 , ..., t N −1 and zero otherwise (VanderPlas 2018).Our definition of W LS (f ) shows the shape of the zerofrequency Lomb-Scargle spectral window but uses the power spectral density normalization of the multitaper.One immediately obvious feature of W (f ) from the missing-data multitapers is the sidelobe suppression: while leakage in the Lomb-Scargle spectral window manifests as a sawtooth pattern in which the highest sidelobes at f = ±15 day −1 have 14% of the amplitude of the main lobe (the central peak of the spectral window), the MDMT lacks the ±15, 30, 45 day −1 sidelobes.The highest secondary peak in the MDMT spectral window at f = ±57 day −1 has only 1.4% the amplitude of the main lobe, greatly decreasing the possibility of misinterpreting an observing cadence artifact as a true signal (e.g.Koopmans 1995, Ch. 3).
The bottom two panels of Figure 9 show the pseudowindows centered at frequencies f 0 = 25 day −1 (blue solid line) and f 0 = 80 day −1 (orange dashed line) in the Lomb-Scargle periodogram (left) and MDMT (right).In both plots, the x-axis variable is ∆f = f − f 0 and the pseudowindows are zoomed in to |∆f | ≤ 20 day −1 in order to focus on the main lobe and the first pair of Lomb-Scargle sidelobes.While the two Lomb-Scargle pseudowindows have similar shapes to W LS (f ) and to each other, they are not identical: at ∆f = 8 day −1 they differ by more than a factor of three, while at ∆f = −20 day −1 they differ by an order of magnitude.In contrast, the MDMT behaves almost as well in reality as it does in theory: the most significant difference between the f 0 = (25, 80) day −1 pseudowindows is 20%.The CHEOPS time series shows the value of reshaping the spectral window with the MDSS, even when the gaps are long and numerous.

Spurious peaks
In this section we briefly summarize a theoretical result concerning significant peaks in nonparametric power spectrum estimators (Thomson & Haley 2014).Restricting ourselves to the unit sampling (∆t = 1) case for the moment, we know from §2.8 that periodograms have a χ 2 2 or exponential distribution and that multitaper power spectra have a χ 2 2K distribution (when frequency-dependent weighting is used, the distribution is χ 2 α where α = 2 K k=1 d 2 k (f ) at frequency f ).One can then assign significance to large peaks in a power spectrum estimate (after one has subtracted the baseline broadband spectrum, or prewhitened ) using the appropriate χ 2 quantile.The contribution in Thomson & Haley (2014) was to show that the tendency of the power spectrum estimate to spuriously upcross a high threshold, z, depends on the second derivative of the serial correlation (dependency across frequency) of the estimator.These upcrossings occur at a rate of for the periodogram and for the multitaper estimator (upcrossing rates are expressed per normalized frequency unit, or Rayleigh resolution).
Equations 20 and 21 confirm that multitaper produces fewer spurious upcrossings per significance threshold by a factor of almost 2N ϖ.This statistical fact is further exemplified in Thomson and Haley (2014) with numerical experiments on synthetic white-noise time series and proton density solar-wind data from the Advanced Composition Explorer spacecraft (their Table 2).Thomson & Haley (2014), Figure 2a shows that when the time-bandwidth product N ϖ is chosen as 5 (and, in fact this holds for all N ϖ > 1), the multitaper power spectrum gives approximately five times fewer spurious peaks than the periodogram (with white noise, all peaks are spurious.)As an aside, we might just find significant telescope time savings when faced with a, say, fivefold decrease in spurious oscillation detections as a consequence of switching from periodograms to the multitaper method!With multitaper, we can go further than just estimating a power spectrum-we can also compute confidence intervals.In particular, we are interested in a confidence interval for the logarithm of the power spectrum.Since the log Ŝ(k) xx (f ) are approximately independent, they can be used to estimate variance.The delete one jackknife estimator of Efron & Stein (1981); Thomson (1990b); Thomson & Chave (1991a); Thomson (1994) is less biased than the sample variance of the individual eigenspectra (Cressie 1981) and is closer to the true value Var{log Ŝxx (f )} even for small K.

Calculating confidence intervals
To jackknife, one computes K delete-one estimators of the log-eigenspectrum log Ŝxx (f ).The logarithmic transformation ensures that the set of jackknife pseudovalues, which we will define shortly, is approximately normally distributed at each frequency (Miller 1974, §5.3).For the simple average multitaper estimator, which we use for time series with missing data (Equation 17), the delete-one log-eigenspectra are log Ŝ\m The jackknifed variance estimate is based on the sample variance of quantities termed pseudovalues, and not the delete-one log-eigenspectra themselves.The pseudovalues are defined as and can be treated as if they come from a normally distributed sample (Cressie 1981), which means their variance is Student t−distributed with K − 1 degrees of freedom.Computing the variance of the p m (f )'s gives the following log power spectrum variance estimate: We have constructed the 100(1 − α)% confidence intervals according to as in Thomson & Chave (1991b, Eq. (2.49)), where t K−1 is the Student-t distribution with K − 1 degrees of freedom and α is the significance level.
Figures 4 and 7 both use light blue shading to show jackknife 95% confidence intervals on the multitaper power spectrum estimators.The confidence intervals provide guidance for interpreting the results of the Ftest for harmonic components ( §3.5): a power spectrum peak may be quite high, but if the confidence interval in and around that peak is wide, the F-test p-value might not be statistically significant.The combination of confidence intervals and the F-test allows the observer to take into account the quality of information at a particular frequency before deciding if an oscillation at that frequency is significant.Confidence intervals are also important when the science goal requires not only identifying oscillations, but precisely estimating their power (see example in §4).

F-test for harmonic components
Adjacent to the topic of power spectral estimation, but separate from it, is harmonic analysis, at which the multitaper technique also excels.Harmonic analysis, which is the detection of oscillations, is the typical astronomical application of periodograms.For example, radial velocity exoplanet hunters will begin often begin an analysis with a generalized Lomb-Scargle periodogram, then fit a time-domain model (for example, a Keplerian for a planet or a Gaussian process for rotation) using the period associated with the highest periodogram peak to initialize the fit.The F-test is the harmonic analysis tool associated with the multitaper technique.We refer the reader interested in technical details to the original source paper (Thomson 1982) and give a basic presentation here, with the relevant extensions for time series with missing data.For an accessible textbook overview of the F-test, see Mudelsee (2014, Chapter 5).
With the Lomb-Scargle periodogram, the usual procedure is to assume the noise is white and use either bootstrapping (Cumming 2004) or extreme-value theory (Baluev 2008) to calculate false alarm probabilities (FAPs) that are constant as a function of frequency.However, both turbulence and signal persistence in a bolometer produce red noise, such that the continuum of Ŝxx (f ) underlying any oscillations has more power at low frequencies than high frequencies.Although Delisle et al. (2020) provide a generalization of the Lomb-Scargle periodogram for noise with arbitrary covariance structure, their method requires the user to estimate a correlation time in order to approximate the covariance matrix.Given that the true correlation time is always unknown, Delisle et al. (2020) recommend (1) examining multiple periodograms computed using different correlation times and (2) using a likelihood maximization to calculate a more realistic covariance matrix.It would be useful to be able to judge the significance of a power spectrum peak without needing to make any assumptions about the noise.In addition, a FAP estimate that is robust against large dynamic ranges in Ŝxx (f ) would reduce the need to iteratively subtract sinusoids or Keplerian models from the time series and examine periodograms of the residuals (e.g.Reinhold et al. 2013;Bourrier et al. 2018), a practice that can be error-prone given uncertainties on the correct sinusoid frequencies.Thomson (1982) gives an effective test for the presence of oscillatory components (i.e.purely sinusoidal signals) in a multitaper power spectrum.Thomson's F-test generalizes easily to the missing-data case (Chave 2019).The test statistic is where x(k) (f ), the eigencoefficients, are given by the complex-valued Fourier transform of the tapered time series x n e −2πif t (unit sampling), ( 27) and wk (0) is the Fourier transform of taper k, evaluated at zero frequency: The estimated complex oscillation component, before squaring, is The test statistic (Equation 26) has an F −distribution where the numerator and denominator have 2 and 2K−2 degrees of freedom, respectively.At frequencies where the F-test p-value is small, there is support for the hypothesis that the time series records an oscillation. 9We choose a stringent F-test p-value cutoff of 1/N (Thomson 1990a), for which the expected number of random sampling-based fluctuations that reach our detection threshold is one per power spectrum estimate.The test statistic is to be interpreted as the estimated power at frequency f compared with the power one obtains when a sinusoid with frequency f is removed from the time series.Some detected oscillatory components may not coincide with large peaks in the power spectrum, especially when signal to noise ratios are low.While power spectrum peaks spaced closer than twice the bandwidth (∆f < 2N ϖ) may be difficult to distinguish if classical multitaper and MDMT estimators are the only tools available, careful inspection of the Ftest statistic reveals that its fundamental frequency unit is the Rayleigh resolution (Equation 15).Thus the Ftest can distinguish signals with small enough frequency separation that their rectangular power spectrum peaks overlap in a plot of Ŝxx (f ).It can also pick out longperiod signals with f ≤ N ϖ that are not visually distinct from zero frequency as long as those signals have f > 2R, meaning they are formally resolvable according to the Rayleigh criterion.The multitaper-F-test combination therefore has the same frequency resolution as Bayesian, least-squares, and generalized Lomb-Scargle periodograms, and is suitable for detecting multiple periodicities (Vautard et al. 1992;Robertson & Mechoso 1998;Zhou et al. 2007, see §5, in which we identify 10 distinct oscillation modes in a "scallop shell" light curve).Furthermore, the F-test is valid no matter the type of noise in the dataset, as demonstrated analytically by Thomson (1982Thomson ( , 1990a) ) and computationally for synthetic and observed seismic datasets by Park et al. (1987), Lees (1995), andPrieto et al. (2007), among others.For accurate estimation of oscillation frequencies, oversampling of the frequency grid (or zero padding in the time domain) is recommended.Note, however, that 9 Formally, the hypothesis being tested at each frequency gridpoint f j is that if we use the Munk & Hasselman (1964) representation for the deterministic, first-moment part of the power spectrum, , and C j is the complex Fourier coefficient at frequency gridpoint j), the value of C j is nonzero (Thomson 1990a).The power spectrum is also assumed to have a non-deterministic noise component given by the second moment of dX.Since S(f )df = E{|dX(f )| 2 }, each nonzero C j corresponds to a delta function in the true power spectrum, which manifests as a narrow, rectangular peak in the multitaper power spectrum estimate due to Equation 5. Physically, delta functions in the power spectrum correspond to sinusoids.
oversampling does not add information to the power spectrum estimate; it merely interpolates the estimate.
Figures 4 and 7 both denote significant oscillations identified by the F-test with green dots.In Sections 4 and 5, we will follow the same convention and discuss the information provided by the F-test about each example dataset.According to the study of stellar surface shear by Reinhold et al. (2013), KIC 6102338 exhibits differential rotation with primary period P 0 = 5.27 days and secondary period P 1 = 4.07 days.While angular rotation speed varies continuously as a function of latitude, it is common for two dominant periods to emerge in the photometry, as spots tend to favor certain latitudes (e.g.Frasca et al. 2011;Fröhlich et al. 2012).Reinhold et al. (2013) searched for periodicities using iterative prewhitening: from a generalized Lomb-Scargle periodogram of each PDC-MAP light curve,10 they identified the period of maximum power P 0 , subtracted a sinusoid with period P 0 from the light curve, and then computed a generalized Lomb-Scargle periodogram of the residuals.The procedure was repeated four times, and then P 0 plus the four periods of maximum power from each residual periodogram were used as initial guesses for the periods in a truncated Fourier series fit to the light curve: where F is the flux and the P i were allowed to vary.A detection of differential rotation was recorded when any of P 1 . . .P 4 fell within 30% of P 0 .Here we will recover the same results as Reinhold et al. ( 2013), but without iterative sinusoid subtraction.We will also use the power at the rotation harmonics to measure the direction of the stellar surface shear.

Data and tapers
Reinhold et al. ( 2013) selected the Kepler Q3 data for their study because Q3 showed fewer instrument artifacts than Q0-Q2.To perform a controlled comparison between multitaper and their iterative prewhitening procedure, we also base our analysis on the KIC 6102338 Q3 long-cadence PDC-MAP light curve.We remove > 3σ outliers, which may be contaminated by stellar flares.The top panel of Figure 10 shows the resulting light curve x n (N = 4131).The bottom panel of Figure 10 shows the MDSS tapers with N ϖ = 4.5 and K = 7.We chose a wider bandwidth than we used for Q16 (N ϖ = 4, K = 6; §3.1) in order to improve variance suppression: with N ϖ = 4.5, we can safely include w (6) n , which has λ (6) = 0.894-similar to λ (5) for Q16 (Figure 5).By decreasing the variance of our power spectrum estimate, we reduce the uncertainty on the shape of each power spectrum peak, making it more likely that true oscillatory signals reach high statistical significance in the F-test.High dynamic range in W (f ) allows us to simultaneously detect the strong signals at the fundamental rotation frequencies and their weaker harmonics.

Power spectrum estimate
Figure 11 shows the power spectrum estimates from the KIC 6102338 Kepler Q3 data.In the top panel, the blue line shows the MDMT with the 95% confidence interval shaded in light blue.The red stars indicate oscillatory signals at the 1 − 1/N significance level, according to the F-test.The frequency of the strongest signal, f 0 = 0.193 day −1 , matches the primary rotation signal identified by Boisse et al. (2012); Reinhold et al. (2013) to within 0.3R, which is well below the fundamental resolution unit.We also detect an oscillation at f 1 = 0.249 day −1 , which is again within 0.3R of the secondary rotation signal reported by Reinhold et al.The multitaper-F-test combination reproduces their differential rotation detection, but does not require any iterative optimization.
The detections at f 0 and f 1 highlight the fact that the F-test can distinguish signals that are not separated by a full bandwidth: visually, the two rectangluar power spectrum peaks at f 0 and f 1 overlap because they trace two oscillations that are separated by ∆f < 2N ϖ (in this case, 2 π = 0.107 day −1 , but f 1 −f 0 = 0.0553 day −1 ).However, both corners of the f 0 peak and the right-hand corner of the f 1 peak are clearly visible.When examining multitaper power spectrum estimates, it's helpful to remember that power spectrum peaks are approximately rectangular; sharp corners usually belong to oscillatory signals.
Figure 11 shows more than just the two fundamental rotation frequencies.We also find significant oscillations at 2f 0 and 2f 1 , the first harmonics of each rotation frequency.Furthermore, if we lower our statistical significance threshold from 1 − 1/N to 0.999, we pick up an oscillation at 3f 0 , the second harmonic of the primary rotation signal (green dot).While detections of differential rotation without harmonics can be used to measure the magnitude of stellar surface shear, the harmonics provide an additional piece of information: the direction of the shear.
Measuring the power ratio of each first harmonic compared with its fundamental frequency, we find that P(f 0 ) = 0.194, while P(f 1 ) = 0.0904.Following Reinhold & Arlt (2015), we assume that spot groups near the stellar equator produce stronger rotation harmonics than spot groups far from the equator.
Our results suggest that f 0 is associated with regions nearer the equator than f 1 , which means KIC 6102338 might exhibit antisolar differential rotation.Indeed, the possibility that P(f 1 ) > P(f 0 ) does not fall within the 95% confidence intervals on Ŝxx (f ).However, the straightforward relationship between P(f ) and spot latitude breaks down when there are large numbers of spot groups.Since the star's rapid rotation suggests that it is young (Barnes 2010) and therefore active (Mamajek & Hillenbrand 2008), the spot patterns may be complex.Further information about the stellar surface (i.e. from Doppler imaging) would be useful for verifying the direction of the shear.Our analysis of the KIC 6102338 Kepler Q3 data demonstrates the value of accuracy and precision in power spectrum estimation: we cannot measure a differential rotation direction using Equation 33 if (a) the confidence intervals are wide or (b) the bias is large.With its bias bounded by Equation 14, the MDMT can be constructed to meet a particular accuracy requirement in Ŝxx (f ); the bias in Ŝxx (f ) in Figure 11 is 443 (e-/s) 2 , which is < 1% of the power in our weakest signal under analysis (harmonic at 2f 1 ).The MDMT should be the preferred Fourier analysis tool for all science goals that require accurate, precise power estimates in addition to period searches.

EPIC 203354381: A COMPLEX ROTATOR WITH
A "SCALLOP-SHELL" LIGHT CURVE From K2 observations of the Upper Scorpius starforming region, Stauffer et al. (2017) identified a new type of stellar variability that produces phased light curves with "scallop shell" shapes.Their hypothesis was that the complex light-curve morphologies result from orbiting material with the same period as the star's rotation (usually < 1 day), possibly due to "slingshotting" of gas in the stellar corona.Stauffer et al. (2018) fol- lowed up with the discovery of more scallop-shell stars in the Upper Sco, ρ Oph, and Taurus star-forming regions.However, both papers confined their light-curve investigations to the time domain.
Here we use multitaper to analyze K2 Campaign 2 observations of EPIC 203354381, an Upper Sco scallopshell star identified by Stauffer et al. (2018).In §5.1 we discuss the time series and taper computation.In §5.2 we identify oscillation modes using the multitaper-F-text combination and discuss their possible physical origins.

Data and tapers
EPIC 203354381 was observed by the Kepler spacecraft in its K2 incarnation from 23 August 2014 to 10 November 2014.For our analysis, we use the detrended light curve produced by the EVEREST pipeline (Luger et al. 2016) with > 3σ outliers clipped.We drop the first 50 observations, which record a stellar flare, from our analysis.The top panel of Figure 12 shows the light curve after outlier and flare rejection.At time BJD = 2456935, marked with a vertical black line, the variability amplitude and light curve morphology suddenly change.We will present a multitaper power spectrum estimate from the entire dataset pictured in Figure 12, plus separate power spectrum estimates from the early data (before BJD = 2456935) and the late data (after BJD = 2456935).
For the full dataset, which has N = 3737, we selected N ϖ = 6.5 and K = 12.Unlike the KIC 6102339 ro- tation signals, the EPIC 203354381 oscillatory modes are well separated in frequency, so we can increase the time-bandwidth product from N ϖ = 4, which we used in §2 and 4, and still be able to distinguish the different modes in a plot of the power spectrum estimate.Recall, however, that the resolution of the F-test does not change with bandwidth.The bottom panel of Figure 12 shows the MDSS for N ϖ = 6.5 and K = 11, with the time series gaps highlighted in gray.Note the tendency of low-order MDSS (small k) to approach zero at the gap edges, a behavior which suppresses spurious high-frequency power produced by sudden observation stoppages ( §3.1).
When we cut the time series in two using the vertical line in Figure 12 and analyze the early data (N = 1944) and late data (N = 1793) separately, we revert to N ϖ = 4, K = 6.The abrupt change in variability amplitude is highlighted in Figure 13, which shows phase-folded versions of the early and late data subsets using the rotation period inferred from the multitaper power spectrum estimates (see §5.2 below).The phased light curve shapes are similar, but the "sharpness" of the late-data light curve suggests that more high-frequency modes are present in the late data than the early data.We will see this prediction borne out in the next section.

Power spectrum estimates
Figure 14 shows the multitaper power spectrum estimate obtained from the entire dataset (top), the F-test p-value (middle), and the multitaper confidence intervals overplotted on the Lomb-Scargle periodogram with the 1/N white-noise false alarm threshold calculated using the Baluev ( 2008) method (bottom).Green stars denote candidate signals for which the null hypothesis probability is < 1/N according to the multitaper / Ftest combination.The multitaper reveals a fundamental mode with f 0 = 1.675 day −1 , which is within 0.1% of the rotation frequency reported by Stauffer et al. (2018).The power spectrum estimate also contains a set of harmonics with f = mf 0 , where m = 2, 3, . . ., 9. Assuming the coronal slingshotting hypothesis is correct, the high-order harmonics reveal that the magnetic structure of the stellar corona is complex, and gas ejections have occurred at many longitudes.The magnetic field must be strong enough to retain some of the ejected gas inside a magnetosphere instead of allowing it to be expelled as a stellar wind (e.g.Morin et al. 2008).There are also two oscillation candidates identified by the F-test that are not part of the rotation harmonic sequence.Their origins are unclear, though we roughly expect one false positive per power spectrum estimate with an F-test pvalue cutoff of 1/N .
Figure 14 highlights the value of a statistical test for harmonic components in the power spectrum that does not assume any particular noise model.Although spikes in the Lomb-Scargle periodogram are visible at the rotation frequency and harmonics m = 2, 3, . . ., 8, only the signals at m = 1 and m = 3 rise above the 1/N false alarm level (though m = 4 is close).Multitaper allows the observers to meaningfully test the hypothesis that oscillatory components are present at high frequencies, even though the redness of the power spectrum ensures that only low-frequency signals can achieve statistical significance according to the white noise assumption.As an aside, although red noise is not obvious on a periodogram plot with a linear y-axis (e.g.Murgas et al. 2023;Stalport et al. 2023), it is easily revealed as a negative slope on a semilog-y plot with frequency on the xaxis (see Appendix A of Dodson-Robinson et al. (2022) for a more detailed explanation of why power spectrum plots should have logarithmic y-axes).
Figure 15 shows multitaper power spectrum estimates from the early (top) and late (bottom) parts of the dataset.While the fundamental rotation frequency plus harmonics m = 2, 3, 4, 7 are detected in the early data, the entire set of harmonics present in Figure 14 does not appear.On the other hand, the later observations taken after BJD = 2456935 contain harmonics m = 2, 3, . . ., 8, 10, 12.All oscillations in the rotation and harmonic sequence have more power than in the early data.There are also candidate oscillations that are not associated with the rotation at f = 5.48 day −1 and f = 5.60 day −1 , plus high-frequency signals with f > 22 day −1 .If these oscillations are real, their origin is unclear.
Figures 12, 13, and 15 could indicate that at BJD = 2456935, the stellar corona slingshotted a parcel of gas at a longitude that was formerly clear.The added complexity of the circumstellar environment manifested in the high-order harmonics present in the power spectrum of the late data.

CONCLUSIONS
We have shown how to apply the missing-data multitaper power spectrum estimator to spacecraft time series with regular observing cadence that have missing observations.The MDMT has the following advantages over the Lomb-Scargle periodogram: • Fewer false positives: By averaging together K independent power spectrum estimates, multitapering reduces the variance of the estimator by a factor of almost 1/K.In addition, stronger serial correlations for classical and missing-data multitaper estimators produce fewer false positives than the (windowed) Schuster periodogram (Thomson & Haley 2014).We find similar improvement over the generalized Lomb-Scargle periodogram in the science examples we show here ( §4 and 5).
• Bias suppression: The broadband bias of a multitaper power spectrum estimate is limited to (1 − λ (K−1) )s 2 x .We urge caution when interpreting Lomb-Scargle estimates: there is no guarantee that leakage from a large (true) oscillation is not causing spurious signals very far away from the main feature.We stress that this fact implies that when the data actually contain one or more large peaks, the number of spurious peaks may be more than is expected for white noise, which can lead to false positives.Furthermore, missing-data multitaper performs better than a periodogram at reproducing the shape and dynamic range of the power spectrum, which is crucial for study- ing broad-spectrum phenomena such as granulation and turbulence.
• Minimal pseudowindowing: Formally the MDMT spectral window is the same for all frequencies (Chave 2019), and in practice we have found that spectral window variation as a function of frequency is limited to ∼ 20%.Our results show that Lomb-Scargle pseudowindows can differ by an order of magnitude ( §3.2).
• True large peaks take on the rectangular shape of the spectral window, which is given by a linear combination of tapers that are optimally large on the interval (−ϖ, ϖ) and small everywhere else.The rectangular shape serves as a visual indicator of a purely sinusoidal oscillation in the power spectrum.
• F-testing: To find statistically significant oscillations, an observer can use the F-test for harmonic components, which works directly on the complex eigencoefficients.The F-test can identify oscillations embedded in a red noise background that would not rise above a white noise-based false alarm threshold, such as the high-frequency modes in the EPIC 203354382 K2 dataset ( §5).
• Confidence intervals: The precision of the MDMT can be quantified, which is important when the science goal requires estimating oscillation power ( §4).
In summary, the missing-data multitaper power spectrum estimator is the best Fourier analysis tool for space missions such that return regularly spaced time series with gaps, such as Kepler, K2, and the upcoming Twinkle mission (Stotesbury et al. 2022).

ACKNOWLEDGEMENTS
The authors would like to thank Alan Chave for discussion on missing-data power spectrum estimators.This material is based upon work supported by the U.S. Department of Energy, 529 Office of Science, under contract number DE-AC02-06CH11357.Additional funding was provided by the Bartol Research Institute.We thank the contribution of an anonymous reviewer for their remarks on an earlier version of this work.
The submitted manuscript has been created by UChicago Argonne, LLC, Operator of Argonne National Laboratory ("Argonne").Argonne, a U.S. Department of Energy Office of Science laboratory, is operated under Contract No. DE-AC02-06CH11357.The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan.http://energy.gov/downloads/doe-public-access-planSOFTWARE All multitaper calculations in this article were performed with the Julia package Multitaper.jl(Haley & Geoga 2020b, https://github.com/lootie/Multitaper.jl).The github repository contations tutorial Jupyter notebooks that demonstrate all multitaper functionality, the notebooks that perform the calculations in this article, and a quickstart guide to Julia for astronomers.For python users, the wrapper multitaperpy is available at https://github.com/lootie/multitaperpy/.This work also made use of LombScargle.jl(Giordano & contribu-tors 2017), Astropy (Astropy Collaboration et al. 2013, 2018), and FINUFFT (Barnett 2020).

DATA AVAILABILITY
The K2 observations of EPIC 203354381 used in §5 can be accessed via 10.5281/zenodo.8436150.
KIC 6102338 Kepler Q0, Q3, and Q16 PDCSAP flux data were accessed from the Mikulski Archive for Space Telescopes (MAST) at the Space Telescope Science In-

Figure 1 .
Figure 1.Kepler Q0 long-cadence observations of KIC 6102338 after replacement of missing data.

Figure 2 .
Figure 2. Top: The six DPSS tapers with the highest spectral concentrations given N = 1639, N ϖ = 4 (ϖ = 0.002441, π = 0.1194 day −1 ).Bottom: The positive half of the K spectral windows W (k) (f ) corresponding to each taper, with a vertical dotted line marking the bandwidth π.The legend at the lower right shows the values of k and λ that correspond to each taper.

Figure 3 .
Figure 3. Left: Distribution of absolute flux errors in the Kepler Q0 observations of KIC 6102338.Right: Distribution of fractional flux errors in the same dataset.

(
Rayleigh 1903).In Figure2, N ϖ = 4-a standard choice in the multitaper literature for time series with N ∼ 1000 (e.g.Thomson 1982;Ojeda & Whitman 2002)-and the main lobe of the spectral window is 2N ϖ = 8 Rayleigh units wide.Next the observer chooses the optimum integer value of K in the range N ϖ ≤ K ≤ 2N ϖ.For example, Figure2could include two more bandlimited DPSS's than are pictured (k = 6, k = 7), but observers are not required to use all DPSS's.The spectral concentration of w (7) t zero at the beginning and end of x n plus the

Figure 4 .
Figure 4. Top: Multitaper power spectrum estimate from the Kepler Q0 observations of KIC 6102338, using the tapers shown in Figure 2. The light blue shading shows the 95% confidence interval.The vertical crossbar shows the expected value of the jackknife variance estimate Var{ Ŝxx(f )}, while the horizontal crossbar has width π. Green dots indicate signals that exceed the 1 − 1/N significance threshold ( §3.5).Middle: The p-value given by the F-test for harmonic components, which is performed at each frequency ( §3.5).Green dots mark the frequencies for which p < 1/N .Bottom: The Schuster periodogram (black) overplotted on the multitaper 95% confidence intervals (orange).The Schuster periodogram has larger variance than the multitaper, which is consistent with the distribution of each estimator (χ 2 2 vs. χ 2 2K ).
Figure 5.Our choice of N ϖ = 4 for both Q0 and Q16 yields a narrower bandwidth ϖ for Q16, which has more observations (N = 3534 for Q16 vs. N = 1639 for Q0).The middle panel shows the p-value associated with the F statistic as a function of frequency ( §3.5).As in Fig-ure 4, the green dots show oscillations significant at or above the 1−1/N level identified by the F-test.The bottom panel shows the Lomb-Scargle periodogram (gray) overplotted on the multitaper 95% confidence interval (orange).
Figure 7 highlights the variance suppression ability of the missing-data multitaper: the 95% confidence interval does not bound the vertical excursions of the Lomb-Scargle periodogram.

Figure 5 .
Figure 5. Top: PDCSAP fluxes from the Kepler Q16 observations of KIC 6102338.Bottom: Missing-data Slepian sequences w (k) tfor the 3534 observations plotted above, constructed with N ϖ = 4 and K = 6.The legend shows the spectral concentration of each W (k) (f ).Gray shading highlights timestamps with missing data.

Figure 6 .
Figure 6.Spectral windows W (k) (f ) from the gpss tapers applied to the Kepler Q16 observations of KIC 6102338 (Figure 5).The legend indicates the spectral concentrations and the vertical dashed line marks the bandwidth.

Figure 7 .
Figure7.Top: The missing-data multitaper power spectrum estimate calculated using the tapers in Figure5.The shaded blue region shows the 95% confidence interval, while the green dots show oscillatory signals that are significant at the 99.9% level according to the F-test.The vertical crossbar has height E{ Var{ Ŝxx(f )}}, while the horizontal crossbar has width π = 0.05539 day −1 .Middle: The p-value given by the F-test for harmonic components, which is performed at each frequency ( §3.5).Green dots mark the frequencies for which p < 1/N .Bottom: The Lomb-Scargle periodogram (gray) overplotted on the multitaper 95% confidence intervals (orange).The Lomb-Scargle periodogram has larger variance than the multitaper, which is consistent with the distribution of each estimator (χ 2 2 vs. χ 2 2K ).

4.
KIC 6102338: DIFFERENTIAL ROTATION, HARMONICS, AND POSSIBLE ANTISOLAR SHEAR Our first missing-data multitaper demonstration uses the Kepler Q3 observations of KIC 6102338, an active M dwarf at a distance of 307 pc (Gaia Collaboration 2020).

Figure 10 .
Figure 10.Top: xn, the Kepler Q3 PDC-MAP light curve of KIC 6102338 after removal of 3σ outliers.Bottom: Missing-data Slepian sequences w (k) t applied to xn.Gray shading highlights missing data and spectral concentrations λ (k) are given in the legend.For this set of tapers, N ϖ = 4.5 and π = 0.0533 day −1 .

Figure 11 .
Figure 11.Power spectrum estimates from the time series xn plotted in Figure 10.Top: The blue line shows the missing-data multitaper with w (k) t from Figure 10, with shaded 95% confidence intervals.Red stars show signals with significance at least 1 − 1/N according to the F-test.The green circle shows an additional signal with 99.9% significance.Annotations give the frequencies of all significant oscillations.Middle: The p-value associated with the F statistic at each frequency.The red stars and green circle have the same meaning as above.Bottom: The Lomb-Scargle periodogram plotted on top of the multitaper confidence intervals, with the 1/N false alarm level indicated by the dotted line.

Figure 12 .
Figure 12.Top: Detrended K2 observations of scallop-shell star EPIC 203354381 with > 3σ outliers clipped.Bottom: Missingdata multitapers for N ϖ = 6.5, K = 11.Spectral concentrations are shown in the legend and gaps are highlighted in gray.

Figure 13 .
Figure 13.Phase-folded versions of the early (green circles) and late (purple stars) parts of the EPIC 203354381 light curve shown in the top panel of Figure 12.The early and late parts of the dataset are before and after BJD = 2456935, respectively.

Figure 14 .
Figure 14.Power spectrum estimates from the time series xn plotted in Figure 12.Top: The blue line shows the missing-data multitaper with w (k) t from Figure 12, with shaded 95% confidence intervals.Green stars show signals with significance at least 1 − 1/N according to the F-test.Middle: The p-value associated with the F statistic at each frequency.Green stars have the same meaning as above.Bottom: The Lomb-Scargle periodogram plotted on top of the multitaper confidence intervals, with the dotted line showing the 1/N false alarm level.

Figure 15 .
Figure 15.Multitaper power spectrum estimates from the early (top) and late (bottom) parts of the EPIC 203354381 K2 light curve.Green stars mark signals with F-test p-value < 1/N .