The Test Case of HD 26965: Difficulties Disentangling Weak Doppler Signals from Stellar Activity^*

Planets orbiting the nearest stars to the Sun are the most highly prized of all exoplanets, as they represent the most accessible targets for follow-up characterization studies. The measurement of precision radial velocities has allowed us to begin to build up a collection of planets orbiting the nearest stars, while also characterizing their orbital parameters. In particular, discoveries like 51 Peg b (Mayor & Queloz 1995), 47 UMa b (Butler & Marcy 1996), 70 Vir b (Marcy & Butler 1996), HD 143361 b and HD 154672 b (Jenkins et al. 2009), HD 86226 b, HD 164604 b, HD 175167 b (Arriagada et al. 2010), HD 128356 b, HD 154672 b and HD 224538 b (Jenkins et al. 2017), GJ 876 b, c, d, e (Rivera et al. 2010), and υ And b, c, d (Wright et al. 2009; Curiel et al. 2011), among others, have allowed us to explore the wide diversity of gas giant planetary systems.

In the last few years, the advances in radial velocity precision that have been driven by technology improvements and better analysis methods have allowed the discovery of the first batch of low-mass planets orbiting nearby stars, e.g., GJ 876 d (Rivera et al. 2010), HD 40307 b, c, d, e, f and g (Mayor et al. 2009; Tuomi et al. 2013a), GJ 581 d (Vogt et al. 2010), GJ 667C b, c and d (Anglada-Escudé et al. 2012, 2013), the candidates orbiting τ Ceti, planets b, c, d, e and f (Tuomi et al. 2013b; Feng et al. 2017), and more recently Proxima Centauri b (Anglada-Escudé et al. 2016a) represent a new population of super-Earth planets not witnessed in the solar system, and are defined as being small planets with masses ∼2–10 M_⊕ that can either be primarily rocky objects or more fluffy, atmosphere dominated worlds (Valencia et al. 2007; Kaltenegger et al. 2011).

In comparison to the gas giants, super-Earths seem to have some dramatically different characteristics, likely related to their formation and early evolution. They generally appear to be orbiting on mostly circular orbits (Tuomi & Anglada-Escudé 2013), come in tightly packed planetary systems (Latham et al. 2011; Lissauer et al. 2011), and do not seem to follow the same metallicity bias as the gas giants (Buchhave et al. 2012; Courcol et al. 2016). In fact, there may be a lack of low-mass planets orbiting nearby and super metal-rich Sun-like stars (Jenkins et al. 2009, 2013). Models that invoke core accretion as the dominant planet formation scenario predict some of these trends, with mass functions rising heavily toward the lowest masses (Mordasini et al. 2008), also shown by analysis of the radial velocity sample of detected planets (Lopez & Jenkins 2012). Planetary formation models also predict a damping of the metallicity bias in planet fraction for low-mass objects, as the stellar metallicity is an observational proxy of the dust content in the inner disk when the planets were undergoing formation. However, the picture may be less clear, as Mulders et al. (2016) have shown that there might be an increase in the occurrence of small rocky planets around host stars with super-solar metallicities and orbital periods <10 days.

Although the radial velocity method has been very successful at planet detection, it is an indirect method and therefore care must be taken when trying to confirm any signal with an amplitude at the few m s⁻¹ level (like many super-Earth signals), as this is the domain where stellar activity effects that are correlated with the rotation of the star can impact the data (Boisse et al. 2011, 2012; Vanderburg et al. 2016). In numerous cases, both large and small planet candidates have been challenged as being due to the effects of stellar activity (e.g., HD 166435 Queloz et al. 2001; HIP 13044, Jones & Jenkins 2014; HD 41248, Santos et al. 2014; GJ 581 d, Robertson et al. 2014; Kapteyn b Robertson et al. 2015; α Cen B b, Rajpaul et al. 2016), with most of these challenges leading to counter-claims (e.g., HD 41248, Jenkins & Tuomi 2014; GJ 581 d, Anglada-Escudé & Tuomi 2015; Kapteyn b and c, Anglada-Escudé et al. 2016b). Therefore, due care must be taken to ensure any signal has been well inspected for the effects of stellar activity and/or stellar rotation.

Once a planet has been confirmed orbiting a nearby star, there exists the ability to perform detailed secondary follow-up studies, like measuring accurate stellar atomic and molecular abundances (Schuler et al. 2015; Meléndez et al. 2017) that could be sign-posts of planetary systems, or searching for transits and secondary eclipse measurements (Baskin et al. 2013; Chen et al. 2014; von Paris et al. 2016). The combination of mass from the velocities and radius from any detected transit allows the bulk density of the planet to be measured (e.g., BD+20594 b, Espinoza et al. 2016; GJ 1214 b, Valencia et al. 2013; GJ 436 b von Braun et al. 2012; Lanotte et al. 2014; 55 Cnc e, de Mooij et al. 2014; Winn et al. 2011) and from there, model comparisons can be made to infer the bulk composition. Therefore, gaining a better understanding of the population of low-mass planets requires the detection of more of these worlds orbiting bright stars in the solar neighborhood.

Here, we present data from a 16-year precision radial velocity monitoring campaign, using multiple high-resolution optical spectrographs, of the nearby (∼5 pc) K0.5 star HD 26965.

HD 26965 (HIP 19849, GJ 166A) is classified as a K0.5V star (Gray et al. 2006) with a visual magnitude of V = 4.43 and an optical color of B − V = 0.82. An activity index of $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }\,=-4.99$ is reported by Jenkins et al. (2011). This value is also consistent with measurements found in other sources in the literature (e.g., −5.09, Gray et al. 2006; −4.97, Murgas et al. 2013) and a comparison with the Sun's mean activity value of $\mathrm{log}{R}_{\mathrm{HK}\odot }^{{\prime} }=-4.91$ (Mamajek & Hillenbrand 2008) tells us that HD 26965 is a chromospherically quiet star.

The remaining stellar parameters were estimated using the Spectroscopic Parameters and atmospheEric ChemIstriEs of Stars code (SPECIES; Soto & Jenkins, submitted). SPECIES derives the effective temperature, surface gravity, metallicity, microturbulence, macroturbulence and rotational velocity, mass, age, and chemical abundances for 11 elements in a self-consistent and automatic manner. The parameters were derived using high-resolution, high signal-to-noise spectra as input for the code, where in this particular case, we have used spectra from HARPS to derive the stellar parameters with SPECIES. The first four parameters were found by measuring the equivalent widths (EWs) for a set of iron lines using the ARES code (Sousa et al. 2007). These values, along with a stellar atmosphere model (Kurucz 1993), were then input to MOOG (Sneden 1973), which solves the radiative transfer equation by imposing excitation and ionization equilibrium. Following on from this, we then derived the chemical abundances, measuring the EWs of a set of lines for each element. Macroturbulence and rotation velocity were computed by measuring the broadening of spectral lines using a Fourier analysis. Finally, mass and age were found by fitting isochrones (Dotter et al. 2008) using the luminosity and temperature of the star. Figure 1 shows the final distributions for the stellar mass, age, and log g using a Markov chain Monte Carlo (MCMC) method within SPECIES. More details about this code, in particular on the treatment of correlations between parameters and uncertainties, can be found in Soto & Jenkins (submitted).

Zoom In Zoom Out Reset image size

We found HD 26965 to have a metallicity [Fe/H] of −0.29 ± 0.13 dex, consistent with previously reported values (e.g., −0.28 dex Valenti & Fischer 2005; Gray et al. 2006; Turnbull 2015), and significantly poorer in metals than the Sun. SPECIES finds a rotational velocity of v sin i = 1.23 ± 0.28 km s⁻¹, which is in agreement with the values of 1.1 ± 1.0, 1.4 ± 0.8, 1.6 ± 0.8 km s⁻¹ reported by Glebocki & Gnacinski (2000) calculated via cross-correlation, calibrated line full width at half maximum (FWHM), and convolution with rotational broadening profiles, respectively.

The $v\,\sin \,i$ value we report is consistent with the old age of the star found by SPECIES if we consider that stars on the main sequence spin-down with time due to the loss of angular momentum from winds and the increase in stellar radius with time that is required to maintain hydrostatic equilibrium as the core changes due to nuclear burning. The Sun has a rotational velocity of only 1.6 ± 0.3 km s⁻¹ (Pavlenko et al. 2012) and we classify it as a slow rotator. In summary, the values found for the parameters make HD 26965 a good candidate for radial velocity planet search, as it can be considered a quiescent and slowly rotating star. The properties and derived parameters for HD 26965 are summarized in Table 1.

Table 1.  Stellar Parameters of HD 26965

Parameter	Value	Source
R.A. (J2000)	04:15:16.32	SIMBAD
Decl. (J2000)	−07:39:10.34	SIMBAD
mV	4.43	SIMBAD
B–V	0.82	SIMBAD
Distance (pc)	4.98 ± 0.01	van Leeuwen (2007)
Spectral type	K0.5V	Gray et al. (2006)
Mass (M⊙)	0.76 ± 0.03	This work (SPECIES)
Age (Gyr)	9.23 ± 4.84	This work (SPECIES)
Luminosity (L⊙)	0.44	Anderson & Francis (2012)
T${}_{\mathrm{eff}}$ (K)	5151 ± 55	This work (SPECIES)
[Fe/H]	−0.29 ± 0.12	This work (SPECIES)
log g	4.54 ± 0.04	This work (SPECIES)
v sin i (km s−1)	1.23 ± 0.28	This work (SPECIES)
log ${R}_{\mathrm{HK}}^{{\prime} }$	−4.99	Jenkins et al. (2011)

Download table as:  ASCIITypeset image

High-precision Doppler measurements of HD 26965 were carried out using four different spectrographs: The High Resolution Echelle Spectrograph (HIRES) installed on the 10 m Keck Telescope in Hawaii, the Carnegie Planet Finder Spectrograph (PFS) mounted on the 6.5 m Magellan II (Clay) telescope at Las Campanas Observatory, CHIRON mounted on the 1.5 m telescope from the Small to Moderate Aperture Research Telescopes (SMARTS) consortium in Cerro Tololo Interamerican Observatory, and the High Accuracy Radial velocity Planet Searcher (HARPS) installed on the 3.6 m ESO telescope at La Silla Observatory.

3.1. HIRES Observations

3.2. PFS Observations

3.3. CHIRON Observations

3.4. HARPS Observations

We started to examine the radial velocity data by using the traditional periodogram analysis approach to look for any periodicities embedded in the data. We used the generalized version (Zechmeister & Kürster 2009) of the Lomb–Scargle periodogram (hereafter GLS; Lomb 1976; Scargle 1982) where we set up a minimum period of 1 day and a maximum period of 10000 days for the search, with 80000 trial periods evenly spaced in the frequency domain.

Figure 2 (top panel) shows the GLS periodogram of the combined radial velocities of HIRES, PFS, CHIRON, and HARPS. The velocities have been mean subtracted. A maxima at 42 days (marked with a red arrow) clearly exceeds the power threshold of 0.1% significance level. There are also two power maxima close to the 0.1% significance threshold at ∼38 days and at ∼360 days. In the bottom panel of Figure 2, we show the periodogram of the sampling (window function) of the combined radial velocities. The secondary power spectrum peak at 360 days found in the periodogram of the velocities is also present here, and therefore can be attributed to the frequency of sampling. However, the peak at 38 days is not present and therefore further investigation is required to determine the origin of this possible signal, which we discuss below as being due to the rotation period of the star.

Zoom In Zoom Out Reset image size

In addition to the traditional periodogram analysis, we have performed a Bayesian analysis to search for periodic signals embedded in the data. We modeled the radial velocities of HD 26965 following the statistical model defined in Tuomi et al. (2014b) and also applied in Jenkins & Tuomi (2014) where we include the following elements:

• (1)  
  A function describing a k-Keplerian planet model.
• (2)  
  A linear trend term.
• (3)  
  A red-noise model consisting of a p-th order moving average—MA(p)–model with an exponential smoothing.
• (4)  
  Linear correlations with the stellar activity indicators.

We write the statistical model as follows

$$\begin{eqnarray}{y}_{i,j} & = & {\gamma }_{j}+\dot{\gamma }{t}_{i}+{f}_{k}({t}_{i})+{\epsilon }_{i,j}+\displaystyle \sum _{n=1}^{{q}_{j}}{c}_{n,j}{\xi }_{n,i,j}\\ & & +\displaystyle \sum _{l=1}^{p}{\phi }_{j,l}\ \exp \ \left\{\displaystyle \frac{{t}_{i-l}-{t}_{i}}{{\tau }_{j}}\right\}{r}_{i-l,j}\end{eqnarray} \tag{ 1 }$$

where ${y}_{i,j}$ corresponds to the observation at time t_i for the j-th instrument, ${\gamma }_{j}$ is the velocity offset for the j-th data set, $\dot{\gamma }$ is a linear trend term, and ${r}_{i,j}$ denotes the residuals after subtracting the model from the measurement. The function f_k is a superposition of k-Keplerian signals,

$$\begin{eqnarray}&&{f}_{k}({t}_{i})=\displaystyle \sum _{m=1}^{k}{K}_{m}[\cos ({\omega }_{m}+{\nu }_{m}({t}_{i}))+{e}_{m}\cos ({\omega }_{m})]\end{eqnarray} \tag{ 2 }$$

where K_m is the velocity semi-amplitude, ${\omega }_{m}$ is the longitude of pericenter, ${\nu }_{m}$ is the true anomaly and e_m is the eccentricity. ${\nu }_{m}$ is also a function of the orbital period and the mean anomaly ${M}_{0,m}$. Hence, f_k is fully described by K_m, ${\omega }_{m}$, e_m, ${M}_{0,m}$ and P_m, $m\in \{1,\,\ldots ,\,k\}$. The white noise term is denoted by the additive random variable ${\epsilon }_{i,j}$. We assume that there is an excess white noise in each data set with a variance of ${\sigma }_{j}$ such that

$$\begin{eqnarray}&&{\epsilon }_{i,j}\sim { \mathcal N }(0,{\sigma }_{i}^{2}+{\sigma }_{j}^{2})\end{eqnarray} \tag{ 3 }$$

with ${\sigma }_{i}$ the uncertainty associated with the measurement ${y}_{i,j}$ and ${\sigma }_{j}$ is the excess white noise or jitter for the j-th data set, that is treated as a free parameter in the model.

The remaining terms define the rest of the noise model, including the red-noise component: the first term with parameters ${c}_{n,j}$ describes the linear correlations with q stellar activity indicators ${\xi }_{n,i,j}$ for the n-different instruments. The second term is the MA(p) component with smoothing over a timescale ${\tau }_{j}=4$ days and ${\phi }_{j,l}$ with a value between −1 and 1 to quantify the correlation between measurements. The smoothing timescale is set to 4 days for simplicity (Tuomi et al. 2013b). We assume the noise is correlated in this timescale although with higher cadence smaller timescales would likely be more appropriate (Tuomi et al. 2013b; Feng et al. 2016).

5.1. Posterior Samplings and Signal Detection

To estimate the posterior probability of the parameters in the model given the observed data, we use Bayes rule that states

$$\begin{eqnarray}&&P(\theta | y)=\displaystyle \frac{P(y| \theta )P(\theta )}{\,\int P(y| \theta )P(\theta )d\theta }\end{eqnarray} \tag{ 4 }$$

where $P(y| \theta )$ is the probability density of the measurements given the parameters (likelihood function) and $P(\theta )$ corresponds to the prior, i.e., what is known about a given parameter and its constraints before making the measurement. The denominator in Equation (4) is a normalizing constant such that the posterior must integrate to unity over the parameter space.

Table 2.  Summary of Instrumental and Observational Parameters for the Different Instruments

Instrument	Resolution	$\langle \mathrm{SNR}\rangle /\mathrm{Resol}.\,\mathrm{Element}$	$\langle \mathrm{Exposure}\ \mathrm{Time}\rangle$	$\langle \mathrm{RV}\ \mathrm{error}\rangle$	N${}_{\mathrm{obs}}$	Time baseline
			(s)	(m s−1)		(years)
HIRES	45000	270 × 4 exp	11	1.2	230	12
PFS	80000	235 × 4 exp	40	1.0	65	5
CHIRON	95000	120 × 3 exp	300	1.6	259	2
HARPS	115,000	150 × 4 exp	100	0.4	437	10

Download table as:  ASCIITypeset image

In our model, we chose the priors for the orbital and instrumental parameters as listed in Table 3.

Table 3.  Prior Selection for the Parameters

Parameter	Prior Type	Range
Semi-amplitude	Uniform	$K\in [\,0,{K}_{\max }]\,$
Period	Jeffrey's	$P\in [\,1,2{P}_{\mathrm{obs}}]\,$
Eccentricity	${ \mathcal N }(0,{\sigma }_{e})$	$e\in [0,1)$
Long. of Peric.	Uniform	$\omega \in [\,0,2\pi ]\,$
Mean Anomaly	Uniform	${M}_{0}\in [\,0,2\pi ]\,$
Jitter	Uniform	${\sigma }_{J}\in [0,{K}_{\max }]$
Smoothing timescale	Constant	${\tau }_{j}=4$

Download table as:  ASCIITypeset image

In order to investigate the signal initially found with GLS periodogram of the combined radial velocities, we use our Bayesian detection method where we sample the parameter space using the Delayed-Rejection Adaptive-Metropolis (DRAM) algorithm (Haario et al. 2006) based on the Adaptive-Metropolis (AM) algorithm (Haario et al. 2001), applied in Tuomi et al. (2014a) and Jenkins & Tuomi (2014). DRAM and AM are both methods for improving the efficiency of the Metropolis–Hastings algorithm (Metropolis et al. 1953; Hastings 1970). The idea behind using DRAM is that when the posterior of a parameter is multimodal, such as the orbital period in the case of Keplerian fits to radial velocity data sets, and a new state for the chain is rejected (see full details in Tuomi et al. 2014b), a new proposed state is drawn centered on the last one. Up to three rejections are allowed before that part of the posterior is finally discarded as a region of low probability. This has the benefit of sampling more heavily the posterior phase space, at the cost of a longer run-time.

Tempered samplings are also performed when searching for signals. We include a β parameter following Tuomi et al. (2014b), such as $\beta \in (0,1)$, meaning we use $P{(\theta | y)}^{\beta }$ instead of the standard posterior probability density, $P(\theta \,| \,y)$. In this way, we can define the "temperature" of the chain simply as $T=1/\beta$ and so a "hot" chain is defined when $T\gt 1$ and a "cold" chain is where T = 1. When $T\gt 1$ the relative height of the maxima in the posterior probability density are decreased to prevent the chains from getting stuck in regions of high probability, allowing them to visit the entire period parameter space. The typical length of a chain is set to be between ${10}^{6}\mbox{--}{10}^{7}$ for the search run and 10⁶ for the initial burn-in period.

We performed a first run for a zero-planet model to determine the observational baseline, and the instrumental noise and stellar noise parameters for each set of radial velocities. We then searched for a signal in the radial velocity data considering a one-planet model. The search was initially done by setting the temperature for the chain hot enough to let the chain explore the entire parameter space. This is especially helpful when the parameter space is highly multimodal. Our tolerance threshold for the acceptance rate is based on the optimal acceptance rate of the Metropolis–Hastings algorithm, which is ∼0.234 (Roberts et al. 1997). A lower threshold for the chain to be accepted was set to 10%, so hot chains with lower acceptance rates were discarded. From these runs, we found a strong signal was present with a period of 42 days. The signal identified from the maximum of the posterior probability density distribution is shown in Figure 3, left panel.

Zoom In Zoom Out Reset image size

We repeated this process by adding additional signals to the model, but we found no more statistically significant periods in the distribution of the posterior probability densities. Finally, to constrain the detected signal, we performed parameter estimations via the AM algorithm by setting a cold chain (β = 1) with the parameters initially set as a small ball around the parameters found previously by the hot chain run with DRAM. This gave rise to the posterior histograms shown in Figure 4, where the period, amplitude, and minimum mass distributions show nice Gaussian forms centered on their respective values, and the eccentricity distribution is consistent with zero. Table 4 summarizes the final set of values for the parameters from our analysis.

Zoom In Zoom Out Reset image size

Table 4.  Solutions for HD 26965

Parameter	HD 26965  b
P (days)	42.364 ± 0.015
K (m s−1)	1.59 ± 0.15
e	0.017 ± 0.046
ω (rad)	0.31 ± 1.93
M0 (rad)	4.92 ± 1.92
a (au)	0.215 ± 0.008
$m\sin i$ (M${}_{\oplus }$)	6.91 ± 0.79
${\gamma }_{\mathrm{PFS}}$ (m s−1)	0.13 ± 0.80
${\gamma }_{\mathrm{HIRES}}$ (m s−1)	0.50 ± 0.57
${\gamma }_{\mathrm{CHIRON}}$ (m s−1)	0.43 ± 0.79
${\gamma }_{\mathrm{HARPS},\mathrm{old}}$ (m s−1)	0.17 ± 0.50
${\gamma }_{\mathrm{HARPS},\mathrm{new}}$ (m s−1)	0.45 ± 0.78
$\dot{\gamma }$ (m ${{\rm{s}}}^{-1}\,$ year−1)	−0.031 ± 0.037
${\sigma }_{\mathrm{PFS}}$ (m s−1)	1.54 ± 0.20
${\sigma }_{\mathrm{HIRES}}$ (m s−1)	2.38 ± 0.15
${\sigma }_{\mathrm{CHIRON}}$ (m s−1)	1.78 ± 0.15
${\sigma }_{\mathrm{HARPS},\mathrm{old}}$ (m s−1)	1.11 ± 0.05
${\sigma }_{\mathrm{HARPS},\mathrm{new}}$ (m s−1)	0.69 ± 0.07
${\phi }_{\mathrm{PFS}}$	0.82 ± 0.10
${\phi }_{\mathrm{HIRES}}$	0.61 ± 0.07
${\phi }_{\mathrm{CHIRON}}$	0.62 ± 0.06
${\phi }_{\mathrm{HARPS},\mathrm{old}}$	0.81 ± 0.04
${\phi }_{\mathrm{HARPS},\mathrm{new}}$	0.90 ± 0.08
c${}_{{\rm{S}}{\rm{PFS}}}$ (m s−1)	61.1 ± 14.0
c${}_{{\rm{S}}{\rm{HIRES}}}$ (m s−1)	53.1 ± 15.0
c${}_{{\rm{BIS}}{\rm{HARPS}},\mathrm{old}}$	0.086 ± 0.046
c${}_{{\rm{FWHM}}{\rm{HARPS}},\mathrm{old}}$	1.8 ± 4.5
c${}_{{\rm{S}}{\rm{HARPS}},\mathrm{old}}$ (m s−1)	1.6 ± 2.6
c${}_{{{\rm{H}}}_{\alpha }\mathrm{HARPS},\mathrm{old}}$ (m s−1)	−12.1 ± 8.9
c${}_{\mathrm{HeI}{\rm{HARPS}},\mathrm{old}}$ (m s−1)	−76 ± 34
c${}_{{\rm{BIS}}{\rm{HARPS}},\mathrm{new}}$	0.27 ± 0.11
c${}_{{\rm{FWHM}}{\rm{HARPS}},\mathrm{new}}$	0.089 ± 0.018
c${}_{{\rm{S}}{\rm{HARPS}},\mathrm{new}}$ (m s−1)	105 ± 36
c${}_{{{\rm{H}}}_{\alpha }\mathrm{HARPS},\mathrm{new}}$ (m s−1)	76 ± 125
c${}_{\mathrm{HeI}\mathrm{HARPS},\mathrm{new}}$ (m s−1)	23 ± 141

Note. Final set of orbital and instrumental parameters. 1σ errors.

Download table as:  ASCIITypeset image

5.2. Model Selection

It is important to define a robust methodology that allows us to compare the results for two given models in order to address the statistical significance of one model with respect to the other.

The probability of a model ${ \mathcal M }$, containing the best-fit parameters for the observed data y, is given by

$$\begin{eqnarray}&&P({ \mathcal M }| y)=\displaystyle \frac{P(y| { \mathcal M })P({ \mathcal M })}{{\displaystyle \sum }_{i=1}^{k}P(y| { \mathcal M })P({ \mathcal M })}.\end{eqnarray} \tag{ 5 }$$

In particular, we want to know if the model containing one planet is more probable than a zero-planet model and so on for the k-planet model with respect to a k-1-keplerian model. To solve this, we compute the probability of a given model by using the corresponding value of the Bayesian Information Criterion (BIC). A complete and detailed discussion can be found in Tuomi & Jones (2012) and Feng et al. (2016). To compare the ${{ \mathcal M }}_{k}$ model with a previous ${{ \mathcal M }}_{k-1}$ model, we simply compute the logarithm of the Bayes factor, ln ${B}_{k,k-1}$, defined via

$$\begin{eqnarray}&&\mathrm{ln}\,{B}_{k,k-1}=\mathrm{ln}\,P(y| {{ \mathcal M }}_{k})-\mathrm{ln}\,P(y| {{ \mathcal M }}_{k-1}).\end{eqnarray} \tag{ 6 }$$

Furthermore, the model containing the best-fit parameters that support the signal has to fulfill the detection criteria described in Tuomi (2012). It must hold that

$$\begin{eqnarray}&&P(y| {{ \mathcal M }}_{k})={sP}(y| {{ \mathcal M }}_{k-1})\end{eqnarray} \tag{ 7 }$$

where $s\gt {10}^{4}$. Hence, using the Bayes factor defined in Equation (6), we require that the ${{ \mathcal M }}_{k}$ model describing the k-keplerian signal has to be more statistically probable than the ${{ \mathcal M }}_{k-1}$ model associated with the k-1-keplerian signal. Following the conservative threshold from Tuomi et al. (2014b), we define that the evidence ratio should be

$$\begin{eqnarray}&&\mathrm{ln}\,{B}_{k,k-1}\gt 9.2\end{eqnarray} \tag{ 8 }$$

which translates posterior odds of 10000:1 that the k model is selected over the k-1 model, in order to satisfy our detection criteria. Table 5 shows the Bayes factors for $k=0,1,2$-planet models with and without activity correlation terms.

Table 5.  Logarithm of Bayes Factors Comparing a k = 0, k = 1 and k = 2 Keplerian Models With and Without Activity Correlations

Bayes Factor	Activity	No Activity
$\mathrm{ln}{B}_{k,k-1}$	Correlations	Correlations
ln ${B}_{\mathrm{1,0}}$	43.38	35.44
ln ${B}_{\mathrm{2,1}}$	0.61	2.96

Download table as:  ASCIITypeset image

Figure 5 shows the minimum mass detection thresholds for additional planets orbiting around HD 26965. The green-filled area represents the liquid-water habitable zone estimated according to Kopparapu et al. (2013a, 2013b). The thresholds are calculated following the methods in Tuomi et al. (2014b). From this figure we can say that planets with minimum masses in excess of Neptune in the habitable zone can be ruled out meaning if there are HZ planets orbiting HD 26965, they would likely be super-Earths or smaller. The red circle represents the planet candidate, barely, but significantly above the detection threshold.

Zoom In Zoom Out Reset image size

5.3. Signal Injection

As an additional test to investigate if the signal was supported for each instrument, we performed a signal injection on the individual data sets. We use the best-fit parameters (i.e., P, K, ω, M₀, e) of the putative signal by using a Keplerian function described in Equation (2). The hypothesis is the following: if the 42-day signal is injected in a given data set and we run our usual Bayesian analysis, we should, in principle, easily detect it. If the signal is indeed present in the data set, we should recover approximately twice the best-fit radial velocity semi-amplitude, as it has been boosted by the injection. On the other hand, if the recovered velocity semi-amplitude is significantly lower than our best-fit values, that would suggest the actual data is not supported by the instrument, or in other words, the precision of the instrument plus the current number of observations do not allow the signal to be detected. When boosting our signal, we recover the candidate period for HARPS OLD, HIRES and CHIRON data sets. For CHIRON, however, the expected peak at 42 days in the GLS is not unique, although it is above the 0.1% significance level, as can be seen from the periodograms shown in Figure 6. In the case of PFS data, we did not recover the candidate period. Instead, we found a strong power at ∼5 days. This could be caused by the sparse sampling and lower number of observations available from this instrument.

Zoom In Zoom Out Reset image size

To investigate the nature of the detected 42-day signal, we perform a similar analysis as in Santos et al. (2014) on the activity indices available from each instrument. First, we searched for periodicities present in the activity indices themselves, again using the GLS, and we show these results in Figure 7. There are no statistically significant peaks associated with the 42-day signal we detect in the radial velocities. However, the periodogram of the HARPS S-indices shows an emerging peak at  38 days, which is very close to the signal we detected in the radial velocities. Interestingly, this was the period found for the rotation of the star from previous analysis of Ca ii lines Saar & Osten (1997), which also agrees with the period inferred from ROSAT measurements (37.1 days; Pizzolato et al. 2003).

Zoom In Zoom Out Reset image size

Figure 8 shows the correlations between radial velocity and activity indicators for HARPS, PFS, and HIRES. The combined S-indices we show have been mean subtracted and then combined together. We have computed the Pearson Rank test coefficients to determine the correlation between these quantities. Results are listed in Table 6 where we also list the uncertainties associated with each coefficient. To calculate these uncertainties, we ran 10000 bootstraps and created a distribution of r coefficients for every activity index, where the standard deviation of the distribution gave us a measurement of the uncertainty on the coefficients. We note that the correlations are not significant within these uncertainties, given the standard statistical limits for claiming a weak ($r\lt 0.5$), a moderate ($0.5\leqslant r\leqslant 0.7$), and a strong correlation ($r\gt 0.7$); therefore, we can conclude that the stellar activity indicators do not argue against a Doppler origin for the signal, yet the correlations indicate we must consider them in our full statistical model. Indeed, the correlations suggest there is a weak impact of the stellar noise on the velocities, and we confirm this, as the probability of our statistical model is higher when we include these correlations, compared to when we exclude them (see Table 5).

Zoom In Zoom Out Reset image size

Table 6.  Pearson Rank Test Coefficients. Correlation between Activity Indicators and Radial Velocities

Activity Indicator	r
BIS${}_{\mathrm{HARPS}}$	−0.02 ± 0.04
FWHM${}_{\mathrm{HARPS}}$	0.74 ± 0.04
S${}_{\mathrm{HARPS}}$	0.14 ± 0.04
H${}_{\alpha \mathrm{HARPS}}$	0.05 ± 0.04
He ${{\rm{I}}}_{\mathrm{HARPS}}$	−0.04 ± 0.04
S${}_{\mathrm{PFS}}$	0.48 ± 0.12
S${}_{\mathrm{HIRES}}$	0.44 ± 0.07
S${}_{\mathrm{COMBINED}}$	0.23 ± 0.03

Download table as:  ASCIITypeset image

Following a similar approach as in Jenkins & Tuomi (2014), we tested the variability and stability of the signal of our candidate. For this analysis, we only considered the HIRES and HARPS OLD data sets, as both of them have a fairly continuous sampling of Doppler measurements along the ∼16 years of observational baseline. The measurements include a total of 662 unbinned velocities, and we chose JD_s = 2,454,600 as the point to split the data, as this was close to the center of the time baseline of the observations and also produced a well balance between HARPS and HIRES data (i.e., not biased to an instrument in particular). The data prior to JD_s contained 408 data points and the data set after the split point contained 254 measurements.

We performed the Bayesian analysis on these two subsets of velocities independently, running cold chains to constrain the orbital parameters of a 1-Keplerian model. We found the signal is detected with values in agreement within uncertainties for the two baselines tested, as well as for the full data set described above. This shows us that the signal is not varying in time and thus the period and amplitude of our planetary candidate is stable over the tested observational baseline, another strong argument against an activity origin since activity processes should be quasi-static, varying over a few rotation periods of the star.

To complement the analysis, we gathered photometric data available from the ASAS Catalog (Pojmanski 1997) to investigate if any periodic signal could be seen in light curves, particularly the rotational period of the star. As mentioned above, the literature values were reported by Noyes et al. (1984) and Saar & Osten (1997), where they found ${P}_{\mathrm{rot}}\,\sim 37.10$ days for this old star. We show the ASAS photometric measurements in Figure 9 and the values are tabulated in Table 12. From the five different apertures available, we selected aperture 1 as its MAD¹⁰ value of 0.219 mag was the smallest. The mean uncertainty in the V-band photometry is ${\sigma }_{\mathrm{ASAS}}=0.036$ mag. Gray circles correspond to the entire photometry set of 568 useful points acquired from 2000 to 2009. However, we excluded the data with poor quality (those not marked "A" or "B" in the catalog) and also those measurements that deviated more than 3σ with respect to the mean value of the time series. The highest quality data (316 points) are shown as black circles in the top panel of Figure 9. The bottom panel in Figure 9 again shows the GLS for the ASAS photometry. We sample the period space starting at a minimum period of 1 day and up to 10,000 days, performing 80,000 period samples. Considering just the data before JD = 2,452,300 tends to favor peaks with higher power toward high frequencies (periods ∼1 day) but without any significant period (or an integer multiple) near the period associated with the 42-day signal of the reported planetary candidate. We also ran the periodogram analysis on the full photometric data set with no significant periods found. Following the relations in Hatzes (2002), we found that a filling factor of f = 0.15 would be required to induce the RV amplitude of 1.6 m s⁻¹ of the signal found in the combined data. If we consider the spots on the surface of the star to be opaque, for the sake of simplicity, the ratio between the stellar flux and the flux considering spots covering 0.15% of the surface of the star would be 0.9985. This means, the loss of light due to spots on the stellar surface can be translated into a Δm = 1.64 mmag. Given the precision of the ASAS photometry for this star, we conclude it is insufficient to be informative.

Zoom In Zoom Out Reset image size

Table 10.  HARPS OLD Radial Velocities of HD 26965

BJD	RV	σ RV	BIS	FWHM	S	H${}_{\alpha }$	He i
	(m s−1)	(m s−1)	(m s−1)	(m s−1)	(dex)	(dex)	(dex)
2452939.80613	−0.434	0.519	1.254	5896.904	0.1753	0.1161	0.5067
2452939.80685	0.231	0.512	1.279	5897.505	0.1742	0.1163	0.507
2452939.80756	−0.367	0.572	1.387	5901.748	0.1747	0.1163	0.505
2452939.80827	0.176	0.524	1.720	5898.172	0.1759	0.1153	0.5061
2452939.80899	2.305	0.954	3.537	5899.064	0.1726	0.1152	0.5039
2452939.80969	1.191	0.612	1.192	5899.194	0.1722	0.1163	0.5075
2452940.76906	−3.630	0.378	0.184	5898.594	0.175	0.1157	0.5076
2452945.76432	−3.375	0.353	2.281	5896.436	0.171	0.1153	0.5067

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Table 11.  HARPS NEW Radial Velocities of HD 26965

BJD	RV	σ RV	BIS	FWHM	S	H${}_{\alpha }$	He i
	(m s−1)	(m s−1)	(m s−1)	(m s−1)	(dex)	(dex)	(dex)
2457274.86039	−0.1517	0.1828	18.5315	5951.1924	0.1777	0.1170	0.5055
2457274.86304	−0.0462	0.1852	18.4779	5951.7017	0.1768	0.1171	0.5059
2457274.86566	0.4185	0.2041	18.6366	5951.8066	0.1777	0.1174	0.5059
2457274.86845	0.3169	0.2026	18.1923	5951.7324	0.1778	0.1173	0.5051
2457277.84809	0.4107	0.2189	18.9974	5950.3081	0.1763	0.1179	0.5056
2457277.85019	0.1997	0.2209	18.6706	5950.1743	0.1777	0.1176	0.5049

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

Given that we find some moderate correlations between the spectral activity indicators and the radial velocities, we supplemented our activity analysis by studying the original Ca ii H&K data from the Mount Wilson Observatory HK Project (Wilson 1978). The Project data are publicly available from the NSO archive¹¹ and include more than 2000 stars observed from 1966 to 1995. There are 1155 HK observations for HD 26965 from JD = 2,439,787.8 to 2,449,771.7. The processed data do not include associated uncertainties to the calibrated S-values. According to Duncan et al. (1991), the uncertainties in the Mount Wilson S-values can be calculated using the weights, W, included in the data that are derived from the photon counts of the measurements. The uncertainty in the S-index measurements is simply defined as ${\sigma }_{S}=S{(\sqrt{W})}^{-1}$. We applied this formula to the reported weights to provide proper uncertainties for the measurements of this star. All HK values for HD 26965 can be found in Table 13.

Table 13.  Mount Wilson HK Project measurements of HD 26965

S${}_{\mathrm{MW}}$	JD—2444000	W
0.223	−523.2	15956.8
0.207	−522.2	16877.2
0.234	−457.2	7264.7
0.206	−219.2	17620.2
0.218	−175.2	17201.9
0.213	−172.2	17475.1
0.219	−107.2	34054.6
0.205	172.8	8716.2
0.227	190.8	17363.2
0.221	200.8	10277.2
0.214	224.8	8639.2
0.205	260.8	8593.6
0.217	278.8	8700.7
0.194	482.9939	3409.7

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We proceeded to run the same periodogram analysis as for the radial velocities and the photometric time series from ASAS. A clear long-term variability of ∼4100 days can be identified (see top and second top panels in Figure 10), providing evidence for a long-period magnetic activity cycle, similar to the long-period solar cycle. After removing this signal by modeling it with a sinusoidal function, a second period of ∼715 days is found in the periodogram (second bottom panel in the figure), likely representing another, shorter period magnetic cycle. Finally, the bottom panel in the figure shows the periodogram of the Mount Wilson S-values after removing the 4100 and 715 day period signals from the data. In this residual periodogram, a signal of ∼42.3 days remains in the data. This is most likely the value reported in Baliunas et al. (1996). The peak is clearly not unique, casting some doubt on its reality, but given that it matches the detected signal in our radial velocity data sets, and rotation periods are known to be quasi-period due to differing spot patterns, changing stellar activity levels, and differential rotation, we must entertain the real possibility that this is actually the rotational period of the star, and not the 38-day period that we found in the measured spectroscopic activity indices. If this is the reality, then HD 26965 represents a case where most of the current suite of tests that we employ to detect planets using radial velocity analyses, fail to remove the noise introduced by the rotation of the star, meaning that we now require better methods to be employed on stars where there are clear correlations between the radial velocities and various activity indicators if we want to detect planets that induce amplitudes at the ∼1 m s⁻¹ level.

Zoom In Zoom Out Reset image size

Through the application of generalized Lomb–Scargle periodograms and tempered MCMC samplings, we conclude that there is a strong periodic signal in the radial velocities of the quiescent and slow rotating K dwarf HD 26965. If interpreted as the Doppler signal induced on the star by an orbiting planet, our best solution explains these variations by the presence of a low-mass, super-Earth planetary candidate that has a minimum mass of 6.92 ± 0.79 ${M}_{\oplus }$ orbiting the host star with a period of 42.364 ± 0.015 days and at a distance of 0.215 ± 0.008 au.

Analysis of V-band photometry from ASAS does not show any significant periodic signal. However, since the amplitude of the signal is small, the precision of the data is not sufficient to detect the signal within the noise of the photometry.

Analysis of the stellar activity indicators does not show statistical evidence supporting a chromospheric origin for the periodic variations in the radial velocities, although we have found correlations between the radial velocities and the activity indices from the different spectrographs. However, when we analyze the independently acquired chromospheric calcium S-indices from the Mt. Wilson HK project, and after removing two long-period activity cycles, we find evidence for the rotation of the star closely matching the period of the radial velocity detected signal.

Regarding this last point, we note that although it is important to properly include activity correlations into any global model of radial velocity data, which when done for this data set, we find a higher statistical probability for the given Keplerian model supporting the planetary signal, if there are statistically significant correlations between activity indicators and the velocity measurements, then additional external activity indicators should be acquired, where possible. Also, moving away from linear correlation models between current activity indices and the radial velocities may be necessary, particularly if the data suggests more complex models, such as quadratics, might be favored. In any case, it is clear that the inclusion of multiple sources of external data that also rule out possible magnetic cycles and rotation periods as the source of any radial velocity signal, can help to maintain the lowest false-positive rate for any given Doppler survey.

In summary, despite all of the evidence favoring a Doppler signal present in this radial velocity data set, the methods described in this paper do not seem to be able to disentangle weak planetary signals from residual photospheric noise, at least when the orbital periods are close to the rotation period of the star and there are correlations present between the velocities and the measured activity indicators.

M.R.D. acknowledges the support of CONICYT-PFCHA/Doctorado Nacional-21140646, Chile. J.S.J. acknowledges support by Fondecyt grant 1161218 and partial support by CATA-Basal (PB06, CONICYT). M.G.S. is supported by CONICYT-PFCHA/Doctorado Nacional-21141037, Chile. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of MaunaKea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.