Magnitude-squared Coherence: A Powerful Tool for Disentangling Doppler Planet Discoveries from Stellar Activity

Welch (1967) built on the work of Bartlett (1948), who proposed estimating the power spectrum of a stationary, regularly spaced time series x_t by dividing the series into segments, computing a periodogram of each segment, and averaging the periodograms:

$$\begin{eqnarray}&&{\overline{S}}_{{xx}}(f)=\displaystyle \frac{1}{K}\sum _{k=0}^{K-1}{\hat{S}}_{{xx}}^{(k)}(f)=\displaystyle \frac{1}{K}\sum _{k=0}^{K-1}\displaystyle \frac{| {\tilde{x}}^{(k)}(f){| }^{2}}{{N}^{(k)}},\end{eqnarray} \tag{ 16 }$$

where K is the number of segments and ${\tilde{x}}^{(k)}(f)$ denotes the NFFT of the data subsequence ${x}_{j}^{(k)}$ of length N^(k). ⁷ For bivariate time series observed at the same time stamps t_j , one uses the same segmentation scheme for both x_t and y_t (i.e., two observables measured from the same astronomical spectrum taken at time t_j are assigned to the same segment k) and computes ${\hat{S}}_{{xy}}^{(k)}(f)$, ${\hat{S}}_{{xx}}^{(k)}(f)$, and ${\hat{S}}_{{yy}}^{(k)}(f)$ for each segment. Averaging together the K different estimates of the cross spectrum and power spectra yields ${\overline{S}}_{{xy}}(f)$, ${\overline{S}}_{{xx}}(f)$, and ${\overline{S}}_{{yy}}(f)$—all the ingredients needed to compute ${\hat{C}}_{{xy}}^{2}$ using Equation (14). Figure 4 shows the Dumusque et al. (2012) α Cen B $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ time series divided into nonoverlapping segments in preparation for using Equation (16).

Zoom In Zoom Out Reset image size

Bartlett (1948) and Welch (1967) explored time-series segmenting, not because they wanted to compute magnitude-squared coherences, but because they were looking for a power spectrum estimator whose variance decreases as the number of observations increases. This is not true of either the standard periodogram (Schuster 1898) or the Lomb–Scargle periodogram (Lomb 1976; Scargle 1982): for both estimators, the variance is independent of the number of observations. Welch (1967) pointed out that nonoverlapping segments (aka Bartlett's method) are ideal for reducing the variance of any $\hat{S}(f)$, as the various ${\hat{S}}^{(k)}(f)$ are fully independent of one another. But there is a catch. Although both the standard periodogram and the Lomb–Scargle periodogram are asymptotically unbiased—i.e., $\hat{S}(f)\to S(f)$ as N → ∞ —the bias can be severe for small or even not-so-small N ("tragedy of the periodogram"; Percival 1994). To understand bias, suppose an RV time series y_t traces a planet in a circular orbit with period P_pl . Since the planet's time domain signature is a perfect sinusoid, ${\hat{S}}_{{yy}}(f)$ should have an infinitely thin delta function at f_pl = 1/P_pl. But the finite duration of y_t creates spectral leakage: in ${\hat{S}}_{{yy}}(f)$, the planet's signal will land mostly on the frequency f_j nearest to f_pl , but there will be nonzero contributions to every frequency in the grid (e.g., Harris 1978). This leakage is responsible for periodogram bias. Thus there is a tradeoff between bias and variance: we want a large number of segments K in order to improve the consistency of ${\hat{S}}_{{xx}}(f)$, ${\hat{S}}_{{yy}}(f)$, and ${\hat{C}}_{{xy}}^{2}(f)$, but if N^(k) is too low our power spectrum and coherence estimates will be consistently biased—that is, ${\hat{C}}_{{xy}}(f)$ will be far away from C_xy (f)—especially at high frequencies (e.g., Podesta 2006; Bronez 1992; Percival & Walden 2020).

To reduce the bias of each ${\hat{S}}^{(k)}(f)$ while retaining most of the variance suppression associated with large K, Welch (1967) proposed using tapered, overlapping segments (tapering is discussed in Section 3.3). In the 50% overlap scheme, we segment x_t as follows:

$$\begin{eqnarray*}\begin{array}{rcl}{x}_{j}^{(0)} & = & {x}_{0}\ldots {x}_{{N}^{(k)}}\\ {x}_{j}^{(1)} & = & {x}_{{N}^{(k)}/2}\ldots {x}_{3{N}^{(k)}/2-1}\\ {x}_{j}^{(2)} & = & {x}_{{N}^{(k)}}\ldots {x}_{2{N}^{(k)}-1}\\ {x}_{j}^{(3)} & = & {x}_{3{N}^{(k)}/2}\ldots {x}_{5{N}^{(k)}/2-1}\\ & \ldots & \\ {x}_{j}^{(K-1)} & = & {x}_{N-{N}^{(k)}}\ldots {x}_{N-1}.\end{array}\end{eqnarray*}$$

Each 50% overlapping segment has 2N/(K + 1) data points, as opposed to N/K for nonoverlapping segments. But because of the overlap, the resulting spectral estimates ${\hat{S}}^{(k)}(f)$ are not independent. The variance reduction associated with overlapping segments is

$$\begin{eqnarray}&&\mathrm{Var}\{\bar{S}(f)\}=\displaystyle \frac{1}{\tilde{K}}\mathrm{Var}\{\hat{S}(f)\},\end{eqnarray} \tag{ 17 }$$

where $\tilde{K}\lt K$ is an effective number of segments. The left panel of Figure 5 depicts the allocation of data points in a sample x_t among segments ${x}_{j}^{(k)}$ in Welch's 50% overlapping segment method.

Zoom In Zoom Out Reset image size

The value of $\tilde{K}$ depends on the type of taper (also called a window) applied to each segment. Tapers are functions w_t that are premultiplied with x_t and y_t to minimize spectral leakage, and are especially valuable for detecting weak signals in the neighborhood of much stronger signals. They are normalized such that ${\sum }_{t=0}^{N-1}{w}_{t}^{2}=1$ so as to conserve power. Figure 6 shows a synthetic RV data set from a star with two unequal-mass planets in circular orbits:

$$\begin{eqnarray}&&{y}_{t}=\cos (2\pi {{tf}}_{1})+0.017\sin (2\pi {{tf}}_{2}),\end{eqnarray} \tag{ 18 }$$

where f₁ = 1.7 days⁻¹ and f₂ = 1.2 days⁻¹. ⁸ When no taper is applied (i.e., the data set retains its rectangular or "boxcar" taper associated with the fact y_t = 0 before and after the observing run such that ${\hat{S}}_{{yy}}(f)=\tfrac{1}{N}| {\sum }_{t=0}^{N-1}({y}_{t}-\bar{y}){e}^{-i2\pi {ft}}{| }^{2}$), one does not detect the signal associated with planet 2 in ${\hat{S}}_{{yy}}(f)$. But when a minimum four-term Blackman–Harris taper (Equation (33) of Harris 1978) is applied to y_t , so that the power spectrum estimate becomes

$$\begin{eqnarray}&&{\hat{S}}_{{yy}}^{w}(f)=| \sum _{t=0}^{N-1}{w}_{t}({y}_{t}-\bar{y}){e}^{-i2\pi {ft}}{| }^{2},\end{eqnarray} \tag{ 19 }$$

planet 2 is detected despite being responsible for only 0.028% of Var{y_t }. Here ${\hat{S}}_{{yy}}^{w}(f)$ is not computed with Welch's algorithm—it is a standard Schuster (1898) periodogram. See Section A.1 for reasons why Figure 6 and all subsequent power spectrum plots have logarithmic y-axes.

Zoom In Zoom Out Reset image size

Viewing tapering from a frequency-domain context, when the time series is evenly spaced with Δt_j = 1, we have the following (Percival & Walden 2020, p. 186):

$$\begin{eqnarray}&&{\rm{E}}\{{\hat{S}}_{{yy}}^{w}(f)\}={\int }_{-1/2}^{1/2}{\left|\sum _{t=0}^{N-1}{w}_{n}{e}^{-i2\pi t(f-g)}\right|}^{2}S(f){df},\end{eqnarray} \tag{ 20 }$$

which is a convolution between the true power spectrum and the spectral window W(f), defined as

$$\begin{eqnarray}&&W(f)={\left|\sum _{t=0}^{N-1}{w}_{t}{e}^{-i2\pi {ft}}\right|}^{2}.\end{eqnarray} \tag{ 21 }$$

That is, the bias of the power spectrum estimator comes from convolving the true power spectrum with the spectral window. Since spectral leakage is created by the "smearing" effect of the Fourier transform of the spectral window, it is best if W(f) resembles a delta function as closely as possible (Harris 1978). For RV data sets, we calculate W(f) using the adjoint NFFT by replacing x_j with w_j in Equation (10). Note that, when generalizing Equation (20) to unevenly spaced time series, one does not strictly obtain a convolution between S(f) and W(f), since, in general, the NFFT cannot simply be inverted (Scargle 1982, Appendix D). However, we will see in Section 4 that peaks in the RV power spectra are shaped like W(f), and will follow Scargle (1982) in referring W(f) as the "spectral window."

When applying Welch's algorithm, tapering the overlapping segments not only mitigates spectral leakage—it also increases the independence of the various ${\hat{S}}^{(k)}(f)$. To see why, we examine the right panel of Figure 5, which shows tapers ${w}_{j}^{(k)}$ applied to each segment depicted in the left panel of Figure 5. At data point t = 5, where ${w}_{j}^{(0)}$ is highest, ${w}_{j}^{(1)}$ is nearly zero. Where ${w}_{j}^{(1)}$ is at its maximum, overlapping tapers ${w}_{j}^{(0)}$ and ${w}_{j}^{(2)}$ are near zero, and so on. The tapers ensure that very little information contained in segment k gets repeated in segments k − 1 or k + 1. The effective number of segments $\tilde{K}$ is

$$\begin{eqnarray}&&\tilde{K}=\displaystyle \frac{K}{1+2{c}^{2}-2{c}^{2}/K},\end{eqnarray} \tag{ 22 }$$

where c is a constant that depends on the type of taper applied (Welch 1967). The boxcar taper belonging to otherwise untapered segments ${x}_{j}^{(k)}$ has c = 0.5, while the minimum four-term Blackman–Harris taper has c = 0.038, yielding $\tilde{K}\approx K$. Figure 7 shows the GJ 3998 RV data of Affer et al. (2016) broken into two 50%-overlapping segments with the associated minimum four-term Blackman–Harris tapers.

Zoom In Zoom Out Reset image size

If all RV planet-search data were observed at evenly spaced time intervals, as in the spacecraft example in Section 2, the benefits of applying tapers to the Welch's segments when computing ${\hat{S}}_{{xx}}(f)$, ${\hat{S}}_{{yy}}(f)$, and ${\hat{C}}_{{xy}}^{2}(f)$ would far outweigh the drawbacks. ⁹ Reexamining Figure 6, we see that ${\hat{S}}_{{yy}}^{w}(f)$—the tapered power spectrum estimate of y_t from Equation (18)—has a dynamic range of over 13 orders of magnitude. This dynamic range would allow planet hunters to identify terrestrial planets, hot Jupiters, rotation, and activity cycles all in the same RV power spectrum estimate ${\hat{S}}_{{yy}}^{w}(f)$: there would be no need to fit the strongest signal, subtract it out, examine a periodogram of the residuals, and keep iterating until no more signals were found. When the temporal cadence is only mildly uneven—as is common in paleoclimatology—the Blackman–Harris window and other similar tapers retain most of their bias-suppression ability (e.g., Ólafsdóttir et al. 2016).

But the wildly uneven observing cadence of RV time series often destroys the tapers (e.g., Scargle 1989). Figure 8 shows the effect of tapering the 51 Peg b RV data set of Butler et al. (2006). The top left shows the published data y_t , while the bottom left shows w_t y_t , where w_t is the minimum four-term Blackman–Harris taper evaluated at the observation time points t_j . The top right shows that the planet's signal (black dotted line) is clearly visible in ${\hat{S}}_{{yy}}(f)$ (blue dashed–dotted line), while on the bottom right, ${\hat{S}}_{{yy}}^{w}(f)$ shows nothing but noise. Almost all of the 51 Peg observations took place at the very beginning and the very end of y_t , when the interpolated w_t is near zero, so tapering removes nearly all the information contained in the time series. When the observing cadence is extremely uneven, we revert to the boxcar tapers when computing ${\hat{C}}_{{xy}}^{2}(f)$. A good rule of thumb is that all tapers should retain the bell-like shapes shown in the right panel of Figure 5. If the "bell" is missing huge chunks or its shape is not recognizable when plotted, tapering may do more harm than good. ¹⁰

Zoom In Zoom Out Reset image size

Another consideration when applying Welch's algorithm to RV data is that large gaps in the middle of a segment should be avoided when possible. Gappy segments are difficult to handle because the resolution limit of ${\overline{S}}_{{xx}}(f)$, ${\overline{S}}_{{yy}}(f)$, and ${\overline{S}}_{{xy}}(f)$ is set by the time duration of the segments:

$$\begin{eqnarray}&&{{ \mathcal R }}^{(k)}=\displaystyle \frac{1}{{t}_{{N}^{(k)}}-{t}_{{0}^{(k)}}},\end{eqnarray} \tag{ 23 }$$

where ${t}_{{N}^{(k)}}$ is the final timestamp in segment k, ${t}_{{0}^{(k)}}$ is the first timestamp in segment k, and $2{{ \mathcal R }}^{(k)}$—called the Rayleigh resolution in analogy to Rayleigh's criterion in optics—is both the smallest oscillation frequency that can be detected and the frequency separation Δf of two barely resolved peaks (e.g., Godin 1972). For evenly spaced data, ${{ \mathcal R }}^{(k)}$ is the same for all segments, but RV data sets always yield varying ${{ \mathcal R }}^{(k)}$. A segment ${x}_{j}^{(k)}$ that contains one or more large gaps has a low apparent Rayleigh resolution, which makes one optimistic that low-frequency information missing from other segments might be available in ${x}_{j}^{(k)}$. But ${\hat{S}}_{{xx}}^{(k)}(f)$ is often misleading at small integer multiples of ${{ \mathcal R }}^{(k)}$ because the gaps make for poor phase coverage of low-frequency oscillations. Here we use the segment with the longest time duration ¹¹ to define the Rayleigh resolution ${ \mathcal R }$ of ${\overline{S}}_{{xx}}(f)$, ${\overline{S}}_{{yy}}(f)$, ${\overline{S}}_{{xy}}(f)$, and ${\hat{C}}_{{xy}}^{2}(f)$,but we urge caution when examining the low-frequency end of each statistic if one or more segments contains a large gap.

In practice it is difficult to keep the Welch's segments gap-free, as seasons and telescope scheduling often combine to isolate a handful of data points from the rest of the time series. The α Cen B $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ time series of Dumusque et al. (2012) shown in Figure 4 would ideally be broken into four segments, but we use only three segments because the first observing season (modified Julian dates 54,550–54,650) recorded just 42 observations—not enough to yield a low-bias periodogram (Pukkila & Nyquist 1985; Springford et al. 2020). No tapers are applied to the α Cen B data because of the gap in segment 1. With nonoverlapping, boxcar-tapered segments, the α Cen B segmenting scheme differs from Bartlett's method only in that each segment has a different number of data points. For all segmenting patterns, the average cross spectrum is weighted by the number of data points in each segment:

$$\begin{eqnarray}&&{\bar{S}}_{{xy}}(f)=\displaystyle \frac{{\sum }_{k=0}^{K-1}{N}^{(k)}{\hat{S}}_{{xy}}^{(k)}(f)}{{\sum }_{k=0}^{K-1}{N}^{(k)}}.\end{eqnarray} \tag{ 24 }$$

Average periodograms ${\bar{S}}_{{xx}}(f)$ and ${\bar{S}}_{{yy}}(f)$ are likewise weighted by N^(k).

To mitigate the worst manifestations of small-sample bias, we recommend that all segments have N^(k) ≥ 100 (Hannan & Nicholls 1977; Pukkila & Nyquist 1985), which requires N ≥ 150 for K ≥ 250% overlapping segments. But many published RV data sets have fewer than 100 observations (a practice we do not endorse). Based on the simulations of Das et al. (2021), who generated small-sample realizations of AR and ARMA processes and compared the resulting periodograms with the processes' analytically known power spectra, we consider N^(k) = 50 a hard lower limit. This means the minimum number of astronomical observations required for Welch's algorithm is 75 (two 50%-overlapping segments each with N^(k) = 50), though more is much better. (In fact, Thomson & Haley (2014) present a time series with N = 1000 for which the periodogram, which records a power-law turbulent cascade, is still biased by more than seven orders of magnitude.) The transformed coherence z(f) and its noise properties can only be described analytically by Gaussian statistics when $\tilde{K}\geqslant 20$ (Enochson & Goodman 1965; Jenkins & Watts 1968), which requires N > 1000. With small N, the false-positive risk may depart from theoretical expectations in unpredictable ways (see below).

While minimizing bias in ${\hat{S}}_{{xx}}^{(k)}(f)$, ${\hat{S}}_{{yy}}^{(k)}(f)$, and ${\hat{C}}_{{xy}}^{2}(f)$ is always an important consideration when deploying Welch's algorithm, an approximate bias correction to ${\hat{C}}_{{xy}}^{2}(f)$ is possible. The bias on the magnitude-squared coherence measurement is

$$\begin{eqnarray}&&\mathrm{bias}\left[{\hat{C}}_{{xy}}^{2}(f)\right]\approx \displaystyle \frac{{\left[1-{\hat{C}}_{{xy}}^{2}(f)\right]}^{2}}{\tilde{K}}\end{eqnarray} \tag{ 25 }$$

(Carter et al. 1973; Bendat & Piersol 2010). The debiased coherence estimate is therefore

$$\begin{eqnarray}&&{\hat{C}}_{{xy}}^{{\prime} 2}(f)={\hat{C}}_{{xy}}^{2}(f)-\mathrm{bias}\left[{\hat{C}}_{{xy}}^{2}(f)\right].\end{eqnarray} \tag{ 26 }$$

All analyses of RV data presented in Section 4 use debiased coherence estimates.

One can also calculate analytical false alarm levels (FALs) for ${\hat{C}}_{{xy}}^{{\prime} 2}(f)$:

$$\begin{eqnarray}&&\mathrm{FAL}=1-{\alpha }^{1/(\tilde{K}-1)}\end{eqnarray} \tag{ 27 }$$

(Carter 1977; Schulz & Stattegger 1997), where FAL is the ${\hat{C}}_{{xy}}^{{\prime} 2}(f)$ threshold associated with false-alarm probability α. Equation (27) gives false alarm thresholds for ${\hat{C}}_{{xy}}^{{\prime} 2}(f)$ given true coherence ${C}_{{xy}}^{2}(f)=0$, i.e., the two time series trace completely unrelated physical phenomena. If a broad-spectrum random process (such as granulation) manifests in both x_t and y_t , then Equation (27) gives artificially low FALs for periodic signals because the two time series have some underlying nonzero coherence at all frequencies. The NWelch software package has an option for calculating frequency-dependent bootstrap FALs for ${\hat{C}}_{{xy}}^{{\prime} 2}(f)$ in addition to using Equation (27). All FALs presented here are calculated using 10,000 bootstrap iterations.

While the Rayleigh resolution (Equation (23)) is the theoretical frequency separation between two barely resolved peaks in ${\hat{S}}_{{xx}}(f)$, ${\hat{S}}_{{yy}}(f)$, or ${\hat{C}}_{{xy}}^{{\prime} 2}(f)$, the true frequency resolution of each statistic is determined by the width of the main lobe in the spectral window W(f). Resolution therefore depends on the choice of taper. In Figure 6, we see that the power spectrum peak associated with planet 1 is quite narrow when the boxcar taper is retained. Although the boxcar taper has poor statistical properties when it comes to leakage and bias, it works well for separating closely spaced signals of similar power. Thus boxcar tapers are useful for asteroseismologists who are trying to resolve modes separated by the small spacing, though they are not optimal for planet hunters who might have to deal with signals of widely varying power. The planet hunter's penalty for using a taper to increase the dynamic range of ${\hat{S}}_{{yy}}(f)$ is lower resolution in all $\hat{S}(f)$. Following Harris (1978), we quantify the half-width of the main lobe in W(f) as the frequency interval ${ \mathcal B }$ over which a sinusoidal signal declines from its peak value by 6 decibels (dB), or a factor of 3.981( ≈ 4). One resolution unit is $2{ \mathcal B }$ wide.

Harris (1978) calculated ${ \mathcal B }$ as a function of ${ \mathcal R }$ for a variety of tapers applied to evenly spaced time series. For the boxcar taper, ${ \mathcal B }=1.21{ \mathcal R }$. ¹² The minimum four-term Blackman–Harris taper has ${ \mathcal B }=2.72{ \mathcal R }$, while the Kaiser-Bessel window (another taper option available in NWelch; see Section Appendix) has ${ \mathcal B }=2.39{ \mathcal R }$. The fact that ${ \mathcal B }$ is a function of ${ \mathcal R }$ means harmonic analysis has a resolution-variance tradeoff in addition to the bias-variance tradeoff discussed in Section 3.1. Users of Welch's method can either prioritize high resolution by constructing a small number of long-duration segments or emphasize false-positive suppression by deploying a larger number of shorter-duration segments. (Of course, small samples and seasonal gaps can constrain the segmenting scheme, making it difficult to find an optimal balance between resolution and variance.) Since almost all RV data sets require segments of varying duration (e.g., Figure 4), we estimate the main lobe half-width of the Welch's estimator as the mean of the main lobe half-widths from the individual segments, weighted by number of points per segment:

$$\begin{eqnarray}&&\hat{{ \mathcal B }}=\displaystyle \frac{{\sum }_{k=0}^{K-1}{N}^{(k)}{ \mathcal B }{\left({ \mathcal R }\right)}^{(k)}}{{\sum }_{k=0}^{K-1}{N}^{(k)}},\end{eqnarray} \tag{ 28 }$$

where ${{ \mathcal R }}^{(k)}$ is defined in Equation (23).

In practice, the Welch's estimator for unevenly spaced data sets has spectral resolution that depends on the exact timing of the observations, not just the segment durations. (The dependence on observation timing applies to the main-lobe width of the Lomb–Scargle spectral window as well.) It is therefore possible for the actual resolution unit to differ from $2\hat{{ \mathcal B }}$ given by Equation (28). NWelch allows the user to examine the specific W(f) associated with any Welch's segmenting and tapering scheme. The software will then empirically calculate ${ \mathcal B }$ by finding the frequency at which $W(f)=\max [W(f)/3.981]$, where $\max [\cdot ]$ is the maximum value.

We recommend examining W(f) not just to find the resolution of a given Welch's estimator, but because changing the segmenting scheme can strongly alter the shape of the spectral window. For example, Figure 9 compares W(f) from the Lomb–Scargle periodogram and the three-segment Welch's estimator applied to the Dumusque et al. (2012) α Cen B data set (Figure 4). As we know from Section 3.3, a periodic signal does not yield a delta function in ${\hat{S}}_{{xx}}(f)$—it instead creates a copy of W(f) centered at the signal frequency. Since the Lomb–Scargle spectral window is beset by "ringing," there can be no clean, isolated peaks in the Lomb–Scargle periodograms of RVs and activity indicators, as noted by (Rajpaul et al. 2016, see also Section 4.2). The Welch's spectral window, which has a well-defined main lobe, is much better suited to identifying periodic signals. Figure 9 suggests that spectral window optimization using Welch's method will be a productive avenue for future research.

Zoom In Zoom Out Reset image size

The last computational method we will describe applies only to power spectrum estimates, not to magnitude-squared coherences, but it is useful for deciding whether a data set contains periodic signals or just noise. Planet hunters often struggle with unrealistic-looking bootstrap false alarm levels (see Section 2.2 of Cumming 2004 for information on how to calculate bootstrap FALs). For example, the top left panel of Figure 10 shows the Hα-index time series ${I}_{{{\rm{H}}}_{\alpha }}$ measured by Robertson et al. (2015) from the Kapteyn's star spectroscopic data set of Anglada-Escude et al. (2014). The time series yields a Lomb–Scargle periodogram with ≥ 10 signals that exceed the bootstrap 1% FAL (Figure 10, top right; periodogram computed with NWelch). But should we really believe that a time series with only 112 measurements can record 10 distinct oscillations? It's more likely that the FALS are misleading (indeed, the same false-alarm threshold problem can be seen in the Lomb–Scargle periodogram of the same data set in Figure 1 of Robertson et al. (2015)), with small-sample statistics, spectral leakage, and aliasing all contributing to the bootstrap failure (e.g., Chernick 2008). The spectral window shown in the bottom panel of Figure 10 has a broad main lobe and a large number of spurious spikes, suggesting that ${I}_{{{\rm{H}}}_{\alpha }}$ periodogram peaks are window function artifacts.

Zoom In Zoom Out Reset image size

In situations like this, we can deploy Siegel's test for compound periodicity. Siegel (1980) developed an extension of Fisher's test, which rejects the null hypothesis of white noise when the maximum power in the normalized periodogram exceeds the critical value g_r . Percival & Walden (2020) approximate g_r as

$$\begin{eqnarray}&&{g}_{r}\approx 1-{\left(\displaystyle \frac{\alpha }{{N}_{f}}\right)}^{1/({N}_{f}-1)},\end{eqnarray} \tag{ 29 }$$

where N_f is the number of entries in the frequency grid and α is the false alarm risk (e.g., α = 0.05 for 5% FAP). While Anderson (1971) notes that Fisher's test is the most powerful identifier of simple periodicity (oscillation at one period only), the test will not work when there are multiple oscillations—such as planet and rotation, rotation and long-term magnetic activity, or rotation with significant power at one or more harmonics. Siegel's test uses a reduced threshold g = λ g_r , where λ < 1, and sums all periodogram power in excess of g to compute the test statistic T_λ :

$$\begin{eqnarray}&&{T}_{\lambda }=\sum _{j=0}^{{N}_{f}-1}\max [0,({\hat{S}}_{{xx}}({f}_{j})-\lambda {g}_{r})].\end{eqnarray} \tag{ 30 }$$

The value of T_λ is then compared with a threshold that depends on the number of entries in the frequency grid. If T_λ exceeds the threshold, the null hypothesis of white noise is rejected and the time series is considered to be periodic. When λ = 0.6, Siegel's test is conservatively optimized for two periodicities, whereas λ = 0.4 is sensitive to three or more periodicities but less robust against noise peaks. For the Kapteyn's star periodogram in Figure 10, Siegel's test does not reject the null hypothesis of white noise, even with λ = 0.4: there is no evidence for periodicity in the time series. We will use Siegel's test in Section 4 to see whether the alternative hypothesis of periodicity is supported for certain time series before using those time series to measure rotation periods. Note that Siegel's test does not differentiate between a smooth power spectrum with a large dynamic range and a power spectrum consisting of white noise plus a single large oscillation—red noise is a failure mode for both bootstrapping and Siegel's test. We are currently incorporating FALs calculated against a red noise model into NWelch and will discuss these in a future publication.

GJ 581 became the third M dwarf known to host a planetary system after Bonfils et al. (2005) used HARPS to discover a Neptune-mass planet with P = 5.36 days. Udry et al. (2007) followed up with reports of planets c (P = 12.93 days) and d (P = 83.6 days), both in or near the habitable zone. Mayor et al. (2009) presented another compelling discovery: planet e (P = 3.15 days, $M\sin i=1.9{M}_{\oplus }$), one of only a handful of super-Earths discovered in multiplanet systems at the time. They also revised the period of planet d to 66.8 days, arguing that the one-year alias introduced by Earth's orbit caused some confusion in the Udry et al. (2007) study. Next, Vogt et al. (2010) used Keck HIRES data to add planets f (P = 433 days) and habitable-zone dweller g (P = 33.6 days) to the inventory, with further follow up by Vogt et al. (2012). The six-planet system with three planets in or near the habitable zone became the object of intense climate modeling efforts (e.g., Heng & Vogt 2011; Pierrehumbert 2011; von Bloh et al. 2011; Wordsworth et al. 2011).

But papers casting doubt on the existence of one or more of the planets began to emerge almost immediately after the report of planets f and g. Anglada-Escudé & Dawson (2010) suggested that planet g was an alias of an eccentricity harmonic of planet d, while Gregory (2011) argued that a Bayesian multiplanet Kepler periodogram could only reliably detect planets b and c. Forveille et al. (2011) questioned the statistical significance of planets f and g after obtaining 121 new HARPS observations. Baluev (2013) demonstrated that planets d, f, and g could be artifacts of red noise with a correlation timescale of 10 days. Finally, Robertson et al. (2014) demonstrated that "planet" d was actually a stellar signal that, when incorrectly modeled in the time domain, created the artifact interpreted as "planet" g. Today the Extrasolar Planets Encyclopedia ¹³ states that GJ 581 hosts only three planets: b, c, and e.

Since Robertson et al. (2014) identified signals d and g as rotation artifacts, we will start by measuring the star's rotation period. Figure 11 shows Welch's power spectrum estimates of the time series used in the Robertson et al. (2014) analysis: ${\hat{S}}_{{xx}}^{w}(f)$ given ${x}_{t}={I}_{{{\rm{H}}}_{\alpha }}$ (top) and ${\hat{S}}_{{yy}}^{w}(f)$ given y_t = RV (bottom). Both power spectrum estimates were computed using three 50% overlapping segments with the minimum four-term Blackman–Harris taper applied to each segment. As with the α Cen B data set, the Welch's estimates have a much cleaner spectral window than the generalized Lomb–Scargle periodograms (Figure 12). Robertson et al. (2014) found a primary peak in the ${I}_{{{\rm{H}}}_{\alpha }}$ Lomb–Scargle periodogram at P₁ = 125 days plus a secondary peak at P₂ = 138 days and attributed the split peak to phase changes in the rotation signal. However, the signals at f₁ = 1/P₁ and f₂ = 1/P₂ are not quite separated by $2{ \mathcal R }$ in the ${I}_{{{\rm{H}}}_{\alpha }}$ generalized Lomb–Scargle periodogram, so cannot be said to be truly distinct. The Welch's power spectrum estimate shows a strong single peak that exceeds the 0.1% FAL at P = 132 days (f = 0.00758 days⁻¹), which we take to be the true rotation period. There are also significant peaks at the first two rotation harmonics, 2f_rot and 3f_rot. Although spectral window of the Welch's estimator is wide, the rotation signal and its harmonics are each separated by more than one resolution unit. The Welch's RV power spectrum has a statistically significant peak at the orbital frequency of planet b (f = 0.186 days⁻¹) plus a local maximum at the frequency of planet c(f = 0.0774 days⁻¹), but no obvious signals at the rotation frequency or its harmonics. The conservative Siegel's test with λ = 0.6 finds a ≥95% chance that both time series are periodic, suggesting that the bootstrap FALs are realistic.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Robertson et al. (2014) used several lines of reasoning to argue that "planets" d and g were really stellar signals: they created a separate fit to RV as a function of ${I}_{{{\rm{H}}}_{\alpha }}$ for each observing season, identified correlations between ${I}_{{{\rm{H}}}_{\alpha }}$ and bisector inverse slope, and calculated a new RV model after subtracting off the best-fit seasonal straight-line models RV(${I}_{{{\rm{H}}}_{\alpha }}$). But a single calculation of ${\hat{C}}_{{xy}}^{{\prime} ;2}(f)$ is enough to reveal the stellar origins of the two signals. Figure 13 shows transformed magnitude-squared coherence z(f) given ${x}_{t}={I}_{{{\rm{H}}}_{\alpha }}$, y_t = RV. Gray bands of width $2{ \mathcal B }$ are centered at the orbital frequencies of planets b and c. The yellow bands, also of width $2{ \mathcal B }$, are centered at the frequencies of "planets" d and g reported by Vogt et al. (2010). Vertical black dotted lines show rotation harmonics f = nf_rot (where n is an integer) near which z(f) exceeds the 5% FAL. "Planet" d sits right atop the first rotation harmonic (f = 2f_rot) at the center of a resolution unit that includes two statistically significant coherence peaks, one of which exceeds the 1% FAL. "Planet" g is within half a resolution unit of the third rotation harmonic (f = 4f_rot), which is also the location of a coherence signal that almost reaches the 0.1% FAL. Significant ${I}_{{{\rm{H}}}_{\alpha }}$-RV coherence can also be seen at other rotation harmonics: at both f = 5f_rot and f = 7f_rot, z(f) exceeds the 0.1% FAL. In contrast, the band centered on planet c is clean. The band centered on planet b includes a signal that rises above the 5% FAL. That signal is likely to be a false positive—indeed, we expect to see more than one false positive above the 5% FAL given that our coherence estimate has a frequency range far greater than 20 resolution units—but further study of the GJ 581 system that incorporates other activity indicators might be useful.

Zoom In Zoom Out Reset image size

The strong signals at the rotation harmonics in Figures 11 and 13 suggest that GJ 581 has a complex rotation signal. Spectral power and coherence at high-order harmonics may result from the star hosting multiple large spots or spot groups (e.g., Rodono et al. 1986; Günther et al. 2020; Perger et al. 2021; Perugini et al. 2021). In the next section, we will investigate coherent stellar signals that do not appear to be associated with the dominant rotation period, but which may be artifacts of differential rotation, giant cells, or supergranulation.

The α Cen B data set assembled by Dumusque et al. (2012) has proven to be useful for identifying spectroscopic signatures of stellar activity. While the star is generally quiet (Cincunegui et al. 2007), it developed one or more large spot groups that caused obvious rotational modulation in the $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ time series from 2010 March 23 to 2010 June 12 (Thompson et al. 2017). This modulation is visible in Segment 2 of Figure 4. Wise et al. (2018) used the 2010 March–June data to find activity-sensitive absorption lines with modulation in either half-depth range or core flux, while Thompson et al. (2017) used the same data to identify pseudoemission features with rotationally driven RV changes. Here we will search for stellar signals using Welch's power spectra of the activity indicators and magnitude-squared coherence between activity indicators and RV.

Figure 14 shows Welch's power spectra of full width at half maximum of the cross correlation between the spectrum and a digital mask (FWHM; Pepe et al. 2000, top left), bisector velocity span (BIS; Toner & Gray 1988, top right), $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ (bottom left), and RV (bottom right) calculated using the segmenting scheme in Figure 4. Before applying Welch's algorithm to the RV data, we removed the binary motion by subtracting the best-fit quadratic model (Endl et al. 2016). To suppress spectral leakage from the long-period activity cycle, each ${x}_{j}^{(k)}$ of FWHM, BIS, and $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ and each ${y}_{j}^{(k)}$ of RV had a linear trend removed before ${\hat{S}}_{{xx}}^{(k)}(f)$ was calculated. The combination of segmenting and detrending mimics the low-pass filtering applied by Dumusque et al. (2012) and the Gaussian process model constructed by Suárez Mascareño et al. (2017). Figure 14 also shows Lomb–Scargle periodograms in light blue. The Welch and Lomb–Scargle periodograms differ at low f partly because the Welch's estimator contains no information about signals with periods longer than the longest segment duration, and partly because of the detrending. ¹⁴

Zoom In Zoom Out Reset image size

Although the Welch's power spectrum of BIS is the only one that has a signal that exceeds the bootstrap white-noise 5% false-alarm threshold, Siegel's test finds that Welch's power spectrum estimates of all of the activity-indicator time series (FWHM, BIS, $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$) show periodicity at the 95% significance level (Section 3.7). We therefore average the frequencies of maximum power in the BIS, FWHM, and $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ periodograms to find the rotation period that best describes the Dumusque et al. (2012) data set: P_rot = 37.7 days (f_rot = 0.0265 day⁻¹), which is consistent with rotation period estimates from the literature (DeWarf et al. 2010; Brandenburg et al. 2017; Suárez Mascareño et al. 2017). ¹⁵ We will use our measured rotation period as a benchmark when searching for stellar signals. Vertical black dotted lines in Figure 14 show shared oscillations identified in the magnitude-squared coherences (see below).

Before we examine any magnitude-squared coherence estimates, we pause to consider Figure 15, which shows scatter plots y_t versus x_t for all $\displaystyle \left(\genfrac{}{}{0em}{}{4}{2}\right)=6$ pairs of observables. RV is not well described as a straight-line function of any activity indicator. We can optimistically hope that's because activity signals are not manifesting in RV, but it is possible that (a) the relationships between RV and activity indicators are nonlinear, phase-lagged (Figure 1), and/or noisy, or (b) $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$, FWHM, and BIS do not provide a complete description of stellar activity. On the other hand, the activity indicators show obvious straight-line relationships with each other—particularly $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ and FWHM—which suggests that they all trace the same underlying physical processes. In a magnitude-squared coherence analysis with two activity indicators as x_t and y_t , we should expect to see statistically significant values in ${\hat{C}}_{{xy}}^{{\prime} 2}(f)$ and z(f).

Zoom In Zoom Out Reset image size

In Figure 16, we see our prediction of high coherences between the activity indicators borne out. The plots show z(f), the $\mathrm{atanh}$-transformed magnitude-squared coherence, given ${x}_{t}=\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$, y_t = FWHM (top); ${x}_{t}=\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$, y_t = BIS (middle); and x_t = FWHM, y_t = BIS (bottom). Solid horizontal lines show false-alarm thresholds calculated with Equation (27), while dotted horizontal lines show bootstrap false-alarm thresholds (Section Appendix). The rotation frequency is marked with a black dashed–dotted line, with the resolution limit indicated in gray. As expected, all three plots show statistically significant z(f) across the entire rotation band. Other frequencies besides rotation at which z(f) exceeds the 1% false-alarm threshold in at least two of three panels in Figure 16 are marked with vertical yellow dashed–dotted lines; if a signal at frequency f is significant in two of three coherences between activity indicators, it means the oscillation is present in all three indicators. Coherent stellar signals are located at f = (0.0460, 0.0985, 0.138, 0.179, 0.275, 0.330, 0.350) days⁻¹, or P = (21.7, 10.2, 7.27, 5.60, 3.64, 3.03, 2.86) days.

Zoom In Zoom Out Reset image size

We know from Figure 9 that the shared signals are not window function artifacts, but their physical origin is unclear: most of them do not fall at simple rotation harmonics (f = nf_rot, where n is an integer), though f = 0.138 days⁻¹ is within half a resolution unit of the n = 4 harmonic. Differential rotation sometimes creates secondary peaks in periodograms of photometry and activity indicators (Reinhold et al. 2013); our best guess is that we are seeing a secondary peak, perhaps along with some beating between differential rotation signals. Supergranulation and giant cells may be in play on the shorter timescales (Nordlund et al. 2009). If we reexamine the Welch's power spectrum estimates of $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$, FWHM, and BIS in Figure 14, we see that the shared signals (again marked by yellow dashed–dotted lines) are all at or near local maxima. Even though none of those local maxima appear statistically significant when compared with our bootstrap white-noise false-alarm thresholds, the underlying oscillations may be real. The downward slopes of ${\hat{S}}_{{xx}}(f)$ for x_t = FWHM, BIS, and $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ suggest that false-alarm thresholds in the α Cen B power spectra should be computed from a red noise model (e.g., Ólafsdóttir et al. 2016). We are developing this functionality and will include it in a future release of NWelch. Red noise has been incorporated in models of RV data by (e.g.,) Baluev (2013); Tuomi et al. (2013), and Feng et al. (2016).

Now we turn to the magnitude-squared coherences between activity indicators and RV. Figure 17 shows z(f) given y_t = RV and x_t = FWHM (top), y_t = RV and x_t = BIS (middle), and y_t = RV and ${x}_{t}=\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ (bottom). The color scheme is the same as in Figure 16. Magnitude-squared coherences in the rotation band are near zero, indicating that no activity indicator is a good tracer of the way rotation manifests throughout the entire duration of the RV data set. In each panel of Figure 17, signals that exceed the 1% FAL in the z(f) estimate shown in that panel only are marked with yellow dashed–dotted lines. The activity indicators' shared signal at f = 0.275 days⁻¹ (Figure 16) also shows up in coherences between FWHM & RV and $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ & RV. FWHM-RV coherence has a second peak at f = 0.167 days⁻¹ (P = 5.99 days, top panel). BIS-RV coherence has peaks at f = 0.115 days⁻¹ and f = 0.288 days⁻¹ (P = 8.70 days and P = 3.47 days, middle panel). The four coherent RV-activity indicator oscillations are marked in the Welch's RV power spectrum in Figure 14 (lower-right panel). The FWHM-RV oscillation at f = 0.167 days⁻¹ coincides with a local maximum of the Welch's RV power spectrum. The $\mathrm{log}{R}_{\mathrm{HK}}^{{\prime} }$ power spectrum also has a local maximum at f = 0.167 days⁻¹, though not the FWHM power spectrum. The coherent RV-BIS oscillation at f = 0.115 days⁻¹ is within one resolution unit of peaks in the FWHM and RV power spectra. All four power spectra have a local maximum at f = 0.275 days⁻¹.

Zoom In Zoom Out Reset image size

Our analysis of the Dumusque et al. (2012) α Cen B data set shows that periodic short-timescale stellar activity occurs at other frequencies besides f_rot and its harmonics. We have identified seven oscillatory signals that are present in all three activity indicators, plus four coherent activity indicator-RV oscillations. Figure 15 demonstrates that the stellar RV signals could not have been diagnosed by fitting a straight-line model to RV as a function of an activity indicator. Accurately identifying the activity signals is an important step toward suppressing false positives in RV planet searches.

Since 2012, the HArps-N red Dwarf Exoplanet Survey (HADES) program has been surveying 78 early-type M dwarfs for Doppler shifts induced by rocky planets (Perger et al. 2017). The first HADES discovery paper featured two planets orbiting GJ 3998 with periods P_b = 2.65 days and P_c = 13.7 days (Affer et al. 2016). Aware of the tendency of M dwarfs to be more active than solar-type couterparts of the same age, the HADES team examined time series of ${I}_{{{\rm{H}}}_{\alpha }}$ and Mt. Wilson S-index measured from the same spectra as the RVs, as well as near-simultaneous EXORAP and APACHE photometry. Affer et al. (2016) identified two signals with periods P₁ = 30.7 days and P₂ = 42.5 days that were present in RV, ${I}_{{{\rm{H}}}_{\alpha }}$, and S-index and posited that P₁ was the true rotation period while P₂ represented modulation due to differential rotation. The photometry was consistent with rotation period P₁. Following Robertson et al. (2014), the HADES team verified that the planetary signals were present in generalized Lomb–Scargle periodograms of RV data from all observing seasons and demonstrated that the Spearman's rank correlation coefficients of RV as a function of S-index and ${I}_{{{\rm{H}}}_{\alpha }}$ were insignificant. They then concluded that the signals at P_b and P_c were planetary in origin.

As with GJ 581 and α Cen B, we begin our analysis by measuring the star rotation period from Welch's power spectrum estimates with the segmentation and tapering scheme shown in Figure 7. Affer et al. (2016) explain how the traditional HARPS and HARPS-N way of measuring RVs, by cross correlating with a mask consisting of a thin rectangle surrounding each spectral line center, is suboptimal for M dwarfs because their spectra feature substantial line blending. Accordingly, in Figure 18 we examine the Welch's and generalized Lomb–Scargle periodograms from the RV, S-index, and ${I}_{{{\rm{H}}}_{\alpha }}$ time series measured with the TERRA pipeline (Anglada-Escudé & Butler 2012). The S-index and ${I}_{{{\rm{H}}}_{\alpha }}$ time series contain statistically significant peaks at 0.0321 days⁻¹ and 0.0311 days⁻¹, respectively, which we average together to find f_rot = 0.0316 days⁻¹/P_rot = 31.7 days. Our rotation period is similar to P₁ identified by Affer et al. (2016). Although the stellar signal with frequency f₂ = 1/P₂ would be more than one resolution unit away from f_rot, the Welch's power spectra do not show any evidence for such a signal.

Zoom In Zoom Out Reset image size

We now turn to magnitude-squared coherence. Although Affer et al. (2016) emphasize that the TERRA pipeline is preferred over the CCF-based HARPS DRS pipeline for M dwarfs, they nevertheless present RVs measured with both methods. They also present a second set of activity indicator time series, with S-index measured according to Henry et al. (1996, H96) and ${I}_{{{\rm{H}}}_{\alpha }}$ measured as in Robertson et al. (2013, R13). Accordingly, we analyze coherences with x_t = TERRA S-index, H96 S-index, TERRA ${I}_{{{\rm{H}}}_{\alpha }}$, and R13 ${I}_{{{\rm{H}}}_{\alpha }}$, and with y_t = TERRA RV and CCF RV—eight measurements of z(f) in all. Figure 19 shows the four transformed coherence measurements with an S-index measurement as x_t . Planets b and c are marked by vertical dashed–dotted black lines surrounded by a shaded region indicating the resolution unit. In three of four S-RV coherences, there is a peak within half a resolution unit of planet c that exceeds the 1% FAL. The fourth, with x_t = H96 S-index and y_t = CCF RV, has a signal at f_c that exceeds the 5% FAL. As with GJ 581, we are also seeing high-order rotation harmonics: all z(f) estimates have high coherence at the fifth harmonic (f = 6f_rot), while z(f) estimates with y_t = CCF RV have high coherence at the third harmonic (f = 4f_rot).

Zoom In Zoom Out Reset image size

If we examine the coherences with an ${I}_{{{\rm{H}}}_{\alpha }}$ measurement as x_t plotted in Figure 20, we find evidence for a stellar signal at the frequency of planet b. When y_t = CCF RV, the ${I}_{{{\rm{H}}}_{\alpha }}$-RV coherence at f_b exceeds the 0.1% FAL. With the TERRA RVs, we find ${I}_{{{\rm{H}}}_{\alpha }}$-RV coherence over the 5% FAL at f_b and nearing the 1% FAL for x_t = TERRA ${I}_{{{\rm{H}}}_{\alpha }}$. However, the band surrounding planet c is clean. ${I}_{{{\rm{H}}}_{\alpha }}$ does not have any coherent oscillations with RV directly at the rotation harmonics, though there is a coherence peak almost within a half resolution unit of 4f_rot for y_t = TERRA RV. The reason Affer et al. (2016) could not see a straight-line relationship between either ${I}_{{{\rm{H}}}_{\alpha }}$ and RV or S-index and RV is because there is more than one oscillation shared between RV and the activity indicators—in addition to the planet candidates, rotation harmonics are present, and so are higher-frequency signals that are not obviously associated with rotation, as in the α Cen B data.

Zoom In Zoom Out Reset image size

While we are not yet ready to state definitively that the GJ 3998 RVs show only stellar signals, the system requires further follow up in light of our results. If coherence between RV and activity indicators at f_b and f_c persists after additional data are taken, that would suggest that the planets are not real. It's also crucial to measure multiple activity indicators, since they may not all trace the same underlying physical phenomena. Indeed, Robertson et al. (2013) highlight the fact that Ca H&K and Hα emission come from different chromospheric depths and point out that signal-to-noise ratio is often problematic in the Ca H&K lines in M dwarf spectra.

What should observers look for when applying bivariate frequency-domain techniques to their own planet-search data sets? In our analyses of GJ 581, α Cen B, and GJ 3998, we roughly followed the following procedure:

• 1.  
  Examine the spectral window of the Welch's estimator. Look for sidelobes that could yield false positives in the power spectra. Calculate the Rayleigh resolution limit, the analytical resolution unit ${ \mathcal B }$ (e.g., Harris 1978), and the empirical value of ${ \mathcal B }$ for a data set's specific observing cadence. NWelch automatically reports these numbers for each segmenting/tapering scheme. Be aware that two signals with a frequency separation of less than $2{ \mathcal B }$ are not statistically distinguishable. It's especially important to consider resolution when searching for planets near the star rotation period or differential rotation.
• 2.  
  Use Welch's power spectra of the activity indicators to estimate the rotation period. While many planet-search targets will have previous rotation period measurements in the literature, the quasiperiodic nature of the rotation signal means one might recover a somewhat different rotation period depending on the observational epoch and activity cycle phase (Robertson et al. 2014). It's important to determine how rotation is manifesting in the particular data set under analysis. For the α Cen B data set, in which three activity indicators yielded three slightly different rotation period estimates, we averaged together the various estimates (Section 4.2). It's possible that the community will come up with a more sophisticated way to handle differing rotation period measurements.
• 3.  
  Search for coherence between RV and activity indicators within one resolution unit of the rotation period, the periods of planet candidates, and the harmonics of all of the above. For any planet that is not large enough or close enough to the star to trigger stellar activity, ${\hat{C}}_{{xy}}^{{\prime} 2}(f)$ and z(f) should be impervious to eccentricity harmonics. Observers who find coherent RV-activity indicator signals at either their planet candidate frequency or its harmonics should use extreme caution before reporting a planet discovery. It's unclear how high up the overtone sequence we should expect to see rotation harmonics—for example, if there's a peak in z(f) near the eighth harmonic (9f_rot), is it truly rotation-related, or is its proximity to a harmonic just a coincidence? We hope stellar physicists will explore the complexity of rotation signals and the extent to which their harmonics should be traceable by RV data sets. We also hope stellar physicists will weigh in on possible sources of coherent RV-activity signals that are not related to rotation.
• 4.  
  When coherent RV-activity indicator signals are found, check to see whether they line up with local maxima in the power spectra. Sometimes a power spectrum peak that's not significant when judged against a white noise model indicates a real signal. The signal processing literature features more sophisticated ways of checking the significance of power spectrum peaks, such as F-testing (e.g., Thomson 1994), false-alarm thresholds determined from red noise models (e.g., Ólafsdóttir et al. 2016), and prewhitening followed by checking against the exponential quantiles. But these techniques are mostly unexplored for RV data. When a peak in z(f) corresponds to a local maximum in the RV and/or activity-indicator power spectrum, it lends support to the hypothesis that the coherence is real and not spurious. However, given that the geophysics literature features coherent signals that do not correspond to power spectrum local maxima (see example in Pardo-Igúzquiza & Rodríguez-Tovar 2012), one should not be too quick to dismiss any signals in z(f) as false positives.

NWelchdefaults to a logarithmic y-axis for all power spectrum plots. The user can select a linear y-axis (see the demo notebooks in the githubrepository for instructions), but there are good reasons to use logarithmic scaling (e.g., Thomson 1994):

• 1.  
  Power spectrum peaks that are not significant when judged against a white noise model can correspond to coherence peaks that are significant, so it is good to be able to see all the local maxima in the power spectrum (see the α Cen B analysis in Section 4.2);
• 2.  
  Human vision and hearing respond to the intensity of stimuli on a logarithmic scale;
• 3.  
  A downward slope in ${\mathrm{log}}_{10}[{\hat{S}}_{{xx}}(f)]$, which is easy to spot on a semilog-y plot, is an indicator of red noise;
• 4.  
  Multiplicative effects appear additive on a logarithmic scale, so convolution in the time domain—which is equivalent to multiplication in the frequency domain—yields simple addition in ${\mathrm{log}}_{10}[{\hat{S}}_{{xx}}(f)]$.

Finally, the convention in the signal processing literature is to use semilog-y plots for power spectra because engineering processes often use linear filters such as moving averages or Gaussian smoothing. A moving average applied to time-domain data is obvious on a semilog-y plot of the power spectrum because the plot will show the filter sidelobes. We recommend that all astronomers working in the frequency domain examine semilog-y plots of their power spectra.

The workhorse calculation that underlies all NWelch functionality is the nonuniform fast Fourier transform (NFFT). Here we provide a brief overview of the mathematics behind the NFFT. Lomb (1976) and Scargle (1982) were the first to develop periodogram (Fourier spectrum) analysis techniques that were quite powerful for finding and testing the significance of weak periodic signals in otherwise random, unevenly sampled data. Given a set of data values x_j = 1, ..., N at respective observation times, the Lomb–Scargle periodogram is constructed as follows. First, compute the data's mean and variance by

$$\begin{eqnarray}&&\bar{x}=\displaystyle \frac{1}{N}\sum _{j}{x}_{j};\quad {\hat{\sigma }}_{x}^{2}=\displaystyle \frac{1}{N-1}\sum _{j}{\left({x}_{j}-\bar{x}\right)}^{2}.\end{eqnarray} \tag{ A1 }$$

Second, for each angular frequency ω = 2π f > 0 of interest, compute a time-offset τ by

$$\begin{eqnarray}&&\tan 2\omega \tau =\displaystyle \frac{{\sum }_{j}\sin 2\omega {t}_{j}}{{\sum }_{j}\cos 2\omega {t}_{j}}.\end{eqnarray} \tag{ A2 }$$

Third, the Lomb–Scargle normalized periodogram (spectral power as a function of angular frequency ω = 2π f) is defined by

$$\begin{eqnarray}&&\begin{array}{r}\begin{array}{l}{P}_{N}(\omega )=\displaystyle \frac{1}{2{\sigma }^{2}}\left[\displaystyle \frac{{\left[{\sum }_{j}({x}_{j}-\bar{x})\cos \omega ({t}_{j}-\tau )\right]}^{2}}{{\sum }_{j}({x}_{j}-\bar{x}){\cos }^{2}\omega ({t}_{j}-\tau )}\right.\\ \,+\,\left.\displaystyle \frac{{\left[{\sum }_{j}({x}_{j}-\bar{x})\sin \omega ({t}_{j}-\tau )\right]}^{2}}{{\sum }_{j}({x}_{j}-\bar{x}){\sin }^{2}\omega ({t}_{j}-\tau )}\right].\end{array}\end{array}\end{eqnarray} \tag{ A3 }$$

The constant τ makes p_n (ω) completely independent of shifting all the t_j by any constant. Lomb (1976) showed that this particular choice of offset has another, deeper, effect: it makes Equation (A3) identical to the equation that one would obtain if one estimated the harmonic content of a data set, at a given frequency ω by linear least-squares fitting to the model

$$\begin{eqnarray*}&&\begin{array}{l}x(t)-\bar{x}=A\cos \omega t+B\sin \omega t.\end{array}\end{eqnarray*}$$

Lomb–Scargle periodograms can be approximated by employing an "extirpolation" process (interpolation of the unevenly spaced points onto a regular grid) using the Lagrange interpolation method, and then using the ordinary FFT (Press & Rybicki 1989). However, as described in Leroy (2012; Springford et al. 2020, Section 7), one can employ the adjoint nonuniform FFT (Keiner et al. 2009) for a fast implementation which is not approximate. Briefly, one computes

$$\begin{eqnarray}&&{S}_{y}=\sum _{i}({x}_{i}-\bar{x})\sin \omega {t}_{i},\quad {C}_{y}=\sum _{i}({x}_{i}-\bar{x})\cos \omega {t}_{i},\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&{S}_{2}=\sum _{i}\sin 2\omega {t}_{i},\quad {C}_{2}=\sum _{i}\cos 2\omega {t}_{i},\end{eqnarray} \tag{ A5 }$$

then

$$\begin{eqnarray}&&\sum _{i}({x}_{i}-\bar{x})\cos \omega ({t}_{i}-\tau )={C}_{y}\cos \omega \tau +{S}_{y}\sin \omega \tau ,\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&\sum _{i}({x}_{i}-\bar{x})\sin \omega ({t}_{i}-\tau )={S}_{y}\cos \omega \tau -{C}_{y}\sin \omega \tau ,\end{eqnarray} \tag{ A7 }$$

$$\begin{eqnarray}&&\sum _{i}\cos 2\omega ({t}_{i}-\tau )=\displaystyle \frac{N}{2}+\displaystyle \frac{1}{2}{C}_{2}\cos 2\omega \tau +\displaystyle \frac{1}{2}{S}_{2}\sin 2\omega \tau ,\end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray}&&\sum _{i}\sin 2\omega ({t}_{i}-\tau )=\displaystyle \frac{N}{2}-\displaystyle \frac{1}{2}{C}_{2}\cos 2\omega \tau -\displaystyle \frac{1}{2}{S}_{2}\sin 2\omega \tau .\end{eqnarray} \tag{ A9 }$$

NWelch uses FINUFFT, the Flatiron Institute nonuniform fast Fourier transform, ¹⁶ to compute the adjoint NFFT. Methods for performing the exponential sums are given by Barnett et al. (2019) and Barnett (2020). For all power spectrum estimates (Figures 11, 14 and 18), NWelch uses the power spectral density normalization:

$$\begin{eqnarray}&&\mathrm{Var}\{{x}_{t}\}=\sum _{m=0}^{{N}_{f}-1}{\hat{S}}_{{xx}}({f}_{m}){\rm{\Delta }}f,\end{eqnarray} \tag{ A10 }$$

where m is the integer index of the frequency grid and Δf is the frequency grid spacing.

The most basic NWelch task is to compute a periodogram of an unevenly spaced time series without tapering, detrending, or segmentation. Here we show that this task produces results that are consistent with astropy.timeseries.LombScargle. The top panel of Figure 21 shows two periodograms of the residuals of the Dumusque et al. (2012) α Cen B RV data after subtracting a quadratic model of binary motion. The periodogram plotted with a dark blue, solid line comes from NWelch, while the periodogram plotted with the light blue, dashed line comes from astropy.timeseries.LombScargle. The periodograms are almost identical, as are their bootstrap 5% and 1% false-alarm levels (purple and green, respectively, with dotted lines for NWelch FALs and solid lines for astropy.timeseries.LombScargle FALs). For simple Lomb–Scargle periodograms with bootstrap false-alarm thresholds, NWelch and astropy.timeseries.LombScargle can be used interchangeably. However, NWelch does not include the Baluev (2008) method for calculating false-alarm probabilities based on extreme value theory.

Zoom In Zoom Out Reset image size