Exoplanet Imitators: A Test of Stellar Activity Behavior in Radial Velocity Signals

We simulate magnetic activity RV signals using GP regression with a quasi-periodic kernel, motivated by the work of Haywood et al. (2014). The quasi-periodic kernel has the form

$$\begin{eqnarray}&&k(t,t^{\prime} )={A}^{2}\exp \left[-\displaystyle \frac{{\left(t-t^{\prime} \right)}^{2}}{2{\tau }^{2}}-\displaystyle \frac{2{\sin }^{2}\left(\tfrac{\pi (t-t^{\prime} )}{{P}_{\mathrm{rot}}}\right)}{{\omega }^{2}}\right],\end{eqnarray} \tag{ 1 }$$

where k(t,t') is the correlation weight between observations taken at times t and t'. A is the mean amplitude of the activity signal, τ is related to the evolution timescale of activity features, P_rot is the stellar rotation period, and ω is related to the average distribution of activity features on the surface of the star. The parameter ω describes the level of high-frequency variation expected within a single stellar rotation. Since the level of high-frequency variation defines the number and spacing of local minima or maxima within the timescale of P_rot, it is physically related to the average distribution of magnetic active regions on the stellar surface. For each unique set of GP hyperparameters (A, τ, P_rot, ω), we simulate 100,000 iterations of stellar RV signals by sampling randomly from the GP prior distribution, with each iteration representing a different phase of the activity signal.

Using Equation (1), we model magnetic activity RVs for targets observed by current RV campaigns, assigning observation times and uncertainties from the real RV data with which we later compare our modeled RVs. We apply a bootstrapping method to real RV uncertainties, using them in a different randomized order with each new iteration of modeled RVs. Section 3 details observation times and RV uncertainties for specific test targets.

The values of A, τ, and P_rot for a given target are adopted from the literature when available or are estimated as follows. We set A equal to the standard deviation of the target's real RV residuals (with signals from confirmed exoplanets removed), minus the median value of reported observational uncertainties. With this value, we test a case in which any remaining spread in the RVs, after the removal of known exoplanet signals, can be attributed to a combination of stellar activity and errors associated with observations. Rajpaul et al. (2016) showed how the removal of signals from known exoplanets can lead to spurious periodic signals in RV data sets. However, only the mean values of our modeled RV data sets (A) depend on real RV residuals. The overall structure of each of our modeled data sets are independent of the value of A, and therefore are insusceptible to the spurious periodic signals mentioned above. We set τ and P_rot according to the best estimates from LC and/or RV analyses performed on the data sets. Section 3 details the values of A, τ, and P_rot adopted for specific targets.

For the purpose of our tests, we use two values of τ to probe different relationships between stellar rotation and the evolution timescale of magnetic active regions. In the first, with τ ≈ P_rot, we utilize real τ estimates from the LC and/or RV analyses mentioned above. In the second, with τ = 10 P_rot, we explore the case of highly stable magnetic activity features, compared to measured activity lifetimes on Sun-like stars (Giles et al. 2017). This is the case of an unchanging stellar surface over multiple rotations and timescales of typical RV observations. Some faculae regions fall under this category, with features surviving up to 10 times as long as spots (Collier Cameron et al. 2019).

To explore whether uncertainties associated with hyperparameter estimates affect final maximum peak distributions, we performed additional simulations using A, τ, and P_rot values falling at the high and low limits of their computed uncertainties. In most of these test cases, final distributions and occurrence rates had similar overall trends to simulations using our published hyperparameter values. Any exceptions to this are further discussed in Section 3.

As mentioned above, ω is physically related to the average distribution of magnetic active regions on the stellar surface. Models have demonstrated that even highly complex activity distributions will average to just two to three large active regions in a given rotation (Jeffers & Keller 2009). The distribution ω = 0.5 ± 0.05 is consistent with this behavior, allowing for two to three local minima or maxima per rotation. This prior on ω has been used successfully to determine exoplanet masses from a number of RV data sets (e.g., Haywood et al. 2014, 2018; Grunblatt et al. 2015; López-Morales et al. 2016). While several RV characterizations have used broader priors on ω, the results from Jeffers & Keller (2009) make a strong case for the much tighter Gaussian prior above (e.g., Faria et al. 2016; Mortier et al. 2016; Astudillo-Defru et al. 2017; Cloutier et al. 2017). Broader priors risk overfitting other noise signals, and mistakenly attributing them as part of the magnetic activity signal. We simulate five different values of ω for each target, sampling evenly from the above distribution, i.e., ω = [0.45, 0.475, 0.5, 0.525, 0.55].

In each simulation, we investigate distributions of maximum RV periodogram peak locations over 100,000 iterations. In each iteration, we generate an RV signal and calculate a generalized Lomb–Scargle (GLS) periodogram on the signal (Lomb 1976; Scargle 1982; Zechmeister & Kürster 2009). We calculate the GLS periodogram with a lower limit of 1.5 days (to avoid the 1 day peak due to nightly observations) and an upper limit of half the baseline of the observations' time span (to consider only periods detectable in the simulated data). From each GLS periodogram, we record the period of the maximum statistically significant peak, with a false alarm probability (FAP) rate >1%. We repeat the process of generating modeled RVs and identifying maximum periodogram peaks over a number of iterations, plotting the final distribution in a histogram.

We consult the final distribution of peak periods to calculate how often maximum peaks in the periodogram occur at a series of important periods, detailed below. We define two occurrence rates at any given period: the first is the percentage of iterations with significant (FAP > 1%) peaks in the RV periodogram that have a maximum peak falling within 5% of that period, and the second is the percentage of all 100,000 iterations that have a maximum peak falling within 5% of the period. We base the 5% metric on the range of uncertainties produced by LC estimates of P_rot, used as priors in the RV fits (e.g., Buchhave et al. 2016; López-Morales et al. 2016; Dai et al. 2017).

We calculate occurrence rates at periods related to P_rot, including P_rot itself and integer multiples of P_rot up to the longest period for which the periodograms were calculated. These periods also include rotational harmonics (e.g., P_rot/2, P_rot/3). We calculate occurrence rates for the same number of rotational harmonics as calculated for integer multiples of P_rot. For example, if we calculate occurrence rates for integer multiples up to P_rot × 5, we calculate occurrence rates for rotational harmonics down to P_rot/5. To track the window function signal, we also calculate occurrence rates for the period of the cadence peak, the maximum peak produced by the cadence of observations. To calculate the cadence periodogram, we produce a signal with the same time stamps as the observations and replace RV amplitudes by random values from the uniform distribution 1.0 ± 1 × 10⁻¹⁵. We then calculate a GLS periodogram on the cadence signal with the same period limits used to calculate the simulated RV periodograms.

Finally, we calculate occurrence rates for maximum peaks falling at a specified period of interest (POI), typically the period of an exoplanet candidate in question. Section 3 details POI selections for specific targets. For targets with multiple seasons of simulated RVs, we also determine a rate of time coherence at the POI. In iterations with a maximum peak occurring at the POI, we calculate the GLS periodogram of each independent season by setting unused RVs to zero with an error of 100 m s⁻¹. This method preserves periodic signals inherent to the observational cadence and was first described in Dumusque et al. (2012). We define a maximum periodogram peak at the POI to be long-lived if it remains the maximum peak in each of the periodograms of each individual season. We define the rate of time coherence as the number of iterations with a long-lived maximum peak occurring at the POI divided by the total number of iterations with a maximum peak at the POI.

K2-131 is a solar-type star with P_rot = 9.68 days and one ultra-short period transiting exoplanet companion, K2-131b (P_b = 0.369 days), confirmed by Dai et al. (2017) with combined RV observations from HARPS-N (Cosentino et al. 2012) and the Magellan Planet Finder Spectrograph (PFS; Crane et al. 2010). In addition to the RV signal of K2-131b, Dai et al. (2017) found a statistically significant peak at 3.0 days in the periodogram of their combined PFS/HARPS-N RV observations, reproduced here in Figure 1. The peak raised the possibility of a potential additional nontransiting exoplanet in the system. However, they attributed the signal to magnetic activity due to its close proximity to the second rotational harmonic (P_rot/3 = 3.2 days) and additional detection of the signal in data from two magnetic activity indicators.

Zoom In Zoom Out Reset image size

We simulated combined RV signals of K2-131 and its known transiting exoplanet, K2-131b. We set a value of POI = 3.0 days to investigate whether the combined signal from K2-131's magnetic activity and K2-131b could produce the 3.0 day periodic signal discussed above. We used observation times and RV uncertainties from the combined set of 41 HARPS-N and 32 PFS observations of K2-131, in which the 3.0 day signal was originally detected (Figure 2; Dai et al. 2017). All of the observations were taken within a single season, the HARPS-N data between 2017 January and April and the PFS data over six nights in 2017 March and April. The observations have a baseline of approximately 66 days, so we set the upper limit for detection of a periodic signal at one-half that time span, approximately 33 days.

Zoom In Zoom Out Reset image size

We simulated magnetic activity RVs of K2-131 using GP regression and the quasi-periodic kernel described in Equation (1). In their original RV analysis, Dai et al. (2017) utilized GP regression to fit for magnetic activity. We used the final GP hyperparameters reported in Table 7 of their paper: A = 26.0 m s⁻¹, P_rot = 9.68 days, τ = 8.9 days. We test this set of hyperparameters with each of the five average activity distributions, ω = [0.45, 0.475, 0.50, 0.525, 0.55]. We also ran a second set of simulations with the activity evolution timescale increased to τ = 96.8 days, in order to explore the case of activity features that remain stable over the timescale of observations.

As described in Section 2.2, we performed additional simulations using A, τ, and P_rot values falling at the high and low limits of the uncertainties reported in Table 7 of Dai et al. (2017). Most of these test cases yielded final occurrence rates and overall trends similar to those reported in Section 4. However, cases testing the low limit of P_rot have an overlap in values at the POI and P_rot/3 (3.0 days and 3.18 days, respectively), within the 5% error. In these cases, occurrence rates at the POI increased to resemble values reported for P_rot/3 in Tables 1 and 2.

In all cases, we simulated the planetary RVs of K2-131b with a Keplerian signal, assigning the exoplanet parameters reported in Table 7 of Dai et al. (2017): K = 6.55 m s⁻¹, P = 0.369 days, e = 0, t_c = 3582.9360 (BJD—2,454,000), where K is the semi-amplitude of the exoplanet RV signal, P is the orbital period, e is the orbital eccentricity, and t_c is the central time of transit. With our lower period limit for periodogram calculations set to 1.5 days, we did not track maximum peak occurrences at the orbital period of K2-131b (P = 0.369 days) in a subsequent analysis of simulated magnetic activity signals. However, since the signal at 3.0 days appeared in real RV data for K2-131 before the removal of the signal from K2-131b, we included the planetary signal to investigate the possibility of the combined signals from magnetic activity and K2-131b producing a maximum peak at the POI.

We generated a total of 10 simulations, using the two values of τ and five values of ω given above, with A and P_rot fixed. In each simulation, we generated 100,000 iterations of RVs from K2-131 and its known exoplanet companion. Figures 3 and 4 show example histograms for simulations with τ = 8.9 days and τ = 96.8 days, respectively, both with ω = 0.55. All five values of ω produced similar final maximum peak distributions, but for simplicity, we only show distributions for a single value of ω. Results for the other values of ω are reported in Tables 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Tables 1 and 2 list occurrence rates for simulations with τ = 8.9 days and τ = 96.8 days, respectively. We calculated occurrence rates for maximum peaks located within 5% of P_rot, its rotational harmonics (P_rot/2, P_rot/3), its integer multiples (P_rot × 2, P_rot × 3), the cadence peak (6.0 days; Figure 1), and the POI (3.0 days; Figure 1). With only a single season of RVs available in the original data set, we did not calculate a rate of time coherence at the POI.

Kepler-20 is a solar-type star with P_rot = 27.4 days and five confirmed transiting exoplanet companions, all with orbital periods shorter than 80 days (Table 3, Fressin et al. 2012; Gautier et al. 2012). Buchhave et al. (2016) published a sixth nontransiting companion, with combined HARPS-N and HIRES data after removing RV signals from the largest three known transiting planets. The predicted RV amplitudes of signals from the smallest two planets are too small to be detected with HARPS-N and HIRES precision. Buchhave et al. (2016) included the nontransiting planet in their fit after detecting a maximum peak located at 34.94 days in the HARPS-N RV residual periodogram that remained coherent over two seasons of observation separated by approximately 138 days (Figure 5). They concluded the long-lived signal to be planetary, considering the approximately 7 day difference between the signal in question and P_rot of the star. Buchhave et al. (2016) also observed no correlation between their RVs and the three activity indicators tested, and therefore did not include a model for stellar activity in their final fit.

Zoom In Zoom Out Reset image size

Table 3.  Transit and Orbital Parameters for Kepler-20's Five Transiting Exoplanet Companions

Parameter	Kepler-20b	Kepler-20c	Kepler-20d	Kepler-20e	Kepler-20f
Orbital period (days)	${3.696115}_{-.000001}^{+.000001}$	${10.85409}_{-.000003}^{+.000003}$	${77.6113}_{-.0001}^{+.0001}$	${6.098523}_{-.000014}^{+.000006}$	${19.57758}_{-.00012}^{+.00009}$
Tc (BJD—2,454,000)	${967.5020}_{-.0002}^{+.0003}$	${971.6080}_{-.0002}^{+.0002}$	${997.730}_{-.002}^{+.001}$	${968.932}_{-.001}^{+.002}$	${967.5020}_{-.0002}^{+.0003}$
Orbital eccentricity	${0.03}_{-0.03}^{+0.09}$	${0.16}_{-0.09}^{+0.01}$	<0.6*	<0.28**	<0.32 **
Planet radius (R⊕)	${1.868}_{-0.034}^{+0.066}$	${3.047}_{-0.056}^{+0.084}$	${2.744}_{-0.055}^{+0.073}$	${0.865}_{-0.028}^{+0.026}$	${1.003}_{-0.089}^{+0.050}$
Planet mass (M⊕)	${9.7}_{-1.44}^{+1.41}$	${12.75}_{-2.24}^{+2.17}$	${10.07}_{-3.70}^{+3.97}$	...	${10.07}_{-3.70}^{+3.97}$

Note. The majority of values are taken from results reported in Table 4 of Buchhave et al. (2016). Values marked by a single asterisk are from the fit reported in Table 2 of Gautier et al. (2012). Values marked by a double asterisk are from the fit reported in Table 1 of Fressin et al. (2012).

Download table as:  ASCIITypeset image

We simulated RVs of Kepler-20's magnetic activity signal. We set a value of POI = 34.94 days to investigate whether the stellar activity signal alone could produce the long-lived periodic signal that Buchhave et al. (2016) attributed to a nontransiting planet, as discussed above. We used observation times and RV uncertainties from the 104 HARPS-N observations of Kepler-20, collected over two seasons between 2014 and 2015 (Figure 6). We did not include the 30 available HIRES observations because the peak originally motivating the planetary signal at 34.94 days was detected in the periodogram of HARPS-N RVs only. The signal was not detected in the periodogram of combined HARPS-N and HIRES RVs. The HARPS-N observations have a baseline of approximately 570 days, so we set the upper limit for the detection of a periodic signal at one-half that time span, approximately 285 days.

Zoom In Zoom Out Reset image size

We simulated magnetic activity RVs of Kepler-20 using GP regression and the quasi-periodic kernel described in Equation (1). Unlike in the case of K2-131, the original RV analysis performed by Buchhave et al. (2016) did not utilize GP regression to fit for magnetic activity, and therefore does not provide estimates of A, τ, and P_rot. We estimated those three hyperparameters as described below.

We set A equal to the standard deviation of the HARPS-N RV residuals minus the median value of the associated uncertainties. We removed signals from Kepler-20b, Kepler-20 c, and Kepler-20d to calculate the RV residuals, because it was in these residuals that Buchhave et al. (2016) detected the signal at 34.94 days. As explained before, Kepler-20e and Kepler-20f induce signals smaller in amplitude than HARPS-N RV precision, and therefore cannot be reliably fitted and removed (Buchhave et al. 2016). To estimate τ and P_rot, we performed an ACF analysis on the same Kepler-20 photometric data originally analyzed by Buchhave et al. (2016), consisting of LCs from 15 Kepler campaigns (Q3-Q17) collected between 2009 and 2013. We applied discrete shifts to the LC and cross-correlated the shifted LCs with the original, revealing peaks separated by a timescale related to P_rot, with correlation powers dropping off at a rate related to τ (McQuillan et al. 2014; Giles et al. 2017). Our ACF analysis failed to converge on a final value for τ, but did produce a P_rot estimate of 27.4 ± 0.8 days. We instead used a value of τ = 22.9 ± 0.2 days, obtained from the relationship between τ, stellar effective temperature, and the scatter in the photometric LC, described in Equation (8) of Giles et al. (2017). We used the following final hyperparameter values in our reported simulations of Kepler-20: A = 2.31 m s⁻¹, τ = 22.9 days, P_rot = 27.4 days.

Again, using the method described in Section 2.2, we found that uncertainties associated with our τ and P_rot estimates did not affect final maximum peak distributions. We also tested two additional values of A: the standard deviation of HARPS-N RV residuals without the mean HARPS-N observational uncertainty subtracted (6.07 m s⁻¹) and the reported semi-amplitude of the RV signal attributed to Kepler-20g (4.10 m s⁻¹; Buchhave et al. 2016). In all of these test cases, final distributions and occurrence rates demonstrated similar overall trends to those reported in Section 4.

We ran a second set of simulations with the activity evolution timescale set to τ = 274.0 days, in order to explore the case of activity features that remain stable over the timescale of observations. For both sets of A, τ, and P_rot values, we again tested five average activity distributions, ω = [0.45, 0.475, 0.50, 0.525, 0.55].

We generated a total of 10 simulations, using the two values of τ and five values of ω given above, again with A and P_rot fixed. In each simulation, we generated 100,000 iterations of magnetic activity RVs for Kepler-20. We plotted a histogram of the distribution of maximum peak periods over all iterations. Figures 7 and 8 show example histograms for simulations with τ = 22.9 days and τ = 274.0 days, respectively, both with ω = 0.55. All five values of ω produce similar final maximum peak distributions, as shown in Tables 4 and 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Tables 4 and 5 list occurrence rates for simulations with τ = 22.9 days and τ = 274.0 days, respectively. We calculated occurrence rates for maximum peaks located within 5% of P_rot, its integer multiples (P_rot × 2, ..., P_rot × 10), its rotational harmonics (P_rot/2, ..., P_rot/10), the cadence peak (37.3 days; Figure 5), and the POI (34.94 days; Figure 5).

In iterations with a maximum peak at the POI, we checked whether the maximum peak was long-lived, surviving both seasons of observation. We used the method described in the final paragraph of Section 2.3 to calculate a rate of time coherence at the POI, listed in Tables 4 and 5. Figure 9 shows an example periodogram in which simulated Kepler-20 activity RVs, with ω = 0.55, produce a long-lived peak at the POI.

Zoom In Zoom Out Reset image size

Figures 3 and 7 show final distributions for simulated RVs of K2-131 and Kepler-20 with evolving magnetic active regions, and therefore stellar surfaces that change over the timescale of observations (τ ≈ P_rot). In these cases, a large fraction of simulated active region signals fail to produce a maximum peak at P_rot or a related period. Some simulations fail to produce significant peaks at all. This is rarely the case for K2-131, with only 0.1%–0.2% of iterations lacking significant peaks, but in the case of Kepler-20, a whopping 81.0%–83.1% of simulated signals fail to produce significant peaks. The range in reported rates is a result of varying ω values, while for simplicity, figures are only shown for a single average active region distribution (ω = 0.55). A large fraction of the statistically significant maximum peaks that do occur in the RV periodogram are located at periods unrelated to P_rot: 80.6%–81.0% of maximum peaks for K2-131 and 71.5%–73.0% of maximum peaks in the case of Kepler-20 (Tables 1 and 4). If analyses of real RV data disregard stellar active regions or use inadequate models, these spurious periodic signals could interact with Keplerian signals of known exoplanets and lead to inaccurate RV mass measurements with underestimated errors.

Figures 4 and 8 show final distributions for simulated RVs of K2-131 and Kepler-20 with magnetic active regions, and therefore stellar surfaces, that remain unchanged over the timescale of observations (activity evolution timescale increased to $\tau \gg {P}_{\mathrm{rot}}$). In these cases, where strong rotation signals would typically be expected, a large fraction of simulated activity signals fail to produce a maximum peak at P_rot or a related period. While essentially all iterations of simulated K2-131 signals produced significant peaks, simulated Kepler-20 RVs still fail to produce significant peaks at any period in 71.6%–72.4% of iterations. Again, a considerable fraction of the significant maximum peaks that do occur in the RV periodogram are located at periods unrelated to P_rot: 34.3%–38.2% of peaks in the case of K2-131 and 14.0%–15.0% of peaks for Kepler-20 (Tables 2 and 5). These results demonstrate that spurious peaks are inherent to RVs of many activity distributions, even when magnetic surface features are unchanged over the timescale of observations. Therefore, even with high-cadence observations of exoplanets orbiting relatively inactive stars, spurious periodic signals could still lead to aforementioned errors in RV mass determinations.

All 10 of our final distributions show maximum peaks in the simulated RV periodograms favoring a period located between P_rot and the cadence peak. In cases where τ ≈ P_rot (Figures 3 and 7), this feature in the histogram could be attributed to random variations in the overall distributions. However, cases where τ ≫ P_rot (Figures 4 and 8) show clear features in the histogram between P_rot and the cadence peak. This suggests that maximum signals in RV periodograms have a tendency to occur at a period related to the limited sampling of the signal at P_rot. In the case of Kepler-20, this favored period occurs at the POI, 34.94 days.

For simulations of K2-131 using a real estimate of the activity evolution timescale (τ ≈ P_rot), 1.0%–1.7% of significant maximum peaks in the RV periodogram occurred at the POI, 3.0 days (Table 1). Considering the close proximity of the POI to the second rotational harmonic (P_rot/3 = 3.2 days), these occurrence rates are strikingly low relative to occurrence rates at other period values. However, as mentioned in Section 3.1, simulations including the lowest value of P_rot within its reported uncertainty lead to an overlap between the POI and P_rot/3, within the allowed 5% window. In this case, occurrence rates at the POI increase to 2.5%–3.6%. Due to the large range of potential maximum peak values in a given distribution, maximum peaks are not particularly likely to occur at any given period value. Therefore, while occurrence rates at the POI seem low, they appear more significant when compared to the highest occurrence rates at P_rot, 6.6%–7.7% (Table 1).

Simulations of Kepler-20 using a real estimate of the activity evolution timescale (τ ≈ P_rot) produce a maximum peak in the RV periodogram at the POI in 4.3%–4.8% of iterations. These rates are relatively high when compared to the highest occurrence rates at P_rot, 6.4%–6.9% (Table 4). Occurrence rates at the POI in simulations for Kepler-20 also seem more significant when compared with occurrence rates at the POI in similar simulations for K2-131. Since the POI for K2-131 is relatively close (0.2 days away) to its nearest P_rot relative (P_rot/3), and the POI for Kepler-20 is more than seven days away from P_rot, we would expect to see greater occurrence rates at the POI in the case of K2-131. However, simulations of K2-131 and K2-131b produce much lower relative occurrence rates at the POI, 3.0 days. Our results therefore prove contrary to the assumption that maximum peaks from magnetic activity in RV periodograms usually occur at periods related to P_rot.

Simulated Kepler-20 RVs further defy assumptions about magnetic activity signals with 16.5%–18.7% of maximum peaks occurring at 34.94 days being long-lived, remaining the maximum peak in periodograms of both seasons of observation. These long-lived periodic signals that occur many days from a star's estimated P_rot could be misinterpreted as exoplanet signals, particularly when analyses disregard models for stellar activity. While our results alone cannot rule out or confirm the existence of nontransiting planets around any target for certain, simulations of Kepler-20 demonstrate how spurious signals in RV periodograms from magnetic active regions could appear planetary in nature.