THE RADIAL VELOCITY DETECTION OF EARTH-MASS PLANETS IN THE PRESENCE OF ACTIVITY NOISE: THE CASE OF α CENTAURI Bb

Pre-whitening is a commonly used tool for finding multi-periodic signals in time series data. This method sequentially finds the dominant Fourier components in a time series and removes them. Traditionally it is employed to derive the frequency spectrum of oscillating stars (e.g., García Hernández et al. 2009), but it also has applications as a means of removing the intrinsic stellar variations due to activity. The mathematical foundation for this is that sines and cosines form a set of basis functions. You can represent most functions as a sum of sine waves with different periods and amplitudes. In fact, the DFT of a time series merely gives you the amplitude as a function of frequency of all the sine functions that are present, including those due to noise. The trick is to use enough Fourier components to represent adequately the function of interest (activity in this case) without introducing spurious periods, or altering the amplitude and frequency of a signal you are trying to detect. So long as artifacts introduced by the multi-component sine fit have frequencies and amplitudes different from our signal of interest, pre-whitening can be a useful tool for detecting weak periodic signals. In this case we have the advantage in that we have a priori information about the signal we are trying to detect.

The pre-whitening process is similar to the harmonic analysis used by of D2012 in that both fit the activity signal with multi-sine components. There are, however, two major differences. First, in pre-whitening the time span is not restricted to a few rotational periods, but in our case we use a much longer time span. Second, the frequencies that are removed by pre-whitening are not restricted to just the rotational frequency and its harmonics. The strongest peak in the DFT is removed regardless of its frequency. One can consider harmonic analysis as a more restricted version of pre-whitening.

The case of CoRoT-7b demonstrated the effectiveness of Fourier component analysis in fitting an activity signal. The RV activity jitter for CoRoT-7 was about a factor of 10 higher than for α Cen B. The pre-whitening process resulted in a K-amplitude of 5.5 ± 0.3 m s⁻¹ (Hatzes et al. 2010). This was consistent with the value of K = 5.27 ± 0.81 m s⁻¹ determined using an entirely different filtering approach that did not rely on periodic functions (Hatzes et al. 2011). On the other hand, the initial K-amplitude derived from harmonic analysis was 1.9 m s⁻¹, or almost a factor of three lower than the final value (Queloz et al. 2009). Furthermore, the harmonic analysis failed to detect the presence of a third planetary companion, CoRoT-7d (Hatzes et al. 2010).

In the pre-whitening process one usually picks the highest peak in the Fourier amplitude spectrum of the time series. A least-squares sine fit is made to the data to determine the optimal frequency, amplitude, and phase. This fit is then subtracted from the data, which also removes (or at least minimizes) the effects of alias peaks caused by this signal. A Fourier analysis on the residuals then finds the next dominant peak, and the process is continued until one reaches the noise level of the amplitude spectrum. Each process thus "whitens" the data in frequency space for the next step. The resulting frequencies that are found should represent the dominant frequency components of the time series. In our case we use the sum of the pre-whitened components we have found to provide us with a fit to the activity variations.

The Fourier components one derives (amplitude, frequency, and phase) depend on the length of the time series and the presence of time gaps (sampling window). Some of the effects of these can be explored by analyzing first the complete data set, and then subsets divided into epochs of the measurements. In doing so one may derive slightly different Fourier components in fitting the underlying activity variations. A robust signal should be relatively insensitive as to how we pre-whiten the data.

3.1.1. Fourier Component Analysis of the Full Data Set

The pre-whitening procedure was applied to the full RV data using the program Period04 (Lenz & Breger 2004). For clarity we only show the DFT of the un-whitened data (top panel Figure 1) rather than each step of the process. In the figure we have marked the frequencies found by the pre-whitening process. Note that removing a peak also removes its alias in Fourier spectrum whose amplitude may be higher than a real peak. The amplitude spectrum is dominated by a forest of peaks in the frequency range 0 < ν < 0.1 day⁻¹. These may be due to the rotational frequency of the star, its harmonics, as well as the spectral window (sampling). The highest peaks correspond to a frequency ν = 0.025 day⁻¹ (period, P = 38 days) which is interpreted as the rotational frequency of the star. This frequency has an amplitude of ≈1.5 m s⁻¹, which is about three times larger than the K-amplitude of the purported planet. Note that the orbital frequency at ν = 0.309 day⁻¹ is hardly visible in the original un-whitened data.

Zoom In Zoom Out Reset image size

Table 2 lists the frequencies, corresponding periods, amplitudes, and phases for all signals found in the RV time series. Note that amplitudes found in Table 2 may differ slightly from those seen in Figure 1 due to removal of peaks and their aliases. The lower panel of Figure 1 shows the final pre-whitened amplitude spectrum with only the peak near the orbital frequency of α Cen Bb (f₉) remaining. Note that in fitting the data simultaneously with all frequencies a slightly higher amplitude (≈1.9 m s⁻¹) for the rotational frequency results when compared to the initial DFT. Removing all frequencies in Table 2 from the RV data results in a standard deviation of 1.17 m s⁻¹. Interestingly, this value of σ is consistent with the rms scatter of Epoch 1 for which D2012 noted was a time when the star was relatively inactive and thus should have a lower RV jitter. We thus use σ = 1.2 m s⁻¹ as our "worst-case" estimate of the RV error.

Table 2. Pre-whitening Results for the RV Data

N	Frequency	Period	K-amplitude	Phase
	(day−1)	(days)	(m s−1)
f1	0.02558 ± 0.00002	39.09 ± 0.03	1.89 ± 0.08	0.00 ± 0.01
f2	0.00131 ± 0.00005	763.36 ± 29.12	0.69 ± 0.08	0.79 ± 0.01
f3	0.08163 ± 0.00004	12.25 ± 0.006	1.00 ± 0.08	0.49 ± 0.01
f4	0.10438 ± 0.00005	9.58 ± 0.005	0.75 ± 0.08	0.13 ± 0.02
f5	0.00603 ± 0.00004	165.83 ± 1.10	0.97 ± 0.08	0.78 ± 0.02
f6	0.06633 ± 0.00005	15.79 ± 0.01	0.71 ± 0.08	0.57 ± 0.02
f7	0.03321 ± 0.00005	101.11 ± 0.15	0.67 ± 0.09	0.62 ± 0.02
f8	0.07841 ± 0.00005	12.75 ± 0.047	0.77 ± 0.08	0.05 ± 0.03
f9	0.30906 ± 0.00009	3.2356 ± 0.0001	0.40 ± 0.08	0.33 ± 0.03

Download table as:  ASCIITypeset image

The low-frequency component (f₂) has a period, 763 days, that is shorter than the time span of the observations of 1230 days. This is most likely due to some variation of the activity, but we cannot exclude with certainty that it may be caused by slight differences between the parabolic fit and the true Keplerian orbit. Unfortunately due to the long period this orbit is poorly known. Because the period is over a factor of 200 greater than the planet orbital period this should not effect the subsequent analysis, particularly for the local trend fitting analysis (see below). All errors are uncorrelated errors from Period04; correlated errors are most certainly larger. Many frequencies can be identified with the rotational frequency, ν_rot, or its harmonics: f₁ = ν_rot, f₃, f₈ ≈ 3ν_rot, f₄ = 4ν_rot, and f₆ ≈ 2ν_rot. This gives some justification to the harmonic analysis used by D2012. Note that the last entry (f₉) in the table corresponds to the orbital frequency of α Cen Bb. The pre-whitened period (P = 3.2356 ± 0.0001 days) is identical to the value P = 3.2357 ± 0.0008 days found by D2012. The amplitude, K = 0.40 ± 0.08 m s⁻¹, is a bit lower, but still consistent within the errors with the value of K = 0.51 ± 0.04 m s⁻¹ from D2012. The implied companion mass from the pre-whitened amplitude is m = 0.89 ± 0.18 M_⊕.

The statistical significance of the 3.24 day signal was assessed with a Scargle periodogram analysis (top panel of Figure 2) of the residual RV data produced after removing the first eight frequencies listed in Table 1. The peak at the planet orbital frequency of 0.309 day⁻¹ (P = 3.24 days) appears to be significant due to its high Scargle power (z ≈ 12). The FAP of this peak was assessed using the bootstrap randomization process. The residual RV values were randomly shuffled while keeping the time values fixed. The highest peak in the Scargle periodogram in the frequency range 0.0001 < ν < 0.5 day⁻¹ was found for each random data set. The number of instances where the shuffled data produced power higher than the observed power provided a measure of the FAP. The resulting FAP was ≈0.004. (Note that all FAP values given below are the result of bootstrap analyses performed with 200,000 shuffles.)

Zoom In Zoom Out Reset image size

It is difficult to assess the significance of a periodic signal in the presence of other signals (in this case at least eight) that have been filtered from the data. The FAP calculated with a bootstrap depends on the scatter in the data. After removing a periodic signal the amount of scatter in the data is reduced and a much lower FAP may thus result. However, one cannot be sure that a peak in the amplitude spectrum is a true signal and should be removed (e.g., rotation), or a noise peak that should remain in the data when performing an FAP analysis. The bootstrap analysis may produce an unrealistically low FAP simply because you have "cleaned" the data by lowering the noise floor in Fourier space. An alternative approach is to use the unfiltered amplitude spectrum itself. Kuschnig et al. (1997) established that peaks in the amplitude spectrum that have a height 3.6 times the surrounding noise level correspond to an FAP ≈ 1%.

Period04 calculates a noise level at the orbital frequency of the planet of 0.23 m s⁻¹. The amplitude of the planet signal is ≈0.4 m s⁻¹, which is only a factor of 1.7 above the noise level. This corresponds to an FAP of ≈100%. If one uses the pre-whitened amplitude spectrum with all frequencies in Table 2 removed except for the planet orbital frequency, the noise level is 0.12 m s⁻¹. This is the same value as simply taking the mean amplitude of peaks over the interval 0.25 < ν < 0.35 day⁻¹. This amplitude is 3.3 times the noise level, which corresponds to an FAP ≈ 5%. The noise level of the DFT amplitude spectrum thus indicates an FAP of approximately a few percent, but with a large uncertainty. Given the complexity of the RV variations it is difficult to assess an accurate FAP.

The efficacy of the Fourier procedure was tested on simulated data. First, the 3.236 day period Fourier component (last entry in Table 2) was removed from the HARPS RV data on the assumption that this signal is real. The orbital solution of D2012 was then added back into the data. We also added synthetic planet signals at slightly different orbital periods. For these other simulations we used the RV data without removing the 3.236 day component on the assumption that this is due purely to noise in the data. In this way all simulations preserved the noise characteristics of the real data. As an example, the lower panel in Figure 2 shows a Scargle periodogram of the residual RV data with an artificial planet inserted prior to filtering the data. The planet signal was taken as a simple sine wave with a period of 3.29 days (slightly different from the period of α Cen Bb) and an amplitude of 0.5 m s⁻¹.

Table 3 summarizes the results of the Fourier filtering of the simulated data. The third entry is the simulation using the orbital parameters from D2012. The table lists the input period of the sine function, P_in, the input amplitude, K_in, the period recovered by the pre-whitening procedure, P_out, the output amplitude K_out, and the FAP of the detected signal computed using a bootstrap. In all cases, the input period and amplitude are recovered well and at high level of significance.

Table 3. Tests of the Fourier Procedure on the Full Data Set

Pin	Kin	Pout	Kout	FAP
(days)	(m s−1)	(days)	(m s−1)
3.350	0.50	3.349 ± 0.001	0.54 ± 0.08	<5.0 × 10−5
3.300	0.50	3.299 ± 0.001	0.64 ± 0.09	<5.0 × 10−5
3.236	0.50	3.236 ± 0.001	0.52 ± 0.08	1.5 × 10−5
3.200	0.50	3.203 ± 0.001	0.54 ± 0.10	2.0 × 10−5
3.150	0.50	3.151 ± 0.001	0.58 ± 0.10	<5.0 × 10−5

Download table as:  ASCIITypeset image

3.1.2. Fourier Component Analysis of the Individual Epochs

The pre-whitening procedure was then applied to the individual epochs in Table 1. In cases where you have periodic signals that are evolving with time (e.g., the birth, decay, evolution, and migration of surface spots) pre-whitening of a long time series with large gaps may produce poorer results. A much better fit to the activity could be obtained by using data covering a shorter time interval. Furthermore, since one is deriving a different set of sine functions for filtering the data, this is an independent check on how robust the signal that was found when analyzing the full data set is.

Figure 3 shows the fit to the activity variations to each of the epochs, and Table 4 lists the sine parameters. The first subscript in the frequency identification (ID) refers to the epoch number from Table 2. The DFTs of the RVs from the individual epochs are shown in Figure 4. In the figure we have marked the frequencies found by pre-whitening the data. In Epoch 2 the second dominant peak was found at 0.356 day⁻¹ (P = 2.8 days). This was not removed since it has a frequency near that of the planet orbital frequency.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 4. Pre-whitening Results for the Epoch Data

ID	Frequency	K-amplitude	Phase	Comment
	(day−1)	(m s−1)
f11	0.020 ± 0.001	0.89 ± 0.29	0.91 ± 0.04	Epoch 1
f21	0.0043 ± 0.001	2.01 ± 0.15	0.27 ± 0.02	Epoch 2
f22	0.0692 ± 0.001	0.76 ± 0.15	0.82 ± 0.03	Epoch 2
f23	0.1347 ± 0.001	0.71 ± 0.15	0.07 ± 0.03	Epoch 2
f31	0.0104 ± 0.0002	1.73 ± 0.16	0.33 ± 0.02	Epoch 3
f32	0.0284 ± 0.0003	2.52 ± 0.16	0.88 ± 0.02	Epoch 3
f33	0.0644 ± 0.0016	0.76 ± 0.16	0.38 ± 0.04	Epoch 3
f34	0.1510 ± 0.0016	0.69 ± 0.16	0.60 ± 0.04	Epoch 3
f41	0.0267 ± 0.0005	1.59 ± 0.15	0.21 ± 0.02	Epoch 4
f42	0.0818 ± 0.0005	1.63 ± 0.13	0.65 ± 0.02	Epoch 4
f43	0.1028 ± 0.0014	1.63 ± 0.14	0.73 ± 0.04	Epoch 4

Download table as:  ASCIITypeset image

The upper panel of Figure 5 shows the periodogram of the total RV residuals from all epochs. The epoch pre-whitened RV residuals have a standard deviation of σ = 1.16 m s⁻¹ with the planet signal present. One can still see a peak at the orbital frequency of the planet, but the power is reduced from that found by the analysis on the full data set. A bootstrap analysis yields an FAP of only 0.07 for this peak. Removing the planet signal reduces the standard deviation slightly to σ = 1.13 m s⁻¹.

Zoom In Zoom Out Reset image size

The epoch pre-whitening was also tested on simulated data. A sine wave (P = 3.29 days, ν = 0.304 day⁻¹, K = 0.5 m s⁻¹) was inserted in the data prior to pre-whitening. The lower panel of Figure 5 shows the Scargle periodogram of the residual RV data. The FAP of this signal is <5 × 10⁻⁶ based on a bootstrap analysis. As another test a sine fit to the 3.24 day period was made to the full data set and removed (i.e., removal of f₉ from Table 2) and the orbital solution from D2012 inserted back into the data. The pre-whitening procedure was then applied to the epoch data. The planet signal was detected with an FAP = 0.01%.

The pre-whitening process on the epoch data was able to detect the planet, but with much reduced significance. One would naively think that since essentially the same technique is applied to both the full and epoch data we should arrive at the same answer. Indeed, we removed from the data the 3.24 day period found in pre-whitening the full data set and re-inserted the orbital solution for α Cen Bb from D2012. Pre-whitening of the epoch data showed significant Scargle power of z ≈ 15, which corresponds to FAP ≈ 0.03%. This is consistent with the full data set pre-whitening: z ≈ 18, FAP ≈ 0.002%. Both approaches to filtering the data detect the planet with a much higher significance than was found. The discrepancy is the first hint that the planet signal may depend sensitively on how the data are filtered. Clearly, an independent, non-Fourier-based technique is needed to help resolve this discrepancy.

Although Fourier pre-whitening is a useful tool for getting a quick result, it has its drawbacks. Some functions, like a linear trend, may require a large number of Fourier components (i.e., free parameters) to fit them when simpler functions with fewer free parameters such as low order polynomials could provide a better fit. Furthermore, with multi-sine components, if one uses insufficient Fourier components, there can be mismatch between the true activity variations and the fit. The sampling window complicates matters further. All of these may introduce false peaks in the filtered data, or increase the apparent significance of noise peaks. It is therefore important to check the results of pre-whitening and harmonic analysis by using alternative filtering techniques that rely less on periodic functions and that have fewer free parameters. For a robust signal, different filtering approaches should produce comparable results as was the case for CoRoT-7b.

A better filter should exploit the fact that we know the periods of interest, namely the orbital period of the planet as well as the rotation period of the star. If we can fit the activity variations over a much shorter time span, the signal should be more coherent and stable. We thus should be able to use functions with fewer parameters that fit better the activity variations at that particular time as opposed to using a global fit that requires more free parameters (sine functions). Trend filtering is often employed to remove the stellar variations when searching for transit signals in light curves (Kovács et al. 2005; Grziwa et al. 2012; Bakos et al. 2013).

The time window for fitting the activity variations is bounded by two limits. The lower limit is the orbital period of the planet. Filtering out activity variations over a time span less than this runs the risk of suppressing any real RV variations due to the planet. The upper limit is defined by the rotation period of the star, a time span over which we consider the activity variations likely to be stable and coherent.

The RV data were visually inspected and divided into time chunks for the local trend fitting. The following criteria were used to decide which data went into a specific chunk: (1) the time interval of a chunk, ΔT, should cover as many cycles of the planet orbit as possible but at least a full orbital period, but less than a stellar rotation period: P_orbit < ΔT < P_rot. (2) The time series should have good sampling (preferably nightly) in the interval and with no gaps longer than a day or two. (3) In fitting trends the data should be grouped to avoid large gags in temporal coverage or abrupt changes in the long-term variations in the chunk. (4) If successive measurements were separated by several days, but seemed to follow the overall trend, they were kept in the analysis. If they showed significant departures from the trend that required more complicated fitting functions (i.e., a high order polynomial with more inflections), they were removed. In total 41 data points were removed from the original HARPS data. In short, the data were grouped in subsets that showed smooth variations of the underlying trend. The fit to these should have different Fourier components that are far removed in frequency from those of the planetary signal.

The choice of fitting function for the trend in each chunk was determined visually so that the fit to the underlying activity variations would have a minimal influence on the shorter periodic variations of the planet. If the data in a chunk showed no long-term variations, the average value was subtracted. If it showed a linear trend, a least-squares linear fit was made. In cases where the chunk trend showed curvature a second or third order polynomial was used. Although we tried to avoid using periodic functions, in one chunk the underlying activity variations could best be fit with a multi-sine component (see below).

Two time intervals RJD = 54,935–54,955 and RJD = 55,672–55,692 showed significant periodic variations. In the first of these intervals the RV data centered on RJD = 54,939.68–54,941.85 showed an additional linear trend on top of the sinusoidal variations (look ahead to the lower left panel of Figure 12). Since there was a large gap of five days with the first group of points in this interval these 12 points were removed (in the lower left panel of Figure 12 these are marked by the bracket). This ensured a simpler fitting function with fewer inflections. (In Section 3.4 we will see that this one time chunk can have a large influence on the amplitude spectrum.) The time interval RJD = 54,935–54,955 was thus divided into two chunks. In the first (chunk 8) a second order polynomial was used to fit the trend, and a third order polynomial for the next chunk (chunk 9).

The interval RJD = 55,672–55,692 had good sampling and the data looked periodic. Two sine functions with ν = 0.043 day⁻¹, K = 1.2 m s⁻¹ and ν = 0.08 day⁻¹, K = 2.3 m s⁻¹ values found by pre-whitening the data were used to fit the activity variations. These data composed chunk 20.

Table 5 lists the Julian day of the time chunks, the time span, ΔT, the number of planet orbits during this span, N_orb, the number of data points used in the fit, N_data, the fitting function employed (constant = average value subtracted, linear, second or third order polynomial), and the rms scatter in the chunk, σ, after removing the trend but with the planet signal present. Figures 6–8 show the individual time chunks used in the analysis of the HARPS data. The error bars represent the best-case error of 0.8 m s⁻¹. The solid lines represent the fit to the underlying trend in the chunk. As a comparison, the fit to the activity using the sine parameters found by pre-whitening the full data set (Table 2, but without the planet contribution) is also shown. Note that there are many instances where the LTF provides a much better fit to the underlying activity variations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The residual RVs after removing the underlying trend were then combined. These had a standard deviation of σ = 0.94 m s⁻¹. The Scargle periodogram (top panel Figure 9) shows no significant peak anywhere in the frequency range ν = 0–0.5 day⁻¹. The peak at the planet orbital frequency is weak and has an FAP of ≈0.4 as determined by a bootstrap. As a quick test the Fourier component at 3.24 days (ν = 0.309 day⁻¹, f₉ in Table 2) was removed from the full data set and the orbit of α Cen Bb of D2012 added back to the data. The data with the simulated planet were divided into chunks and trend filtered as the previous data, and calculating new trends for the data. The Scargle periodogram of the residuals (lower panel Figure 9) shows that the planet signal should have been easily detected and with an FAP ≈ 0.01%. The local trend filtering (LTF) method does not confirm the presence of a planetary signal at 3.235 days around α Cen B.

Zoom In Zoom Out Reset image size

The LTF procedure was tested further to see how well it could recover known signals in simulated data generated in different ways. This was done on three different simulated data sets.

• 1.  
  The 3.24 day planet period was removed from the RV data using the sine function parameters found by the pre-whitening process. A sine function with the same amplitude as the planet (0.5 m s⁻¹) but with a slightly different period, P = 3.27 days, was inserted back into the data. A different period was employed simply to avoid the frequency in the amplitude spectrum where signal was removed. This simulation keeps the original noise characteristics of the data.
• 2.  
  A signal with a period of 3.37 days and a K-amplitude of 0.5 m s⁻¹ was inserted into the RV data. A different period to that of the planet was used because this possible signal is still in the data and we want to avoid interference between the two. (As far as LTF is concerned a 3.37 day period is no different from a 3.23 day period.) In this simulation all the noise characteristics of the data are kept, and with the assumption that the signal at 3.24 days is also due to noise.
• 3.  
  A simulated activity signal was generated using the first eight frequencies, amplitudes, and phases listed in Table 2. (Hereafter we will refer to this simulated activity signal as "the activity function.") This was then sampled in the same way as the data, and random noise with standard deviation σ = 0.8 m s⁻¹ was added. The orbit of D2012 was then inserted into the data. The advantage of this simulation is that the underlying trend for each chunk may be slightly different from the cases above. The disadvantage is that the noise characteristics, which are now Gaussian, may be different than for the real data.

The LTF technique was applied to each of these three artificial data sets. Although the time span for each chunk was kept fixed, the fit to RV data was repeated in each case. Scargle periodograms of the trend-removed residuals (Figure 10) show that the local filtering process was able to recover the input signal and at a high significance in all cases. A bootstrap analysis with 200,000 shuffles showed no instance where the periodogram of random data exceeded the real periodogram (FAP < 5 × 10⁻⁶). Note that these signals were detected at a much higher level of significance than the original data. Tests using planet orbital variations with a different phase produced similar results. It appears that the filtering process is not suppressing a possible real signal from a planet.

Zoom In Zoom Out Reset image size

The last simulation assumed our best-case estimate of the noise. It is of interest to explore the detection limits of the LTF procedure as a function of different noise levels. Again synthetic data consisting of the activity function plus the planet orbit parameters of D2012 were used, but this time with different levels of random noise added. Figure 11 shows the Scargle power of the RV residuals as a function of the standard deviation, σ, of the random noise. Shown are simulations for three values of the K-amplitude (0.3, 0.4, and 0.5 m s⁻¹). If the real measurement error is close to the best-case value, a K-amplitude of ≈0.35 m s⁻¹ could be detected with an FAP = 1%. For a K-amplitude of 0.5 m s⁻¹ the planet could be detected at the 1% level even for σ as high as ≈1.4 m s⁻¹. For planets with certain orbital periods it should be possible to overcome the activity noise.

Zoom In Zoom Out Reset image size

The LTF method was also used with different time windows and fitting functions to explore, as in the case of pre-whitening, how different filters could influence the results. A second trend filtered version of the data (hereafter LTF2) was made with the following minor differences compared to the previous version of the trend filtering (hereafter LTF1).

• 1.  
  In LTF1 the data taken during RJD = 54,549–54,572 were divided into two chunks because of a 5 day gap. A linear fit to each was made separately (chunks 2 and 3). In LTF2 a parabola was fit to all the data (upper right panel of Figure 12) across the gap.
• 2.  
  In chunk 6 of LTF1 a parabola was fit to the data, but after removing the last data points that came after a 2 day gap. In LTF2 these points were included, but it forced one to use a higher order polynomial (top right panel of Figure 12).
• 3.  
  For measurements made during RJD = 54,933–54,956 the data were divided into two chunks (chunks 8 and 9) in LFT1 because of a five day gap in the time sampling. In LTF2 data chunks 8 and 9 were combined and a multi-component sine function was fit to the underlying trend throughout the time interval (lower left panel of Figure 12).
• 4.  
  For chunk 20 an additional sine component was used as determined from pre-whitening the data (lower right panel of Figure 12).

Zoom In Zoom Out Reset image size

Figure 12 shows the new trend fits to these time chunks. The total rms scatter of the data (including the other chunks) after removing the trends is about 1.00 m s⁻¹, or slightly worse than for LTF1. For all the other chunks the same time windows and fitting functions were employed as per LTF1. The periodogram of the combined trend-removed residuals shows modest power (FAP ≈ 0.006) at the frequency near the orbital frequency of the planet (top panel of Figure 13). (Also, the highest peak occurs at a slightly different frequency, ν = 0.30615 day⁻¹, or P = 3.26 days). However, it is highly suspicious that a markedly different periodogram is obtained when altering slightly the fit to a small subset of the data. This simulation only reinforces what was found in performing the pre-whitening on the full and epoch data sets—slight differences in the way the RV data are filtered can produce dramatically different results.

Zoom In Zoom Out Reset image size

In the course of filtering the activity signal from the individual chunks we discovered that data in the time window JD−2,400,000 = 54,933–54,956 (hereafter referred to as "Chunk8-9" since it is a combination of chunks 8 and 9 in Table 5) had peculiar frequency characteristics that may influence the outcome when using different approaches to filtering the activity signal. A Fourier analysis of this chunk revealed a significant signal at 3.3 days. Pre-whitening of the data reveals two additional frequencies (Table 6) with frequencies near the rotational frequency, ν_rot, and ≈5ν_rot. The 3.3 day signal is of particular interest since this is uncomfortably close to the planet period. However, it is unlikely that this is due entirely to the planet since its amplitude is too large by a factor of two. This signal is significant as a bootstrap analysis shows that FAP = 3.5 × 10⁻⁴. In producing the residuals from Chunk8-9 that were used in LTF1 the 3.3 day period was kept in the data for the obvious reasons that it nearly coincides with the planetary signal.

As a test we tried different filtering of the data only in Chunk8-9. The first and last half of the data in the chunk were filtered in different ways (e.g., second, third order polynomials, or sine functions). Different fits were performed after deleting points after a large time gap, first and last points from the data, data showing large variations with respect to adjacent values, etc. Ten different versions of filtering the chunk were tried, and in all cases no more than 20 points were removed. The residuals were then added to the rest of LTF2 residuals. The Scargle power of the total residuals from all chunks ranged from as low as z = 7.0 to as high as z = 11.4. The average power was 8.8 ± 1.3. This corresponds to a range in FAP of 3%–30%. We also took the original LTF2 residuals from Chunk8-9 and replaced the corresponding time values in LTF1, and this alone boosted the Scargle power at the planet orbital frequency from z = 6.4 (FAP = 0.4) to z = 9.75 (FAP = 0.05).

In most cases the highest power in the periodogram in the frequency range 0.25 < ν < 0.35 day⁻¹ was not at the planet orbital frequency, but rather P = 2.94 days (ν = 0.3397 day⁻¹). The lower panel of Figure 13 shows one such filtered version of Chunk8-9 where the dominant peaks are not coincident with the planet orbital frequency.

We checked whether the planet signal could be detected in the subset RV measurements without data from Chunk8-9. Local trend filtering (LTF1) clearly shows no significant signal in these subset data (top left Figure 14). However, a 3.15 day periodic signal (K = 0.5 m s⁻¹) that was inserted into the data was found after applying LTF1 even without Chunk8-9. The same results were found when pre-whitening the RV data without Chunk8-9. A weak but insignificant peak is found at the planet orbital frequency (top right panel Figure 14). Applying the pre-whitening to the data with the artificial planet (P = 3.15 days) can recover the input signal at a much higher significance (lower right panel in Figure 14). In summary, both pre-whitening and LTF methods should have been able to detect the planetary signal of α Cen Bb even without the data from Chunk8-9, but they do not.

Zoom In Zoom Out Reset image size

The behavior of the Scargle power as a function of the number of data points is a good way to assess the significance of a real periodic signal (see Hatzes & Mkrtichian 2004). A real signal should have Scargle power that increases in an expected way as you add more data. The behavior of the statistical significance for the complete data was also inconsistent with the expectations of a real signal. Figure 15 shows the power at the orbital frequency of the planet as a function of the number of data points, N, using the residuals from local trend fitting and adding data sequentially in chronological order. We show the results for LTF1 and LTF2.

Zoom In Zoom Out Reset image size

The figure also shows three simulated data sets. For these we added the orbit of α Cen Bb to the activity function generated from the sine components found by the pre-whitening of the individual epochs. The total RV curve was then sampled in the same way as the data, and three different levels of random noise were added (σ = 0.8, 1.0, and 1.4 m s⁻¹). Local trend fitting was then applied using the same time windows as LTF1, but fitting the trends separately for each choice of random noise.

There are several features to note about this figure. The slope of the power versus N function is much steeper for the simulated data, even with σ as high as 1.4 m s⁻¹. LTF1 shows power that is essentially flat except for the slight up-tick after the last data points are added. Even then the FAP is ≈40%. The power from LTF2 behaves more erratically. There is a sharp increase as the first data points are used, followed by just as sharp a decline as more data are added. The "high" significance of the planet detection in LTF2 only occurs after adding the last 100 data points, which is inconsistent with what one expects for random noise. This argues that the LTF1 choice of filtering may be a better approach to filtering out the activity variations.

We conclude that even though LTF2 shows modest power at the orbital frequency of the planet, this is not a significant detection and is consistent with the non-detection of LTF1.

The question naturally arises: "Why do some approaches to removing the activity signal produce such discrepant results?" The sampling coupled with the noise characteristics may give us some insight into this. In this case it is best to compare the unfiltered data in the Fourier domain. The amplitude spectrum of random noise with σ = 1.2 m s⁻¹, our worst-case estimate of the noise level, that is sampled in the same way as the data shows several peaks near the planet frequency and with amplitudes comparable to the K-amplitude of the planet (top panel of Figure 16). Noise in the presence of the activity function shows an amplitude spectrum (middle panel of Figure 16) similar to that of the real data (lower panel of Figure 16). Spectral leakage from the activity signal into the frequency range ν ≈ 0.3–0.31 day⁻¹ may boost power in a noise peak coincident with the planet orbital frequency. The details as to how this noise peak is filtered may explain why the planet is present in some filtering approaches, but not others.

Zoom In Zoom Out Reset image size

To explore further the influence of noise on the amplitude spectrum, synthetic data consisting of only random noise (no activity signal) with σ = 1.2 m s⁻¹ were generated and sampled in the same way as the real data. A total of 100 random data sets were created using different seed values for the random number generator. The top panel of Figure 17 shows the distribution of the strongest peaks in the period range 2.85 < P < 3.8 days (0.26 < ν < 0.35 day⁻¹). This has a peak at P = 3.22 ± 0.01  days. The average amplitude of the peaks is K = 0.24 ± 0.04. In the unfiltered amplitude spectrum the velocity amplitude at the planet orbital frequency is 0.38 m s⁻¹. The amplitude scales approximately linear with σ, so for a noise level of 2 m s⁻¹ the noise peaks would have an amplitude of ≈0.4 m s⁻¹.

Zoom In Zoom Out Reset image size

The epoch subsets of the data show similar noise characteristics. As an example we only show the Epoch 3 data. The lower panel in Figure 17 shows the noise peaks in the same frequency range for random data (σ = 1.2 m s⁻¹) sampled as the Epoch 3 data. The distribution is similar, but in this case the average peak amplitude is slightly higher at K = 0.33 ± 0.13 m s⁻¹.

As seen from Figure 16 the activity signal may also boost the amplitude of the noise peak. The same simulation was performed including the activity function and random noise with σ = 1.2 m s⁻¹. In this case the noise peaks had a mean amplitude of K = 0.39 ± 0.04, essentially the same value as in the unfiltered amplitude spectrum.