© ESO, 2016

1. Introduction

Radial velocity (RV) measurements have demonstrated great successful at discovering extrasolar planets. New instrumentation, improved calibration methods, and innovative analysis techniques have steadily improved the RV precision to the point that we can routinely make measurements with precisions of ~1 m s^-1 or better. At this level of precision the stellar intrinsic noise now represents a significant contribution to the measurement “error”, often referred to as the stellar RV “jitter”. In the best case this RV jitter can hinder the detection of planetary companions, in the worse case it can create a periodic signal that is interpreted as arising from a planetary companion.

M dwarf stars are objects that have become particularly attractive for Doppler surveys because the RV amplitude of the host star caused by an orbiting terrestrial planet in the habitable zone is ~few m s^-1, a value easily measured by current techniques. Unfortunately, M dwarfs can be active and the orbital period of a planet in the habitable zone is days to weeks and this is comparable to the time scales of stellar activity (rotational modulation, spot evolution, etc.). Activity-related RV jitter may produce false planets.

A case in point is the planetary system around the M dwarf star GL 581. Mayor et al. (2009) reported four planet in the system and shortly afterwards Vogt et al. (2010) claimed six planet candidates with periods up to 433 d. Two of these planets, GL 581d and GL 581g, were of particular interest because they had orbital distances that placed them within the habitable zone of the star. Subsequent RV measurements (Forveille et al. 2011) could not confirm the presence of GL 581g, although this has been the subject of debate (Vogt et al. 2012). Hatzes (2014) demonstrated that the signal attributed to GL 581g was real and significant, but most likely due to activity as it was not coherent on long time scales.

Baluev (2013) was the first to question the reality of GL 581d based on an analysis of the red noise in the RV data. This doubt was considerably strengthened through a study by Robertson et al. (2014). They argued that the 66-d orbital period of GL 581d was actually a harmonic of the 130 d rotation period of the star. This was based on an apparent correlation between the RV variations due to GL 581d and the equivalent width of Hα, their so-called I_Hα index. Correcting the RVs for this correlation eliminated the signals attributed to GL 581d and GL 581g, while boosting the significance of the other three planets.

The nature of the 66-d RV period of GL 581 (hereafter referred to as the “orbital period” regardless of whether the planet exists or not) still remains the subject of debate. Anglada-Escudé & Tuomi (2015) questioned the conclusions of Robertson et al. (2014), arguing that their result came from an improper use of periodograms on the residual data. Although Anglada-Escudé & Tuomi agreed that there was a substantial correlation between the RV variations and I_Hα index, they argued that there was no clear evidence of time variability of this index.

It is important to establish the true nature of the RV variations of the purported planet GL 581d as this has important consequences for RV searches for planets in the habitable zone of M dwarf stars. This is especially true since different analyses of the same I_Hα index measurements arrived at different conclusions. Therefore, we investigated the Hα variations using a different approach to the one taken by Robertson et al. (2014), namely a Fourier analysis of the frequency components in the Hα time series. The aim of this work is to confirm or refute the anti-correlation between RV and Hα equivalent width found by Robertson et al. (2014).

Fig. 1 Time series of the Hα equivalent width measurements for GL 581. The line represents the five frequency fit using the entries in Table 1 (f1−f5).

2. The temporal variations of Hα

2.1. Frequency analysis

For this study we used the I_Hα index measurements of GL 581 from Robertson et al. (2014) which are shown in Fig. 1. The data were analyzed using a Fourier approach. We first found the dominant component in the Fourier transform, and this was fit with a sine function of the appropriate period, amplitude, and phase. The contribution of this frequency was then removed from the time series and we proceeded to the next dominant sine term in the residuals. This process is often referred to as “pre-whitening”. The advantage of this technique is that since you are fitting and removing a sine function sampled in the same way as the data, alias frequencies are also removed.

Normally one associates pre-whitening with finding strictly periodic signals in your data, often in cases associated with stellar oscillations. However, pre-whitening can be more versatile in that it finds the dominant Fourier components that describe the overall variations in your time series, even if they do not appear to be periodic. The mathematical foundation for this is that trigonometric functions form basis set and a linear combination of these provide an alternative mathematical representation of the series. Pre-whitening finds those Fourier components in the time series that are clearly above the noise level. The sum of these may provide an adequate representation of the time variations in the RV due to activity.

Fig. 2 Pre-whitening procedure on the Hα data taken between JD = 2 455 200 and 2 455 645. The original amplitude spectrum before (top) and after (bottom) removing the dominant frequency at 0.0078 d-1. The vertical dashed line marks the orbital frequency of GL 581d (ν = 0.015 d-1, P = 66.7 d).

For our pre-whitening analysis we used the program Period04 (Lenz & Breger 2005). Table 1 lists the dominant frequencies (labeled f₁ – f₅), equivalent periods, and amplitudes found by the Fourier analysis. The pre-whitening procedure was stopped when the dominant peak in the residuals had a peak less than four times the surrounding noise level. Peaks of this amplitude generally have a false alarm probability of ~1% Kuschnig et al. (1997). Once the dominant frequency components were found a simultaneous fit was made to the data optimizing the amplitude, period, and phase of the individual sine terms. The multi-sine fit using these frequencies is shown as the curve in Fig. 1. The first two frequencies, f₁ and f₂ are associated with the ~130 d rotation period of GL 581. The third frequency is the orbital frequency of GL 518d.

The orbital frequency can also be easily detected in a subset of the Hα data taken JD = 2 455 200 and 2 455 645. Because of the shorter data set, fewer Fourier components were found, namely one at ν = 0.0079 d^-1 (P = 126.9 d) and another at the orbital frequency of GL 581d, ν = 0.0150 d^-1 (P = 66.4 d). Figure 2 shows the pre-whitening procedure on the subset I_Hα data. The removal of the dominant peak at ν = 0.0077 d^-1 (P = 129 d) results a strong peak at ν = 0.0152 d^-1, or P = 65.8 d, near the orbital period of GL 581d.

Fig. 3 Top: L-S periodogram of all IHα measurements after removing the contribution of all frequencies (f1, f3–f5 in Table 1) except that associated with GL 581d (f2). bottom: L-S periodogram of the IHα measurements taken between JD = 2 455 200 and 2 455 645 after removing the contribution of f6 in Table 1.

2.2. Statistical significance of the variations

The statistical significance of a periodic signal can be assessed using the Lomb-Scargle (L-S) periodogram (Lomb 1976; Scargle 1982). The top panel in Fig. 3 shows the L-S periodogram of the residual I_Hα measurements after removing all frequencies except that associated with GL 581d (i.e. f₁, f₂, f₄, and f₅ were removed). (Henceforth, we shall refer to “residual” I_Hα measurements as those values where the contribution of all frequencies except f₃ have been removed from the data.)

One can estimate the false alarm probability (FAP) from the L-S power, z, and the expression FAP = 1 − (1−e^{− z})^N ≈ Ne^{− z}, where N is the number of independent frequencies (Scargle 1982). For the entire data set this results in FAP ≈ 10^-13, For the subset data this is a FAP ≈ 10^-6. The latter value was confirmed using a bootstrap procedure (Murdoch et al. 1993; Kürster et al. 1997) and 2 × 10⁵ random shuffles of the data. In no instance was the L-S power of the random data greater than the actual value (FAP <5 × 10⁶).

The L-S periodogram of the residual I_Hα data from the subset data (JD = 2 455 200 – 2 455 645) is shown in the lower panel of Fig. 3. Note the dramatic increase in the L-S power when using the full data set. This indicates a long-lived and reasonably coherent signal. However, one should be wary of determining the FAP of a time series that has been modified. In this case the contributions of four frequencies were removed from the data before estimating the FAP. If one removes the dominant peaks in the periodogram the remaining ones will always look more significant. Therefore we took additional approaches to estimate the FAP.

The L-S periodogram of the original time series before pre-whitening shows L-S power of z≈ 11.9 at the orbital frequency of GL 581d (ν = 0.015 d^-1). The above expression gives us the FAP for random noise producing more power than the actual data over a broad frequency range. However, we are interested in assessing the FAP at a known frequency in the data, i.e. the orbital frequency of GL 581d. In this case the FAP is given by FAP = e^{− z} (Scargle 1982), or the previous expression with only one independent frequency (N = 1). This results in FAP ≈ 10^-5.

The FAP was assessed using a revised version of the bootstrap method. The five frequency components in Table 1 were removed from the full I_Hα index data. The resulting residuals represent our “noise” model for the data. The values from this noise model were then randomly shuffled keeping the times fixed. The sine functions from the frequency analysis were then added back into into the noise data, except for the contribution from GL 581d, f₃ (i.e. only f₁, f₂, f₄, and f₅ were added back in). A L-S periodogram was calculated and the maximum power in the frequency range ν = 0.01−0.02 d^-1 found. We chose this frequency range as we are interested in the probability noise would produce significant power at ν = 0.015 d^-1. With 2 × 10⁵ random shuffles of the data there was no instance when the random L-S periodogram produced power higher than the original data.

We also performed the same procedure on the subset of the I_Hα index data. In this case there are fewer Fourier components in the time series, namely the rotational frequency at ν = 0.0079 d^-1 and the orbital frequency of GL 581d, ν = 0.0150 d^-1. Both of these components were subtracted from the data to produce the noise model for the bootstrap procedure. As before, the orbital frequency of GL 581d (ν = 0.0150 d^-1) was not added back into the noise model prior to computing the L-S periodogram. Again, after 2 × 10⁵ shuffles there was no instance of the random periodogram having higher power than the actual data. A bootstrap analysis of the full and subset data indicates that the FAP for the I_Hα variations at the orbital frequency of GL 581d is <2 × 10^-6.

The FAP can also be estimated directly from the Fourier amplitude spectrum and the height of a peak above the background noise. Kuschnig et al. (1997) using Monte Carlo simulations established a relationship between the peak height above background and the FAP (see Fig. 4 in their paper). In the case of the I_Hα amplitude spectrum the peak at ν = 0.015 d^-1 is 4.3 times above the noise level. This results in FAP ≈ 10^-4. Both the Fourier amplitude spectrum and the L-S periodogram analysis both indicate that the I_Hα variations at the orbital frequency of GL 581d are highly significant.

Fig. 4 Phase-binned averages (Δφ≈ 0.1) residual Hα variations phased to the 125-d period for “GL 581d”. The solid curve represents a sine fit. The dashed curve represents the RV orbital solution for GL 581d.

3. I_Hα versus RV variations

Figure 4 shows the residual I_Hα variations phased to the orbital period of GL 581d determined by Hatzes (2014). For clarity the data has been phase-binned on intervals Δφ≈ 0.1. The error bars represent the standard deviation divided by the square root of the number of points in each bin. The solid curve represents a sine fit to the data. Also shown as a dashed line is the RV orbital solution for GL 581d. The I_Hα – RV variations for GL 581d are anti-correlated thus supporting the conclusions of Robertson et al. (2014). The variations can be fit by a pure sine function. If the RV variations attributed to GL 581d are actually due to activity, then these should mimic a circular Keplerian orbit. There appears to be a slight phase shift of ≈0.1 between the RV and I_Hα data but this is not deemed significant and may be an artifact of the binning process.

Over the interval JD = 2 455 200 and 2 455 645 the sampling of the data (both RV and I_Hα) was excellent and almost three cycles of variations were covered. The residual I_Hα data during this time interval are shown in Fig. 5 along with the RV variations due to GL 581d (i.e. the contribution of all other planets and activity have been removed). Note that in this case the RV and I_Hα variations are exactly 180° out-of-phase with each other.

A sine fit to the RV and I_Hα variations over this interval and allowing the period to be a free parameter resulted in consistent values for the period for the two quantities, namely 71.99 ± 1.51 d and 70.64 ± 1.64 d for the Hα and RV, respectively. These periods are slightly longer by about 1.4σ compared to the orbital period of 66.64 ± 0.08 d derived using the full data set. This may be an indication of possible period variations due to differential rotation which may also be consistent with activity-related variations. Unfortunately, we cannot be sure given the short time span of the data.

We also examined the ratio of the amplitudes of the I_Hα and RV variations (A_Hα/A_RV). This ratio should be constant if the two variations stemmed from the same origin. However, we might see temporal variations if the RV amplitude is constant and due to a planet, while the activity signal from I_Hα varied due to the evolution of the stellar active regions.

Figure 6 shows the ratio A_Hα/A_RV for three different epochs. The I_Hα amplitude was calculated in parts per thousand (PPT). Although there appears to be a slight decrease in the ratio with time, this is not significant. To within the errors the ratio A_Hα/A_RV is constant with time. However, given that the active regions may evolve at different time scales to those on sun-like stars the constancy of A_Hα/A_RV be less informative due to the restricted time base of the measurements.

Fig. 5 Top: variations in IHα over the time span JD-2 450 000 = 5272–5464. The curve represents a sine fit with the orbital period of GL 581d. Bottom: RV variations of the purported planed “d” (all other planet signals removed) over the same time span. The curve represents the orbital solution.

Fig. 6 Ratio of the RV to IHα amplitude (IHα amplitude is in PPT) for three different epochs.

4. Discussion

Using a pre-whitening procedure we were able to isolate the variations of I_Hα in GL 581 at a frequency of ν = 0.015 d^-1, P = 66.7 d which is coincident the RV variations attributed to the planet GL 581d. These variations show sinusoidal variations that are 180° out-of-phase with the “orbital” RV variations due to GL 581d. This confirms the anti-correlation between the RV variations of GL 581d and I_Hα found by Robertson et al. (2014). The RV variations attributed to GL 581d are most likely due to stellar activity. GL 581d represents another case where activity-related RV variations can mimic a Keplerian orbit (see Queloz et al. 2001; Bonfils et al. 2007).

The I_Hα variations appear to be long-lived since they are visible in data spanning almost 7 years. This is certainly a cautionary tale for RV programs searching for planets in the habitable zone around M-dwarf stars.

Of course, the presence of the period of GL 581d in I_Hα is no proof that the planet does not exist. There is no physical reason why a planet cannot have the same orbital period as the stellar rotation, or one of its harmonics. There are two arguments against keeping GL 581d as a confirmed exoplanet. First, when one finds significant periodic variations in an activity indicator (photometry, I_Hα, Ca II, bisectors, etc.) with the same period as the RV variations, then this is generally accepted as a “non-confirmation” of the planet candidate. There is no reason to change this criterion simply to “save” a habitable planet. It is better to exclude a questionable exoplanet in the census rather than keeping it when it actually is not there.

Second, Robertson et al. (2014) demonstrated that after correcting the RVs due to their correlation with I_Hα this boosted the significance of the other RV signals that are certainly due to planetary companions. This behavior is also consistent with activity-related RV variations for the 66.7 d signal.

We examined the the ratio of the RV to I_Hα amplitudes. If the RV amplitude was directly tied to the activity then this ratio should remain roughly constant. On the other hand, if the 66.7 d RV variations were due to a planet then the amplitude of these would remain constant, whereas an activity cycle may produce amplitude variations in the I_Hα index. This should result in a variation in the RV to I_Hα amplitude and the planet candidacy of GL 581d could remain. Although our measurements of the RV to I_Hα amplitude are consistent with no variations, the measurement error is too large to exclude with certainty any temporal variations. Unfortunately, to detect any temporal variations in either the I_Hα or RV amplitude would require an inordinate amount of additional measurements. At the present time there is no compelling reason to think that GL 581 is more than a 3-planet system.

References

All Tables

All Figures

	Fig. 1 Time series of the Hα equivalent width measurements for GL 581. The line represents the five frequency fit using the entries in Table 1 (f1−f5).
In the text

	Fig. 2 Pre-whitening procedure on the Hα data taken between JD = 2 455 200 and 2 455 645. The original amplitude spectrum before (top) and after (bottom) removing the dominant frequency at 0.0078 d-1. The vertical dashed line marks the orbital frequency of GL 581d (ν = 0.015 d-1, P = 66.7 d).
In the text

	Fig. 3 Top: L-S periodogram of all IHα measurements after removing the contribution of all frequencies (f1, f3–f5 in Table 1) except that associated with GL 581d (f2). bottom: L-S periodogram of the IHα measurements taken between JD = 2 455 200 and 2 455 645 after removing the contribution of f6 in Table 1.
In the text

	Fig. 4 Phase-binned averages (Δφ≈ 0.1) residual Hα variations phased to the 125-d period for “GL 581d”. The solid curve represents a sine fit. The dashed curve represents the RV orbital solution for GL 581d.
In the text

	Fig. 5 Top: variations in IHα over the time span JD-2 450 000 = 5272–5464. The curve represents a sine fit with the orbital period of GL 581d. Bottom: RV variations of the purported planed “d” (all other planet signals removed) over the same time span. The curve represents the orbital solution.
In the text

	Fig. 6 Ratio of the RV to IHα amplitude (IHα amplitude is in PPT) for three different epochs.
In the text