The HARPS search for southern extra-solar planets XLV. Two Neptune mass planets orbiting HD 13808: a study of stellar activity modelling's impact on planet detection
E. Ahrer, D. Queloz, V. M. Rajpaul, D. Ségransan, F. Bouchy, R. Hall, W. Handley, C. Lovis, M. Mayor, A. Mortier, F. Pepe, S. Thompson, S. Udry, N. Unger
MMNRAS , 1–17 (2020) Preprint 9 February 2021 Compiled using MNRAS L A TEX style file v3.0
The HARPS search for southern extra-solar planets XLV.Two Neptune mass planets orbiting HD 13808: a study of stellaractivity modelling’s impact on planet detection ★ E. Ahrer , † , D. Queloz , , V. M. Rajpaul , D. Ségransan , F. Bouchy , R. Hall ,W. Handley , , C. Lovis , M. Mayor , A. Mortier , , F. Pepe , S. Thompson ,S. Udry , N. Unger Astrophysics Group, Cavendish Laboratory, JJ Thomson Avenue, CB3 0HE Cambridge, UK Department of Physics, University of Warwick, Gibbet Hill Road, CV4 7AL Coventry, UK Departement d’astronomie, Université de Genève, Chemin des Maillettes 51, CH-1290 Versoix, Switzerland Kavli Institute for Cosmology, Cambridge, Madingley Road, CB3 0HA Cambridge, UK
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
We present a comprehensive analysis of 10 years of HARPS radial velocities of the K2V dwarfstar HD 13808, which has previously been reported to host two unconfirmed planet candidates.We use the state-of-the-art nested sampling algorithm PolyChord to compare a wide varietyof stellar activity models, including simple models exploiting linear correlations betweenRVs and stellar activity indicators, harmonic models for the activity signals, and a moresophisticated Gaussian process regression model. We show that the use of overly-simplisticstellar activity models that are not well-motivated physically can lead to spurious ‘detections’of planetary signals that are almost certainly not real. We also reveal some difficulties inherentin parameter and model inference in cases where multiple planetary signals may be present.Our study thus underlines the importance both of exploring a variety of competing models andof understanding the limitations and precision settings of one’s sampling algorithm. We alsoshow that at least in the case of HD 13808, we always arrive at consistent conclusions abouttwo particular signals present in the RV, regardless of the stellar activity model we adopt;these two signals correspond to the previously-reported though unconfirmed planet candidatesignals. Given the robustness and precision with which we can characterize these two signals,we deem them secure planet detections. In particular, we find two planets orbiting HD 13808at distances of 0 . , .
26 AU with periods of 14 . , . , 𝑀 ⊕ . Key words: methods: data analysis – methods: statistical – techniques: radial velocities –stars: activity – stars: individual: HD 13808
The Radial Velocity (RV) method has been an important and pro-ductive tool for discovering exoplanets ever since it led to the dis-covery of the first exoplanet orbiting a Sun-like star (Mayor &Queloz 1995). However, the search for small Earth- and Neptune-like planets orbiting Sun-like stars is very challenging. They giverise to relatively small radial velocity signatures of order 1 m s − or less, whilst stellar magnetic activity can induce RV signals that ★ Based on observations made with HARPS spectrograph on the 3.6-mESO telescope at La Silla Observatory, Chile † E-mail: [email protected] can mimic planetary ones, with amplitudes of order many metersper second (e.g Queloz et al. 2001).Thus, with the increasing precision of RV measurements facil-itated by upcoming extreme-precision Doppler spectrographs suchas ESPRESSO (Pepe et al. 2014), EXPRES (Jurgenson et al. 2016),HARPS3 (Thompson et al. 2016), HIRES (Pasquini et al. 2008)and NEID (Schwab et al. 2016), the characterization of the influ-ence from host stars’ activity on the RV measurements becomesever more important. Hence there is a strong interest in developingtools to disentangle stellar and planetary signals in RV data.A variety of methods exists that attempt to ‘correct’ or modelstellar activity signals in RVs, and thus reduce the possibility offalse-positive planet detections. Models describing this include pre-whitening and red-noise models (e.g Hatzes et al. 2010; Feroz & © a r X i v : . [ a s t r o - ph . E P ] F e b E. Ahrer et al.
Hobson 2013), as well as studying stellar activity indicators suchas the bisector inverse slope (BIS) and Full Width Half Maximum(FWHM) of the cross-correlation function (CCF) between a targetstar and a template spectrum, or measurements of chromosphericactivity in the target star, such as the log 𝑅 (cid:48) HK index (e.g Boisse et al.2009; Queloz et al. 2009; Dumusque et al. 2011b). The presenceof periodic signals in indicators usually suggest an activity-inducedsignal rather than a planetary one, since a genuine planet wouldinduce a periodic Doppler shift in all spectral lines, but would notproduce the same periodic variations in the stellar activity indi-cators. It is straightforward to use linear correlations (e.g. Quelozet al. 2001) and to include harmonic models which use the rotationalperiod of the star and its harmonics (e.g Dumusque et al. 2012) todescribe stellar activity induced RV variations.Other, more computationally expensive methods include mod-elling the stellar surface features directly (e.g Boisse et al. 2012)and predicting activity-related RV variability with stellar activityindicators using the 𝐹𝐹 (cid:48) method and Gaussian processes (GPs) asdescribed in e.g. Aigrain et al. (2012) and Rajpaul et al. (2015);such approaches have already successfully been used to identifyfalse-positive planetary detections as e.g. in Haywood et al. (2014b)or Rajpaul et al. (2016).To quantify the quality of different models for describing ob-served RV data, Bayesian model comparison has proven to be arobust approach (e.g Feroz & Hobson 2013; Faria et al. 2016; Hallet al. 2018), also in combination with GP modelling (Faria et al.2020). However, a significant challenge inherent in this method isvery high computing costs when dealing with high-dimensionalproblems, since computing Bayesian model evidences entails inte-grating over all model parameter posteriors; a related challenge isdegeneracy inherent in combinations of Keplerian and stellar ac-tivity models (see Nelson et al. 2020). Nevertheless, it is widelyaccepted that Bayesian model comparison is a theoretically-soundapproach to answering questions about competing physical mod-els, and that such computational burdens are therefore a priceworth paying (Goodman 1999).This paper presents a comprehensive analysis of a series of246 spectroscopic measurements on HD 13808 carried out withHARPS from 2003 to 2014. This star belongs to a sample of brightstars selected for their low level of RV “jitter” to maximize thesensitivity of the survey to low-mass exoplanets. Series of stellarspectrum of HD 13808 with a signal to noise (SN) ratio large enoughto reach the photon noise RV measurement error of order 50 cm s − are available. In addition, the exposure time was at least 15 minutesto minimize the effect of acoustic mode stellar oscillations (e.g.Dumusque et al. 2011a; Chaplin et al. 2019). Early in the survey ofHD 13808, small amplitude RV variations were detected, suggestinga combination of time variable signal due to a multi-planet systemsand moderate stellar activity. We present a comprehensive analysisof HARPS RVs on that star using Bayesian inference to combinevarious stellar activity models with a set of independent Keplerianorbit solutions.With this paper we aim to examine the importance of compar-ing systematically a number of competing physical models, since theusefulness of Bayesian model comparison is limited by the qualityof the models considered (favouring one inappropriate model overanother does not imply that either model is in any sense ‘correct’). The same can not be said for various other commonly-used approachesto model selection, such as residual minimization, Bayesian informationcriterion testing, etc. (Gelman et al. 2013)
Table 1.
Description of the HARPS RV data of HD 13808.HD13808Number of observations 246Time span (days) 4051mean(RV) (m s − ) 41 095rms(RV) (m s − ) 3.94mean( 𝜎 RV ) (m s − ) 0.76median(log 𝑅 (cid:48) HK ) -4.90 We show that if an inadequate stellar activity model is considered,one may be led astray and end up with spurious conclusions aboutthe number of planets present in one’s system. Models we considerin our analysis include: linear correlations between RV data and stel-lar activity indicators; sinusoidal models; and joint GP regressionof RV measurements and activity. Moreover, we present a thoroughtest of the nested sampling algorithm PolyChord, a state-of-the-art sampler, and discuss its limitations and the importance of theprecision settings in exoplanet applications. The HD 13808 systemis used as a test case for this study as (i) the host star shows stellaractivity variability, and (ii) it has been suggested in the literaturethat the system hosts at least two planet candidates, though the ex-istence of these planets was never securely confirmed (see Mayoret al. 2011 and Gillon et al. 2017; at the time of writing this paper,the NASA Exoplanet Archive lists HD 13808 as having no con-firmed planets). With the analysis in this paper the status of thesecandidates elevates to confirmed planets.This paper is structured as follows. We describe the HARPSobservations in Section 2, followed by the modelling of Keplerianorbits in Section 3. We introduce the stellar activity models ap-plied in this work in Section 4 and outline model comparison usingBayesian inference and PolyChord in Section 5. Afterwards, wepresent the results of our analysis of HD 13808 RVs in Section 6and conclude with our discussion in Section 7. Spectra of HD 13808 were measured with HARPS, a fiber-fed spec-trograph installed in a vacuum vessel, mounted on the 3.6m ESOtelescope of La Silla in Chile (Mayor et al. 2003). Observationswere made using the most precise observing mode which utilizes asimultaneous calibration by the ThAr calibration lamp. In this con-figuration HARPS achieves a long term RV precision better than1 m s − allowing us to detect small stellar RV variations of this order(Lovis et al. 2008).The data obtained by HARPS are automatically processed onsite by a data reduction software – the HARPS DRS – that extractsthe spectra, calibrates it and eventually computes a cross-correlationfunction with a stellar template. The stellar radial-velocity is mea-sured from the CCF as well as other parameters like the FWHM orthe BIS Queloz et al. (2001); Santos et al. (2002). In addition foreach HARPS spectra, the value of the Calcium S activity index isestimated from the chromospheric re-emission in the Ca II H and Klines and converted in standard log 𝑅 (cid:48) HK index.Since the release of the first version of the HARPS DRS in2003, a series of successive versions with improved algorithmshave been developed, leading to a steady gain in the RV precision. Available online at exoplanetarchive.ipac.caltech.edu .MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 Table 2.
HD 13808: Stellar properties, with corresponding references (1)Gray et al. (2006); (2) Kharchenko (2001); (3) Skrutskie et al. (2006); (4)Santos et al. (2013); (5) Delgado Mena et al. (2019); (6) GAIA DR2 Bailer-Jones et al. (2018) HD13808 ReferenceSpectral Type K2V (1)V (mag) 8 . ± .
01 (2)K (mag) 6 . ± .
02 (3)T eff (K) 5035 ±
50 (4)Fe/H (dex) − . ± .
02 (4) 𝑀 ( 𝑀 (cid:12) ) . ± .
022 (5)Age (Gyr) 7 . ± . 𝑅 𝑆 ( 𝑅 (cid:12) ) . + . − . (6)d (pc) 28 . ± . Details about the historical changes in the HARPS DRS algorithmsmay be found in following series of papers: Baranne et al. (1996);Pepe et al. (2002); Lovis & Pepe (2007); Mayor et al. (2009a,b).Measurements used in this paper were reprocessed with version3.5 of the HARPS DRS A description of the HARPS data used in this study is givenin Table 1 while HD 13808’s stellar parameters are summarisedin Table 2. HD 13808 has been observed by HARPS 246 timesover a span of more than 10 years, achieving consistent quality ofmean( 𝜎 RV ) = 0.76 m s − . To model the observed RV of a star at time 𝑡 𝑖 with 𝑁 𝑝 planets,we sum the planets’ individual Keplerian terms, neglecting planet-planet interactions; following Feroz et al. (2011) and the formalismgiven in Balan & Lahav (2009) we arrive at:RV ( 𝑡 𝑖 ) = 𝑉 𝑖 + 𝑁 𝑝 ∑︁ 𝑝 = 𝐾 𝑝 (cid:2) cos ( 𝑓 𝑖, 𝑝 + 𝜔 𝑝 ) + 𝑒 𝑝 cos ( 𝜔 𝑝 ) (cid:3) , (1)where 𝑉 𝑖 is the systemic velocity, 𝐾 𝑝 the RV semi-amplitude, 𝑓 𝑖, 𝑝 ,the true anomaly, 𝑒 𝑝 the orbital eccentricity, and 𝜔 𝑝 the argumentof periastron of the 𝑝 th planet, respectively. In addition, the intro-duction of the mean longitude 𝜆 𝑝 of the 𝑝 th planet is necessary aspart of the computation of 𝑓 𝑖, 𝑝 which requires 𝑒 𝑝 and the period ofthe planet 𝑃 𝑝 as well.A planet’s semi-amplitude 𝐾 𝑝 is related to its mass 𝑀 𝑝 , periodand orbital eccentricity 𝑒 𝑝 , as well as to the mass of the parent star 𝑀 𝑆 and the inclination of the system 𝑖 relative to the observer: 𝐾 𝑝 = . − √︃ − 𝑒 𝑝 (cid:18) 𝑃 𝑝 (cid:19) − / 𝑀 𝑝 sin ( 𝑖 ) 𝑀 Jup (cid:18) 𝑀 𝑆 𝑀 (cid:12) (cid:19) − / . (2)In summary, we have five free parameters per planet ( 𝐾 , 𝜔 , 𝑒 , 𝑃 , and 𝜆 ) plus a white-noise ‘jitter’ term 𝜎 + RV which is added to theobservational error in quadrature (see Section 5.2) and any termsfor describing the systemic velocity, 𝑉 – typically a constant offsetplus a linear or possibly quadratic polynomial term. For the analysisin this paper, a quadratic polynomial term was used for all models.Additionally, our models for describing stellar activity contributions All data used in the analysis will be available via VizieR at CDS.
Table 3.
Priors for the Keplerian, RV trend and additive white noise (‘jitter’)parameters used in all of our models and analyses.* We also require for the periods to be sorted i.e. 𝑃 𝑏 < 𝑃 𝑐 etc.** We also require that the corresponding eccentricity e < 1.*** This prior applies to all polynomial terms.Parameter Prior Lower Bound Upper Bound 𝑃 (d)* Log Uniform 5 100 𝐾 (m s − ) Log Uniform 0 . √ 𝑒 sin ( 𝜔 ) ** Uniform − √ 𝑒 cos ( 𝜔 ) ** Uniform − 𝜆 (rad) Uniform 0 2 π 𝑉 (m s − )*** Uniform − RV max RV max 𝜎 + RV (m s − ) Uniform 0 20 to RVs contain anything from one to sixteen free parameters; thesemodels are described in detail in Section 4.The uninformative priors for our five free parameters per planetare listed in Table 3, as well as the priors for the polynomial termsand for the additive white-noise term. Note that the lower and upperboundary for planet periods was set to 5 and 100 d for computationalefficiency. This choice was supported, in the first instance, by thefact that the Lomb-Scargle power spectrum of RVs revealed nosignificant periodicities above about 60 d (see Fig. 3). Moreover,models allowing planets with periods longer than 100 d or shorterthan 5 d always had lower Bayesian evidences than models withoutsuch planets (when using a GP stellar activity model – cf. Section 6),and the inferred RV semi-amplitudes of these short- or long-periodplanets was always consistent with zero. The priors used for stellar activity model parameters weregenerally also chosen to be uninformative, and are discussed inSection 4 and summarized in Table 4.
In this section we introduce the multiple ways of modelling stellaractivity applied in this work. Simple models like linear dependen-cies are considered, as well as harmonic modelling of the rotationperiod of the star, combining BIS and RV measurements in a si-multaneous fit. We also discuss more sophisticated approaches tomodelling activity-induced RV variations, namely the 𝐹𝐹 (cid:48) approachand simultaneous GP regression over multiple activity indicators. A linear relation between RVs and BIS was considered as seenin equation (3), where RV
Kepler ( t ) are the RV signatures of theplanets at time 𝑡 as described previously in equation (1); BIS ( t ) represents the BIS measurements taken at time 𝑡 and 𝛼 describesthe free parameter for the linear relation. A quadratic polynomialterm including a constant offset is represented by 𝑉 .RV total ( 𝑡 ) = RV Kepler ( 𝑡 ) + 𝛼 BIS ( 𝑡 ) + 𝑉. (3)This method has already been successful in identifying false However, the computational cost of evaluating models with expandedplanet-period priors was typically an order of magnitude greater than whenusing the more restricted prior, and posterior multi-modality was far morepronounced.MNRAS000
Kepler ( t ) are the RV signatures of theplanets at time 𝑡 as described previously in equation (1); BIS ( t ) represents the BIS measurements taken at time 𝑡 and 𝛼 describesthe free parameter for the linear relation. A quadratic polynomialterm including a constant offset is represented by 𝑉 .RV total ( 𝑡 ) = RV Kepler ( 𝑡 ) + 𝛼 BIS ( 𝑡 ) + 𝑉. (3)This method has already been successful in identifying false However, the computational cost of evaluating models with expandedplanet-period priors was typically an order of magnitude greater than whenusing the more restricted prior, and posterior multi-modality was far morepronounced.MNRAS000 , 1–17 (2020)
E. Ahrer et al.
Table 4.
Priors for the parameters of the various stellar activity modelsdescribed in Section 4.* 𝑅 𝑆 is the radius of the star and 𝜎 𝑆 the estimated error. + 𝜎 is the observed standard deviation of the FWHM.Parameter Prior Lower Bound Upper Bound 𝛼 Uniform − RV max + RV max 𝑃 rot Uniform 30 d 42 d 𝑃 magn Uniform 500 d 4500 d 𝜙 Uniform 0 2 π 𝐶 BIS
Uniform − BIS max + BIS max 𝐶 magn Uniform − log 𝑅 (cid:48) HK , max + log 𝑅 (cid:48) HK , max 𝛽 Uniform − RV max + RV max Ψ Uniform FWHM max
FWHM max + 𝜎 + 𝛿𝑉 𝑐 𝜅 Uniform 0 10 000 𝑅 ∗ Uniform 𝑅 𝑆 − 𝜎 𝑆 * 𝑅 𝑆 + 𝜎 𝑆 * 𝑉 𝑐 Uniform − rms ( RV ) + rms ( RV ) 𝑉 𝑟 Uniform − rms ( RV ) + rms ( RV ) 𝐿 𝑐 Uniform − rms ( log R (cid:48) HK ) + rms ( log R (cid:48) HK ) 𝐵 𝑐 Uniform − rms ( BIS ) + rms ( BIS ) 𝐵 𝑟 Uniform − rms ( BIS ) + rms ( BIS ) 𝑃 GP Uniform 10 d 100 d 𝜆 p Jeffreys 0 .
01 10 𝜆 e Jeffreys 10 d 400 d positives (e.g Queloz et al. 2001). However, note that it has beenshown that there can be a temporal offset between activity andassociated RV variations on time scales of few days causing a ‘blur’of any linear correlation between those two parameters, making alinear model a choice with a major caveat (e.g Santos et al. 2014;Collier Cameron et al. 2019). The prior distribution for 𝛼 is shownwith the other parameters related to stellar activity modelling inTable 4.In a similar way, a linear dependency on the log 𝑅 (cid:48) HK activityindicator was considered as we hoped to capture any possible linearlong term relationship withRV total ( 𝑡 ) = RV Kepler ( 𝑡 )+ 𝛽 log 𝑅 (cid:48) HK ( 𝑡 ) + 𝑉 (4)where 𝛽 represents the linear factor and log 𝑅 (cid:48) HK ( 𝑡 ) thelog 𝑅 (cid:48) HK measurements at time 𝑡 . The prior distribution for 𝛽 canbe found in Table 4. Our second activity model entailed fitting a sinusoid to both theRVs and BIS time series, enforcing an identical period but allowingdifferent phases and amplitudes between the two. We assume thatthis captures solar spots and other stellar features which are sensitiveto the rotation of the star. In addition, it accounts for a likely timeshift which is a problem when considering linear correlations asmentioned in Section 4.1.Equation (5) and (6) show the two models with the harmon-ics of the rotation period simultaneously fit to the RV and BISmeasurements, respectively. The parameter 𝑃 rot is the putative rota-tion period of the star fit to both models, while the other parameters 𝐾 𝑖 , 𝜙 RV and 𝐾 𝑘 , 𝜙 BIS and 𝐶 BIS describe the sinusoidal signal of the 𝑖 th and 𝑘 th harmonic for the RV ( 𝑡 ) and BIS ( 𝑡 ) model respectively, with 𝑁 ℎ being the number of harmonics. As before, the parameter 𝑉 represents a quadratic polynomial contribution.RV total ( 𝑡 ) = RV Kepler ( 𝑡 ) + 𝑁 ℎ ∑︁ 𝑖 = 𝐾 𝑖 sin (cid:18) π 𝑖𝑡𝑃 rot + 𝜙 RV (cid:19) + 𝑉 (5)BIS ( 𝑡 ) = 𝑁 ℎ ∑︁ 𝑘 = 𝐾 𝑘 sin (cid:18) π 𝑘𝑡𝑃 rot + 𝜙 BIS (cid:19) + 𝐶 BIS (6)The prior distribution for 𝐾 is as in Table 3; the priors for 𝜙, 𝐶 BIS and 𝑃 rot are displayed in Table 4. Note that the prior on therotation period 𝑃 rot was chosen to cover a narrow range based on thefollowing: (i) the rotation period for HD 13808 was estimated beforeto be ∼
40 d based on the average of log 𝑅 (cid:48) HK measurements by Loviset al. (2011); (ii) signals likely corresponding to the harmonics ofthe rotation period are detected in the periodograms of the stellaractivity indicators, see Section 6.1 (Fig. 3). Further motivation wasprovided by the fact that (iii) when considering zero or 1-planetmodels, the rotation period of the star became locked on to theperiod of one of the two planets, and that (iv) during preliminaryruns, the MAP value for the rotation period never went above 40 d.By limiting the period to below 42 d we reduced the computationalburden of our modelling. In a similar fashion to the harmonic activity model, the ‘magn.cycle’ model fits a sinusoid to two data sets, this time to the RVs andthe log 𝑅 (cid:48) HK time series, enforcing an identical period but allowingdifferent phases and amplitudes. We assume that this captures stellaractivity corresponding to a long-term magnetic activity cycle.Equation (7) and (8) show the two models with the cycle period 𝑃 magn simultaneously fit to the RV and log 𝑅 (cid:48) HK measurements,respectively. The parameters 𝐾 magn , RV , 𝜙 magn , RV and 𝐾 magn , 𝜙 magn and 𝐶 magn describe the RV ( 𝑡 ) and log R (cid:48) HK ( 𝑡 ) model respectively; 𝑉 represents a constant offset and a quadratic polynomial term.RV total ( 𝑡 ) = RV Kepler ( 𝑡 ) + 𝐾 magn , RV sin (cid:18) π 𝑡𝑃 magn + 𝜙 magn , RV (cid:19) + 𝑉 (7)log R (cid:48) HK ( 𝑡 ) = 𝐾 magn sin (cid:18) π 𝑡𝑃 magn + 𝜙 magn (cid:19) + 𝐶 magn (8)The prior distribution for 𝐾 is as in Table 3; the priors for 𝜙, 𝐶 magn and 𝑃 magn are displayed in Table 4. The upper and lowerprior limits for 𝑃 magn were determined by a broad peak around1000—4000 d in the periodogram of log 𝑅 (cid:48) HK as well as by theobvious periodic signal in the log 𝑅 (cid:48) HK measurements (see Fig. 2). 𝐹𝐹 (cid:48) method The FWHM measurements were used for applying the 𝐹𝐹 (cid:48) methodwhich was introduced by Aigrain et al. (2012) as a method for re-lating the photometric brightness and RV variations of a star. Ituses the flux of the star Ψ ( 𝑡 ) and its derivative (cid:164) Ψ ( 𝑡 ) as an indica-tor for spot coverage and predicts RV variations. These variationsinclude the RV perturbation Δ RV rot ( t ) due to the presence of spots MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 on the rotating photosphere and their effect of the suppression ofconvective blueshift Δ RV conv ( t ) .Aigrain et al. (2012) describes these with the following twoequations where 𝑓 represents the drop in flux produced by a spotat the centre of the stellar disc, 𝑅 ∗ is the stellar radius, 𝛿𝑉 𝑐 isthe difference between the convective blueshift in the unspottedphotosphere and that within the magnetized area and 𝜅 is the ratioof this area to the spot surface. Δ RV rot ( t ) = − Ψ ( t ) Ψ (cid:20) − Ψ ( t ) Ψ (cid:21) R ∗ f (9) Δ RV conv ( t ) = (cid:20) − Ψ ( t ) Ψ (cid:21) 𝛿 V c 𝜅 f (10)The total RV variation Δ RV activity created by stellar activity isthen the sum of both terms: Δ RV activity = Δ RV rot ( t ) + Δ RV conv ( t ) (11)The parameter 𝑓 is approximated as 𝑓 ≈ Ψ − Φ min Ψ (12)with Φ min being the minimum observed flux.It has been argued that CCF FWHM andlog 𝑅 (cid:48) HK measurements should both behave, to first order, asthe convective blueshift suppression term in the 𝐹𝐹 (cid:48) method,which in turn is a close proxy for the integrated active regioncoverage on the visible stellar hemisphere (Rajpaul et al. 2015);indeed, it has been shown in practice that the CCF FWHM canbe a good tracer of photometric flux (e.g. Suárez Mascareñoet al. 2020). Thus, as we did not have photometric measurementsof HD 13808 available, we chose to interpret the FWHM as aproxy for the flux Ψ ( 𝑡 ) . The derivative of the ‘flux’ was computednumerically by modelling the FWHM with a GP using the celeritealgorithm by Foreman-Mackey et al. (2017). In total, there werethree free parameters in this stellar activity model, 𝛿𝑉 𝑐 𝜅 , 𝑅 ∗ and Ψ - their respective priors appear in Table 4. We used the GP framework developed by Rajpaul et al. (2015),hereafter R15, to model RVs simultaneously with log 𝑅 (cid:48) HK and BISobservations. In short, this framework assumes that all observed stel-lar activity signals are generated by some underlying latent function 𝐺 ( 𝑡 ) and its derivatives; this function, which is not observed directly,is modelled with a Gaussian process (Rasmussen & Williams 2006;Roberts et al. 2013).Following R15, activity variability in the RV, log 𝑅 (cid:48) HK and BIStime series can be modelled as: Δ RV = 𝑉 𝑐 𝐺 ( 𝑡 ) + 𝑉 𝑟 (cid:164) 𝐺 ( 𝑡 ) , (13)log 𝑅 (cid:48) HK = 𝐿 𝑐 𝐺 ( 𝑡 ) , and (14)BIS = 𝐵 𝑐 𝐺 ( 𝑡 ) + 𝐵 𝑟 (cid:164) 𝐺 ( 𝑡 ) , (15)respectively. The coefficients 𝑉 𝑐 , 𝑉 𝑟 , 𝐿 𝑐 , 𝐵 𝑐 and 𝐵 𝑟 are free param-eters relating the individual observations to the unobserved Gaus-sian process 𝐺 ( 𝑡 ) . In R15’s framework, 𝐺 ( 𝑡 ) itself can be looselyinterpreted as representing the projected area of the visible stellardisc covered in active regions at a given time; the GP describing 𝐺 ( 𝑡 ) is assumed to have zero mean and covariance matrix K , where 𝐾 𝑖 𝑗 = 𝛾 ( 𝑡 𝑖 , 𝑡 𝑗 ) . As in R15, we adopt the following quasi-periodiccovariance kernel function: 𝛾 ( 𝑡 𝑖 , 𝑡 𝑗 ) = exp (cid:34) − sin (cid:2) π ( 𝑡 𝑖 − 𝑡 𝑗 )/ 𝑃 GP (cid:3) 𝜆 − ( 𝑡 𝑖 − 𝑡 𝑗 ) 𝜆 (cid:35) , (16)where 𝑃 GP is the period of the quasi-periodic activity signal, 𝜆 p is the inverse harmonic complexity of the signal (such that signalsbecome sinusoidal for large values of 𝜆 p , and show increasing com-plexity/harmonic content for small values of 𝜆 p ), and 𝜆 e is the timescale over which activity signals evolve. This quasi-periodic covari-ance kernel has been widely used to model stellar activity signalsin both photometry and RVs (e.g. Haywood et al. 2014a; Rajpaulet al. 2015; Grunblatt et al. 2015; Bonfils et al. 2018). The full ex-pressions for the covariance between the three different observablesmodelled are given in R15.We expect, in principle at least, that this sophisticated GP-based approach to modelling RVs jointly with activity indicatorsshould enable more reliable planet characterization, for several rea-sons. Firstly, by modelling multiple activity-sensitive time seriessimultaneously (e.g. log 𝑅 (cid:48) HK , BIS, FWHM, etc., or some subset ofthese), more information can be gleaned on activity signals in RVs,compared to exploiting only simple correlations between RVs and(typically) one of these time series. Additionally, the frameworkuses GP draws and derivatives thereof as basis functions for mod-elling available time series, rather than e.g. sinusoids or other simpleparametric models, the inappropriate use of which could easily leadto the introduction of correlated signals into model residuals. TheGP basis functions could in principle take any form, although inthe GP framework their properties are constrained to some extentby the data itself, and from reasonable prior assumptions about thequasi-periodic nature of stellar activity signals. The GP frameworkalso incorporates the 𝐹𝐹 (cid:48) formalism directly as a special case; theformer approach may be thought of as a generalization of the latter. As we shall use Bayesian inference to evaluate the relative posteriorprobabilities of different models, we summarise briefly here therelevant formalism. Firstly, Bayes’ Theorem, given in equation (17),is used to relate (i) the posterior probability Pr ( Θ | 𝐷, 𝑀 ) = P ( Θ ) of the parameters Θ given data 𝐷 and a model 𝑀 to (ii) the priordistribution Pr ( Θ | 𝑀 ) = 𝜋 ( Θ ) of Θ given 𝑀 , (iii) the likelihoodPr ( 𝐷 | Θ , 𝑀 ) = L( Θ ) of 𝐷 given Θ and 𝑀 and the Bayesian evidencePr ( 𝐷 | 𝑀 ) = Z of 𝐷 given 𝑀 . Following the notation used by Ferozet al. (2009):Pr ( Θ | 𝐷, 𝑀 ) = Pr ( 𝐷 | Θ , 𝑀 ) Pr ( Θ | 𝑀 ) Pr ( 𝐷 | 𝑀 ) , (17)or simply P ( Θ ) = L( Θ ) 𝜋 ( Θ )Z . (18)For model selection, one can compare two models 𝑀 and 𝑀 ,given data 𝐷 , by computing the ratio of their respective posteriorprobabilities; this ratio is also known as the Bayes factor 𝑅 : R = Pr ( 𝑀 | 𝐷 ) Pr ( 𝑀 | 𝐷 ) = Pr ( 𝐷 | 𝑀 ) Pr ( 𝑀 ) Pr ( 𝐷 | 𝑀 ) Pr ( 𝑀 ) = Z Z Pr ( 𝑀 ) Pr ( 𝑀 ) . (19) MNRAS000
40 d based on the average of log 𝑅 (cid:48) HK measurements by Loviset al. (2011); (ii) signals likely corresponding to the harmonics ofthe rotation period are detected in the periodograms of the stellaractivity indicators, see Section 6.1 (Fig. 3). Further motivation wasprovided by the fact that (iii) when considering zero or 1-planetmodels, the rotation period of the star became locked on to theperiod of one of the two planets, and that (iv) during preliminaryruns, the MAP value for the rotation period never went above 40 d.By limiting the period to below 42 d we reduced the computationalburden of our modelling. In a similar fashion to the harmonic activity model, the ‘magn.cycle’ model fits a sinusoid to two data sets, this time to the RVs andthe log 𝑅 (cid:48) HK time series, enforcing an identical period but allowingdifferent phases and amplitudes. We assume that this captures stellaractivity corresponding to a long-term magnetic activity cycle.Equation (7) and (8) show the two models with the cycle period 𝑃 magn simultaneously fit to the RV and log 𝑅 (cid:48) HK measurements,respectively. The parameters 𝐾 magn , RV , 𝜙 magn , RV and 𝐾 magn , 𝜙 magn and 𝐶 magn describe the RV ( 𝑡 ) and log R (cid:48) HK ( 𝑡 ) model respectively; 𝑉 represents a constant offset and a quadratic polynomial term.RV total ( 𝑡 ) = RV Kepler ( 𝑡 ) + 𝐾 magn , RV sin (cid:18) π 𝑡𝑃 magn + 𝜙 magn , RV (cid:19) + 𝑉 (7)log R (cid:48) HK ( 𝑡 ) = 𝐾 magn sin (cid:18) π 𝑡𝑃 magn + 𝜙 magn (cid:19) + 𝐶 magn (8)The prior distribution for 𝐾 is as in Table 3; the priors for 𝜙, 𝐶 magn and 𝑃 magn are displayed in Table 4. The upper and lowerprior limits for 𝑃 magn were determined by a broad peak around1000—4000 d in the periodogram of log 𝑅 (cid:48) HK as well as by theobvious periodic signal in the log 𝑅 (cid:48) HK measurements (see Fig. 2). 𝐹𝐹 (cid:48) method The FWHM measurements were used for applying the 𝐹𝐹 (cid:48) methodwhich was introduced by Aigrain et al. (2012) as a method for re-lating the photometric brightness and RV variations of a star. Ituses the flux of the star Ψ ( 𝑡 ) and its derivative (cid:164) Ψ ( 𝑡 ) as an indica-tor for spot coverage and predicts RV variations. These variationsinclude the RV perturbation Δ RV rot ( t ) due to the presence of spots MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 on the rotating photosphere and their effect of the suppression ofconvective blueshift Δ RV conv ( t ) .Aigrain et al. (2012) describes these with the following twoequations where 𝑓 represents the drop in flux produced by a spotat the centre of the stellar disc, 𝑅 ∗ is the stellar radius, 𝛿𝑉 𝑐 isthe difference between the convective blueshift in the unspottedphotosphere and that within the magnetized area and 𝜅 is the ratioof this area to the spot surface. Δ RV rot ( t ) = − Ψ ( t ) Ψ (cid:20) − Ψ ( t ) Ψ (cid:21) R ∗ f (9) Δ RV conv ( t ) = (cid:20) − Ψ ( t ) Ψ (cid:21) 𝛿 V c 𝜅 f (10)The total RV variation Δ RV activity created by stellar activity isthen the sum of both terms: Δ RV activity = Δ RV rot ( t ) + Δ RV conv ( t ) (11)The parameter 𝑓 is approximated as 𝑓 ≈ Ψ − Φ min Ψ (12)with Φ min being the minimum observed flux.It has been argued that CCF FWHM andlog 𝑅 (cid:48) HK measurements should both behave, to first order, asthe convective blueshift suppression term in the 𝐹𝐹 (cid:48) method,which in turn is a close proxy for the integrated active regioncoverage on the visible stellar hemisphere (Rajpaul et al. 2015);indeed, it has been shown in practice that the CCF FWHM canbe a good tracer of photometric flux (e.g. Suárez Mascareñoet al. 2020). Thus, as we did not have photometric measurementsof HD 13808 available, we chose to interpret the FWHM as aproxy for the flux Ψ ( 𝑡 ) . The derivative of the ‘flux’ was computednumerically by modelling the FWHM with a GP using the celeritealgorithm by Foreman-Mackey et al. (2017). In total, there werethree free parameters in this stellar activity model, 𝛿𝑉 𝑐 𝜅 , 𝑅 ∗ and Ψ - their respective priors appear in Table 4. We used the GP framework developed by Rajpaul et al. (2015),hereafter R15, to model RVs simultaneously with log 𝑅 (cid:48) HK and BISobservations. In short, this framework assumes that all observed stel-lar activity signals are generated by some underlying latent function 𝐺 ( 𝑡 ) and its derivatives; this function, which is not observed directly,is modelled with a Gaussian process (Rasmussen & Williams 2006;Roberts et al. 2013).Following R15, activity variability in the RV, log 𝑅 (cid:48) HK and BIStime series can be modelled as: Δ RV = 𝑉 𝑐 𝐺 ( 𝑡 ) + 𝑉 𝑟 (cid:164) 𝐺 ( 𝑡 ) , (13)log 𝑅 (cid:48) HK = 𝐿 𝑐 𝐺 ( 𝑡 ) , and (14)BIS = 𝐵 𝑐 𝐺 ( 𝑡 ) + 𝐵 𝑟 (cid:164) 𝐺 ( 𝑡 ) , (15)respectively. The coefficients 𝑉 𝑐 , 𝑉 𝑟 , 𝐿 𝑐 , 𝐵 𝑐 and 𝐵 𝑟 are free param-eters relating the individual observations to the unobserved Gaus-sian process 𝐺 ( 𝑡 ) . In R15’s framework, 𝐺 ( 𝑡 ) itself can be looselyinterpreted as representing the projected area of the visible stellardisc covered in active regions at a given time; the GP describing 𝐺 ( 𝑡 ) is assumed to have zero mean and covariance matrix K , where 𝐾 𝑖 𝑗 = 𝛾 ( 𝑡 𝑖 , 𝑡 𝑗 ) . As in R15, we adopt the following quasi-periodiccovariance kernel function: 𝛾 ( 𝑡 𝑖 , 𝑡 𝑗 ) = exp (cid:34) − sin (cid:2) π ( 𝑡 𝑖 − 𝑡 𝑗 )/ 𝑃 GP (cid:3) 𝜆 − ( 𝑡 𝑖 − 𝑡 𝑗 ) 𝜆 (cid:35) , (16)where 𝑃 GP is the period of the quasi-periodic activity signal, 𝜆 p is the inverse harmonic complexity of the signal (such that signalsbecome sinusoidal for large values of 𝜆 p , and show increasing com-plexity/harmonic content for small values of 𝜆 p ), and 𝜆 e is the timescale over which activity signals evolve. This quasi-periodic covari-ance kernel has been widely used to model stellar activity signalsin both photometry and RVs (e.g. Haywood et al. 2014a; Rajpaulet al. 2015; Grunblatt et al. 2015; Bonfils et al. 2018). The full ex-pressions for the covariance between the three different observablesmodelled are given in R15.We expect, in principle at least, that this sophisticated GP-based approach to modelling RVs jointly with activity indicatorsshould enable more reliable planet characterization, for several rea-sons. Firstly, by modelling multiple activity-sensitive time seriessimultaneously (e.g. log 𝑅 (cid:48) HK , BIS, FWHM, etc., or some subset ofthese), more information can be gleaned on activity signals in RVs,compared to exploiting only simple correlations between RVs and(typically) one of these time series. Additionally, the frameworkuses GP draws and derivatives thereof as basis functions for mod-elling available time series, rather than e.g. sinusoids or other simpleparametric models, the inappropriate use of which could easily leadto the introduction of correlated signals into model residuals. TheGP basis functions could in principle take any form, although inthe GP framework their properties are constrained to some extentby the data itself, and from reasonable prior assumptions about thequasi-periodic nature of stellar activity signals. The GP frameworkalso incorporates the 𝐹𝐹 (cid:48) formalism directly as a special case; theformer approach may be thought of as a generalization of the latter. As we shall use Bayesian inference to evaluate the relative posteriorprobabilities of different models, we summarise briefly here therelevant formalism. Firstly, Bayes’ Theorem, given in equation (17),is used to relate (i) the posterior probability Pr ( Θ | 𝐷, 𝑀 ) = P ( Θ ) of the parameters Θ given data 𝐷 and a model 𝑀 to (ii) the priordistribution Pr ( Θ | 𝑀 ) = 𝜋 ( Θ ) of Θ given 𝑀 , (iii) the likelihoodPr ( 𝐷 | Θ , 𝑀 ) = L( Θ ) of 𝐷 given Θ and 𝑀 and the Bayesian evidencePr ( 𝐷 | 𝑀 ) = Z of 𝐷 given 𝑀 . Following the notation used by Ferozet al. (2009):Pr ( Θ | 𝐷, 𝑀 ) = Pr ( 𝐷 | Θ , 𝑀 ) Pr ( Θ | 𝑀 ) Pr ( 𝐷 | 𝑀 ) , (17)or simply P ( Θ ) = L( Θ ) 𝜋 ( Θ )Z . (18)For model selection, one can compare two models 𝑀 and 𝑀 ,given data 𝐷 , by computing the ratio of their respective posteriorprobabilities; this ratio is also known as the Bayes factor 𝑅 : R = Pr ( 𝑀 | 𝐷 ) Pr ( 𝑀 | 𝐷 ) = Pr ( 𝐷 | 𝑀 ) Pr ( 𝑀 ) Pr ( 𝐷 | 𝑀 ) Pr ( 𝑀 ) = Z Z Pr ( 𝑀 ) Pr ( 𝑀 ) . (19) MNRAS000 , 1–17 (2020)
E. Ahrer et al.
Table 5.
Jeffreys’ scale, as introduced by Jeffreys (1983), for interpreting dif-ferences in Bayesian evidences. This scale interprets the strength of evidencefavouring one model over another.|ln R | Odds Probability Remark< 1.0 (cid:46) ∼ ∼ ∼ Note that Pr ( 𝑀 ) Pr ( 𝑀 ) is the relative a priori probability between the twomodels, which is usually set to one.To decide whether the evidence difference is significant to favorone model over the other, we make use of the Jeffreys scale as givenin Table 5, where the logarithmic scale of the evidence difference isused. This simplifies equation (19) to:ln (R) = ln (cid:18) Z Z (cid:19) = ln Z − ln Z (20) It is commonly assumed that an additive white Gaussian noise(AWGN) model is sufficient for describing observational noise, asin probability theory the central limit theorem states that the sum ofindependent random variables tends towards a Gaussian distribu-tion, even if the original ones are not normally distributed (Fischer2011). In this case, the likelihood L( Θ ) of parameters Θ can bewritten as: L( Θ ) = 𝑁 (cid:214) 𝑖 = √︃ π 𝜎 𝑖 exp (cid:32) − [ 𝑣 ( 𝑡 𝑖 ; Θ ) − 𝑣 𝑖 ] 𝜎 𝑖 (cid:33) , (21)where 𝑣 ( 𝑡 𝑖 ; Θ ) describes the model’s predicted RV for parameters Θ ,while 𝑣 𝑖 describes the RV observed at time 𝑡 𝑖 , with correspondingerror estimate 𝜎 𝑖 which contains the observational error and theadditive noise ‘jitter’ term 𝜎 + RV added in quadrature. The logarithmiclikelihood, which is more convenient to work with, is computed as:ln L( Θ ) = 𝑁 ∑︁ 𝑖 = − ln √︃ π 𝜎 𝑖 − 𝜎 𝑖 [ 𝑣 ( 𝑡 𝑖 ; Θ ) − 𝑣 𝑖 ] . (22)Note that the above formalism applies to all of our models except for the GP model, which explicitly generalises the AWGN model byallowing for red (correlated) noise. In this case, the log likelihoodmay be computed via the more general expressionln L( Θ ) = − 𝑁 ln 2 𝜋 − ln det K − r T K − r , (23)where K is the matrix defining the covariance between all pairs ofobservations (see R15), and r is a vector of residuals with the 𝑖 th element given by 𝑣 ( 𝑡 𝑖 ; Θ ) − 𝑣 𝑖 . Note that in the special case where K is a diagonal matrix with the 𝑖 th element given by 𝜎 𝑖 , i.e. wherethe noise is assumed to be white, equation 23 reduces to equation22. As the Keplerian model requires the Kepler equation to besolved, we made use of the efficient CORDIC-like method intro-duced by Zechmeister (2018) where double precision is obtainedwithin 55 iterations. While computing likelihoods we checked if each pair of planetswould be stable following the criterion introduced by Gladman(1993): Δ > √ 𝑅 𝐻 ( 𝑖, 𝑗 ) where Δ = 𝑎 𝑗 − 𝑎 𝑖 is the differencebetween the semi-major axis of the 𝑖 -th and 𝑗 -th planet and 𝑅 𝐻 ( 𝑖, 𝑗 ) the planets’ mutual Hill radius. This has been used before e.g. byMalavolta et al. (2017).If the parameters drawn resulted in a planet system which wasconsidered unstable, the likelihood was set to zero, i.e. 𝑍 =
0; inour case, as we worked with log likelihoods, this was approximatednumerically by ln 𝑍 = − × . We used PolyChord as it is an state-of-the-art nested samplingalgorithm, designed to work with very high dimensional parameterspaces (Handley et al. 2015). A short discussion comparing it withMultiNest (Feroz et al. 2009) can be found in Hall et al. (2018). Itwas developed using C++ and Fortran and can also be called as aPython package.In order to test the reliability of our algorithm with PolyChordwe used two sets of simulated data. One data set contained randomuncorrelated Gaussian noise (with a standard deviation of 1m s − )without any planetary signals, while the other one consisted oftwo planets on elliptical orbits with the same Gaussian noise. Bothdata sets also included an offset of a few m s − and linear trendof 10 − m s − per day. Models of 0–3 planet signals with longterm trends in the form of second-order polynomials were fittedand their evidences compared. Both elliptical and circular orbitalsolutions were computed. We note here that in general, the evidenceuncertainty reported by PolyChord is an underestimate, and it isbetter practice to estimate it via the range of scatter in the evidenceacross multiple runs (Higson et al. 2018; Nelson et al. 2020); weused the latter approach throughout our analyses. (It is probable thatthis behaviour would also be ameliorated by increasing the numberof live points beyond the computational resources available at thetime of writing.)For the first data set without planetary signals, the algorithmsuccessfully favoured the no-planet model consisting of a second-order polynomial. Every other model computed showed lower ev-idences, with at least ln R ≈
R ≈ ∼ ∼ for these sort of exoplanet applications. This resulted ininconsistencies where the sampling in some runs missed transitionsto areas of higher likelihood (see Fig. 1) and showed very differentevidence values and posterior distributions. To ensure consistent The convergence criterion in PolyChord is defined as the point wherethe fraction of total evidence contained in the live points drops below thedefault value of 10 − (Handley et al. 2015). MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 log l o g X e b cos( b )1 0 1 e b cos( b )1 0 1 e b cos( b )0246 K b m s Figure 1.
A typical run with PolyChord. The bottom graph shows thelogarithmic likelihood over the course of the run, with log 𝑋 correspondingto the prior volume i.e. at the start of the run log 𝑋 is at its highest value. Twophase transitions are clearly visible as ’knees’ in the logL-logX curve. On thetop, the live-point distributions of two parameters ( 𝐾 𝑏 and √ 𝑒 𝑏 cos ( 𝜔 𝑏 ) )within a given likelihood contour (indicated by the red vertical lines) areshown, demonstrating the transitions between different phases. Note that ifthe precision criterion is not set low enough, one runs the risk of stoppingnested sampling at too early a log 𝑋 value (i.e. at one of the lower knees),and this is indeed what we observed in preliminary tests. Table 6.
An overview of the HARPS RV data subsets of HD 13808 re-ferred to as low and high activity when log 𝑅 (cid:48) HK is < − .
90 and > − . 𝑅 (cid:48) HK < − .
90 > − .
90N data 123 123 Δ T (d) 3332 3701rms(RV) (m s − ) 3.53 3.78mean( 𝜎 RV ) (m s − ) 0.72 0.78mean(log 𝑅 (cid:48) HK ) − . − . and robust runs, the convergence criterion was lowered to valuesbetween 10 − to 10 − for runs with a large number of planetsand complex stellar activity models. As an historical aside, nestedsampling was in fact invented (Skilling 2006) in order to solveprecisely these kind of phase transition problems. The extent of stellar activity observed from HD 13808 was investi-gated by looking at the activity indicators BIS, FWHM and log 𝑅 (cid:48) HK .None of them showed an obvious linear correlation with the RVmeasurements. Note that this can occur e.g. when there is a relationbut with a simple time lag, as demonstrated with the case of the Sun Figure 2.
Measurement of the log 𝑅 (cid:48) HK activity index on the series of spectraof HD 13808. The black dashed line represents the median value of − . Collier Cameron et al. (2019). However, the log 𝑅 (cid:48) HK measurementsdo show a notable cycle from a low of − .
10 to a high of − .
70, asdisplayed in Fig. 2. A cycle with similar shape and period is alsoseen in the BIS and FWHM time series, although with smaller am-plitude. We thus decided to extend our analysis by splitting the fullRV data into two subsets based on the median(log 𝑅 (cid:48) HK ) = − . 𝑅 (cid:48) HK < − .
90 and log 𝑅 (cid:48) HK > − .
90, respectively. Table 6summarises the statistics of the full data set and those two subsets.Note that although we call one subset ‘high activity,’ HD 13808 isstill classified as an inactive star, as an ‘active’ classification wouldusually require an average of log 𝑅 (cid:48) HK > − .
75 (Henry et al. 1996).This log 𝑅 (cid:48) HK cycle shows great similarities to the one of the Sunwith its recent maxima and minima being at log 𝑅 (cid:48) HK (cid:39) − .
83 andlog 𝑅 (cid:48) HK (cid:39) − .
96; the extrapolated log 𝑅 (cid:48) HK value for the Maun-der minimum period is log 𝑅 (cid:48) HK (cid:39) − .
10 (Mamajek & Hillenbrand2008).
To identify strong periodicities in our time series we computed theBayesian generalized Lomb-Scargle (BGLS) periodogram (Lomb1976; Scargle 1982; Zechmeister & Kürster 2009; Mortier et al.2015) of the two separated data subsets for the RV, the stellar activityindicators BIS, FWHM and log 𝑅 (cid:48) HK and a constant function (CF)with identical time stamps.The resulting BGLS periodograms are shown in Fig. 3 wherepeaks in the RV periodograms are visible at ∼
15 and ∼
55 d. Setsof strong peaks are visible in the BGLS periodograms of the activityindicators log 𝑅 (cid:48) HK , BIS and FWHM between 30 −
40 d and a fewbetween 25 −
30 d. Strong peaks are evident at ∼
19 and ∼
32 d inthe log 𝑅 (cid:48) HK high activity periodogram, with peaks at ∼ ∼ ∼
34 d in the BIS periodogram. The FWHM periodogramsshow peaks in the region of ∼
35 d.The BGLS periodograms of the high- and low-activity subsetsof the data suggest that the RV measurements are affected by stellaractivity: note, in particular, the excess power around the putativestellar rotation period and the first two harmonics in the high-activitysubsets of the BIS and log 𝑅 (cid:48) HK time series. It is also evident thatthe stellar activity indicators differ quite strongly when the star isslightly more active; the difference between the mean log 𝑅 (cid:48) HK valueof the two subsets is Δ (cid:10) log 𝑅 (cid:48) HK (cid:11) = . MNRAS , 1–17 (2020)
E. Ahrer et al.
Figure 3.
BGLS periodograms of, from top to bottom, the RV, BIS, FWHM and log 𝑅 (cid:48) HK measurements and of a constant function (CF) i.e. a dataset with thesame timestamps and RV errors but assuming constant RV. The colours red and black indicate the high activity and low activity phases. The dark grey shadedbox represents the range of the putative rotation period with its second (grey) and third harmonic (light grey). The different models we considered are summarised in Table 7,along with a short description of each. The model including onlycircular orbital solutions is called ‘Circular’, while the ellipticalorbital solutions are referred to as ‘Kepler’. A ‘+’ indicates that the‘Kepler’ solution was complemented with a specific stellar activitymodel, such as linear dependency on the BIS (‘Kepler + linearBIS’). This was done equivalently for a subset of models applied tothe low and high-activity subsets of the data.
The relative Bayesian evidences for each model with different num-ber of Keplerians are shown in Table 8. For ease of comparison, weset the log evidence for models with two Keplerians to zero. Thuspositive values are more favoured and negative ones less favouredthan the 2-planet models, with significance to be interpreted accord-ing to Jeffreys’ scale (Table 5). Table 8 also includes the rms andmedian values for the RV residuals computed using the maximuma posteriori (MAP) values of the various 2-planet models.Note that the absolute evidences of different models cannotbe meaningfully compared in a global sense as the evidence valuesdepend, among other things, on the data fitted; hence an 𝑛 -Keplerian MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 Table 7.
HD 13808: Overview of the models used for the study of the HD 13808.Model name DescriptionCircular Circular orbit modellingKepler Eccentric orbit modelling+ magn. cycle Kepler + long-term magnetic activity cycle; Section 4.3+ linear BIS Kepler + linear dependency on the BIS; Section 4.1, equation (3)+ BIS 𝑃 rot st harm. Kepler + simultaneous fit of 1 st harmonic of 𝑃 rot to the BIS and RVs; Section 4.2+ BIS 𝑃 rot nd harm. Kepler + simultaneous fit of 1 st and 2 nd harmonic of 𝑃 rot to the BIS and RVs; Section 4.2+ BIS 𝑃 rot rd harm. Kepler + simultaneous fit of 1 st , 2 nd and 3 rd harmonic of 𝑃 rot to the BIS and RVs; Section 4.2+ linear log 𝑅 (cid:48) HK Kepler + linear dependency on the log 𝑅 (cid:48) HK ; Section 4.1, equation (4)+ 𝐹 𝐹 (cid:48)
Kepler + FF’ method with FWHM as a proxy for flux; Section 4.4Gaussian process Simultaneous GP modelling of RVs and activity indicators; Keplerian terms for RVs; Section 4.5
Table 8.
HD 13808: Relative Bayesian evidences for different number of planets for each model and their scatter – as well as the residual RV (O − C) rms andmedian absolute deviation (MAD) – for the 2-planet model (or 1-planet model, if favoured) computed using the respective MAP values (see Table 10). Themodel with two planetary signals for each model type was set to be zero. The errors correspond to the standard error on the mean evidence across multipleruns or to the highest individual run error provided by PolyChord, whichever of the two is greater.Model No planets 1 planet 2 planets 3 planets 4 planets RV residual RV residualln (R) ln (R) ln (R) ln (R) ln (R) rms m s − MAD m s − Circular − ± − ± ± + ± + ± − ± − ± ± + ± + ± − ± − ± ± + ± − ±
13 3.89 2.41+ linear BIS − ± − ± ± − ± ± 𝑃 rot st harm. − ± − ± ± − ± − ± 𝑃 rot nd harm. − ± − ± ± − ± − ± 𝑃 rot rd harm. − ± − ± ± − ± − ± − ± − ± ± − ± + ± 𝑅 (cid:48) HK − ± − ± ± − ± − ± 𝐹 𝐹 (cid:48) − ± − ± ± − ± − ± − ± − ± ± − ± − ± − ± − ± ± − ± + ± − ± − ± ± + ± + ± 𝑃 rot rd harm. − ± − ± ± + ± + ± 𝐹 𝐹 (cid:48) − ± − ± ± − ± + ± − ± − ± ± − ± − ± − ± − ± ± + ± + ± − ± − ± ± + ± + ± 𝑃 rot rd harm. − ± + ± ± + ± + ± 𝐹 𝐹 (cid:48) − ± − ± ± − ± − ± − ± − ± ± − ± − ±000
13 3.89 2.41+ linear BIS − ± − ± ± − ± ± 𝑃 rot st harm. − ± − ± ± − ± − ± 𝑃 rot nd harm. − ± − ± ± − ± − ± 𝑃 rot rd harm. − ± − ± ± − ± − ± − ± − ± ± − ± + ± 𝑅 (cid:48) HK − ± − ± ± − ± − ± 𝐹 𝐹 (cid:48) − ± − ± ± − ± − ± − ± − ± ± − ± − ± − ± − ± ± − ± + ± − ± − ± ± + ± + ± 𝑃 rot rd harm. − ± − ± ± + ± + ± 𝐹 𝐹 (cid:48) − ± − ± ± − ± + ± − ± − ± ± − ± − ± − ± − ± ± + ± + ± − ± − ± ± + ± + ± 𝑃 rot rd harm. − ± + ± ± + ± + ± 𝐹 𝐹 (cid:48) − ± − ± ± − ± − ± − ± − ± ± − ± − ±000 , 1–17 (2020) E. Ahrer et al. model that fits only RVs cannot e.g. be meaningfully compared toanother 𝑛 -Keplerian model that fits RV simultaneously with a BISor log 𝑅 (cid:48) HK time series. In passing we note, however, that it was possible to compare the absolute evidence values of the parametricmodel ‘Kepler + BIS 𝑃 rot rd harm. + magn. cycle’ with the GPmodel, as these fitted the same three time series; the 2-planet GPmodel achieved a log evidence an order of magnitude higher thanthe former model ( − ± − ± 𝑅 (cid:48) HK , andBIS series jointly), and moreover that high-quality fits under theGP model were localised to a comparatively small volume of priorspace. Table 8 shows that many of the simple parametric models – thoughnot the GP model – favoured the ‘detection’ of three or even fourplanets. However, we have strong empirical reasons to reject modelswith more than two Keplerian components.First, we found that the periods of the third Keplerian signal(and also the fourth, where applicable), tended to change with everyrun, while the periods of ∼
14 d and ∼
54 d appeared as MAP valuesin every run, for all models with two or more Keplerians. This isdemonstrated in Table 9, where we list the periods of the additionalKeplerian periods found by every model. The most common peri-ods for an ‘extra’ Keplerian were around 12 and 19 d – as it turnsout, these periods out are also strongly visible in the high-activityperiodogram of the log 𝑅 (cid:48) HK and BIS time series (see Fig. 3). Asmentioned in the BGLS periodogram analysis, we interpret theseperiods as harmonics of the stellar rotation period. The fact thatthese periodicities manifested in Keplerian terms in all but the mostcomplex activity models reflects the inadequacy of the simpler mod-els at capturing all the activity variability – neither the 12 d nor the19 d Keplerian periodicities showed up under the GP model.Second, it also turned out that these supposed planetary sig-nals described by the Keplerians either did not show up or theirsemi-amplitudes decrease to below 1 m s − when consideringonly the low-activity subset of observations. Conversely, the semi-amplitudes of the supposed third planets occurring at ∼
12 d or ∼
19 d increased significantly when only the high activity data setwas modelled compared to when using the full data set (see Table 9).Perhaps most significantly, our GP model – which we regarded a priori as being by far the most realistic of our activity models –decisively favoured a 2-planet interpretation of our full set of obser-vations. The GP modelling even favoured a 2-planet interpretationwhen considering only the high-activity subset of observations, andresulted in planet properties consistent with those derived from thelow-activity subset or the full data set (see Table 10). This was nottrue for the simpler parametric models; the 𝐹𝐹 (cid:48) method equivo-cated between 2- and 3-planet solutions, and was not reliably ableto detect the ∼
54 d periodicity.Our GP activity model, though relatively complex in its ownright, did a comparatively good job of fitting variability in RVsunrelated to planets, when compared to most of the other activitymodels (see residuals in Table 8). This is despite our requirementthat the GP fit activity-related variability in RVs simultaneously with BIS and log 𝑅 (cid:48) HK time series, which ensured that the GP wouldbe very unlikely to try to fit planetary variability. Presumably, then,under the GP modelling there was little residual activity-inducedvariability that could be explained using a third or fourth Keplerianterm; hence these more complex models were rejected. Regarding the parametric activity models (i.e. not the GPmodel), the addition of extra Keplerian terms beyond the 2-planetmodel in many cases did enable a non-trivial fraction of extraactivity-related variability to be fitted. This may have been the case,for example, because the amplitudes, phases and periods of thedifferent terms in the harmonic models were by definition fixed,whereas the quasi-periodic GP model is flexible enough to fit bothlow and high activity variability, possibly with changing phases andarising from a combination of rotation periods (as might be thecase due e.g. to differential rotation), rather than a single ‘master’rotation period (Rajpaul 2017). Thus the non-GP activity models inmany cases favoured solutions with three or four Keplerian terms.The upshot, then, is that one’s conclusions about the number ofplanetary signals present in an RV time series is extremely sensitiveto whether one has modelled other, non-planetary (e.g. activity-induced) variability present in the RVs. Given the high dimensionality of the joint stellar activity plus plan-etary models we considered (over 30 parameters for some of the4-planet plus activity models), accurate Bayesian evidence compu-tation required a large number of posterior samples to be drawn –in the case of 4-planet models, several hundred million posteriorsamples were typically required, even with a state-of-the-art sam-pler such as PolyChord. This did not prove too burdensome for theparametric activity models; a typical run for a 4-planet model with ∼ − . × The upshot was thatevaluating our GP models required of order fifty thousand CPUcore hours on a high-performance computing platform. (Attemptsto speed up computation by reducing PolyChord’s precision cri-terion led to unacceptably large scatter in computed evidences.)Accurate evaluation of the Bayesian evidences of these GP modelswould simply not have been feasible on a desktop or even a smallcomputing cluster.At face value this might seem to be a serious shortcoming of theGP model. However, we note that Bayesian evidence computation iscomputationally challenging in general (Nelson et al. 2020), and itwould seem that you ‘get what you pay for’: the parametric activitymodels are certainly far cheaper to evaluate, though as we argued inSection 6.2.2, their over-simplicity inevitably leads to the detectionof spurious planets. Such issues are avoided with the GP model.
The posterior distributions for the inferred semi-amplitude 𝐾 , pe-riod 𝑃 and eccentricity 𝑒 of the planets in all of our 2-planet modelsare summarised in Table 10. Despite the wide array of modelsconsidered, the planetary parameters show a remarkable degree of While a number of techniques for speeding up GP regression do exist, e.g.Foreman-Mackey et al. (2017), we are not aware of any that are well-suitedto cases where covariances need to be formulated over multiple inputs andoutputs simultaneously, e.g. RVs jointly with activity indicators.MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 Table 9.
HD 13808: Periods for a putative third planet found in the various runs by the different models, their corresponding semi-amplitudes and the uncertainty(per Table 8) in the model evidences. Values in bold occur in at least 2 of the runs.Model 𝑁 runs 𝜎 ln Z Periods (d) Semi-amplitudes (m s − )Circular 5 3 , , ∼ . , . , . , , ∼ . , . , .
8+ magn. cycle 3 4 12 , ∼ . , .
8+ linear BIS 3 2 10 , , ∼ . , . , .
9+ BIS 𝑃 rot st harm. 3 7 11 , , ∼ . , . , .
7+ BIS 𝑃 rot nd harm. 3 6 10 , , ∼ . , . , .
2+ BIS 𝑃 rot rd harm. 3 4 , ∼ . , .
9+ magn. cycle 3 8 , ∼ . , .
7+ linear log 𝑅 (cid:48) HK , , ∼ . , . , . 𝐹 𝐹 (cid:48) , ,
43 1 . , . , . , ∼ . , . , , , ∼ . , . , . , . , , ∼ . , . , .
0+ BIS 𝑃 rot rd harm. 3 4 8 , , ∼ . , . , . 𝐹 𝐹 (cid:48) , , ∼ . , . , . , , ∼ . , . , . , , ∼ . , . , . , ∼ . , .
6+ BIS 𝑃 rot rd harm. 3 3 9 , , ∼ . , . , . 𝐹 𝐹 (cid:48) , ∼ . , . , ∼ . , . consistency, with the credible intervals for their semi-amplitudesand eccentricities agreeing within 1 𝜎 across almost all models; thesame holds true when considering only the low-activity subset ofobservations. Differences between orbital periods inferred underdifferent models are never greater than a fraction of a per cent. Thisprovide additional evidence that the detected signals are indeedplanetary and robust. It is worth noting that the GP model favours amarginally larger MAP semi-amplitude 𝐾 𝑏 , and a marginally largerMAP semi-amplitude 𝐾 𝑐 , than the typical parametric models. How-ever, the values favoured by the GP model are bracketed both aboveand below by other parametric models, and so are in no sense ex-treme; concern that a GP might ‘absorb’ some of a planetary signalis clearly unfounded in this case. Under the parametric activitymodels, 𝐾 𝑐 tends to decrease in the low-activity data subset com-pared to the high-activity subset, suggesting that the outer Keplerian(with period broadly similar to the likely stellar rotation period) isactually being used to absorb some RV variability the too-simplisticactivity models cannot account for. 𝐾 𝑏 appears comparatively in- After all, the framework from R15 is designed so that the GP componentof the model is only able to explain variability simultaneously present in RVsand activity-sensitive time series; therefore, except in extremely pathologicalcases, there should be little risk of the GP wrongly fitting a planetary signal. sensitive to stellar activity levels, which may reflect the fact that 𝑃 𝑏 is about 2 . > 𝜎 levels. On this basis we regardour modelling to represent the secure RV detection of two planetsorbiting HD 13808. By way of contrast, it is worth noting thateven transiting planets with tightly constrained periods and orbitalphases may sometimes prove difficult to characterize, with inferredsemi-amplitudes diverging wildly depending on which model isused: Kepler-10 c is an excellent example (Fressin et al. 2011;Dumusque et al. 2014; Weiss et al. 2016; Rajpaul et al. 2017). Werefer, hereafter, to the planets we have detected as HD 13808 b andHD 13808 c with periods of ∼
14 and ∼
54 d, respectively.
The stellar rotation period of 38 . ± .
47 d inferred under theGP model is broadly bracketed by the ∼ MNRAS000
47 d inferred under theGP model is broadly bracketed by the ∼ MNRAS000 , 1–17 (2020) E. Ahrer et al.
Table 10.
HD 13808: Comparison of the orbital parameters semi-amplitude 𝐾 , period 𝑃 and eccentricity 𝑒 for the two planets HD 13808 b and HD 13808 cfor each model. The parameters for the 2-planet models of the run with the highest evidence from all runs are shown here. Note that not all models favour the2-planet case.Model 𝐾 𝑏 (m s − ) 𝐾 𝑐 (m s − ) 𝑃 𝑏 (d) 𝑃 𝑐 (d) 𝑒 𝑏 𝑒 𝑐 Circular 3 . + . − . . + . − . . + . − . . + . − . - -Kepler 3 . ± .
25 2 . ± .
28 14 . + . − . . + . − . . + . − . . ± .
13+ magn. cycle 3 . ± .
21 2 . ± .
23 14 . + . − . . + . − . . + . − . . + . − . + linear BIS 3 . ± .
25 2 . ± .
24 14 . + . − . . ± .
059 0 . + . − . . + . − . + BIS 𝑃 rot st harm. 3 . ± .
25 1 . ± .
25 14 . ± . . ± .
012 0 . ± .
017 0 . ± . 𝑃 rot nd harm. 3 . ± .
25 2 . + . − . . + . − . . + . − . . + . − . . + . − . + BIS 𝑃 rot rd harm. 3 . ± .
24 2 . ± .
25 14 . + . − . . + . − . . + . − . . + . − . + magn. cycle 3 . ± .
21 2 . ± .
20 14 . + . − . . ± .
026 0 . + . − . . + . − . + linear log 𝑅 (cid:48) HK . ± .
21 2 . ± .
23 14 . + . − . . + . − . . + . − . . + . − . + 𝐹 𝐹 (cid:48) . ± .
26 1 . + . − . . ± . . + . − . . + . − . . + . − . Gaussian process 3 . ± .
22 2 . + . − . . ± . . + . − . . + . − . . + . − . Circular low activity 3 . + . − . . ± .
22 14 . + . − . . + . − . - -Kepler low activity 3 . + . − . . + . − . . ± .
016 53 . + . − . . + . − . . + . − . + BIS 𝑃 rot rd harm. 3 . ± .
25 2 . ± .
24 14 . ± . . + . − . . + . − . . + . − . + 𝐹 𝐹 (cid:48) . ± .
29 1 . + . − . . + . − . . + . − . . + . − . . + . − . Gaussian process low activity 3 . + . − . . + . − . . ± . . + . − . . + . − . . + . − . Circular high activity 2 . ± .
40 1 . + . − . . ± .
30 19 . + . − . - -Kepler high activity 2 . ± .
40 1 . + . − . . + . − . . + . − . . + . − . . + . − . + BIS 𝑃 rot rd harm. 2 . + . − . . + . − . . − . − . ±
10 0 . + . − . . + . − . + 𝐹 𝐹 (cid:48) . ± .
33 1 . + . − . . + . − . + − . + . − . . + . − . Gaussian process high activity 2 . ± .
34 1 . + . − . . + . − . . + . − . . ± .
036 0 . + . − . inferred by our various parametric models (Table 11), and agreeswell with the ∼
40 d period estimated by Lovis et al. (2011).The combined RV semi-amplitude ascribed to activity underthe GP model (taking into account both the convective and rota-tional terms in the model, i.e. 𝑉 𝑐 and 𝑉 𝑟 ) is a little under 1 m s − ;meanwhile, the fitted RV jitter term under the GP model is of or-der 2 m s − . By contrast, the ‘Kepler + BIS 𝑃 rot rd harm.’ modelascribes only around 30 cm s − of RV variability to activity, whileabsorbing nearly 6 m s − of variability as white-noise jitter (seeFig. 5). These differences are not surprising, given that our quasi-periodic GP function draws do not need to have constant amplitudesand phases, need not be strictly periodic, etc.: consequently, the GPcan fit more activity-related variability than the more inflexible para-metric model, which must instead absorb the stellar red noise via alarge jitter term.The GP model also allows us to constrain the evolution time-scale of active regions: the inferred value of 𝜆 e = ±
10 d (Table 12)suggests an active region lifetime (or spot evolution time scale) ofabout two rotation periods, which is significantly longer than isobserved for active regions on the Sun (Bradshaw & Hartigan 2014).The hyper-parameter 𝜆 p is harder to interpret physically, though the value 𝜆 p ∼ 𝑃 rot rd harm. + magn. cycle’ estimatethese variations to have a period of 3682 ±
14 d and 3681 ±
19 dwith RV semi-amplitudes of 1 . + . − . m s − and 1 . + . − . m s − ,respectively. The estimate for the RV semi-amplitude of the mag-netic cycle in the first model may be affected by the lack of theshort-period activity, though the MAP values for both period andsemi-amplitude are consistent within their 1 𝜎 credible intervals. In summary, our favoured activity model is the GP model; we mightconsider the ‘Kepler + BIS 𝑃 rot rd harm.’ (or arguably the 𝐹𝐹 (cid:48) )model to be the best of the parametric activity models, even though MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 Table 11.
HD 13808: Rotation periods determined by the different models.Model Rotation period (d)Kepler+ BIS 𝑃 rot st harm. 36 . + . − . + BIS 𝑃 rot nd harm. 35 . + . + . + BIS 𝑃 rot rd harm. 32 . + . − . + magn. cycle 36 . + . − . Gaussian process 38 . ± . 𝑃 rot rd harm. 36 . + . − . Gaussian process low activity 41 . + . − . Kepler high activity+ BIS 𝑃 rot rd harm. 37 . ± . . ± . Table 12.
MAP values and ± 𝜎 credible intervals for the parameters of theactivity-related parameters of the GP model.Parameter Inferred value Units 𝑉 𝑐 . ± .
14 m s − 𝑉 𝑟 . + . − . m s − 𝐿 𝑐 . + . − . - 𝐵 𝑐 . + . − . m s − 𝐵 𝑟 . + . − . m s − 𝑃 GP . ± .
47 d 𝜆 p . + . − . - 𝜆 e ±
10 d it falls short of the GP. Both the GP and the ‘Kepler + BIS 𝑃 rot rd harm.’ models favour a two-planet solution over a one-planetsolution; they both strongly reject solutions containing more thantwo Keplerians, supporting our non-planetary interpretation for any‘extra’ Keplerians. In addition, they both show amongst the low-est residual RV scatters (Table 8). However, the GP model is theonly one to reject under all circumstances solutions containing morethan two Keplerians, and to reliably detect the second planetary sig-nal even when considering the high-activity subset of observations.Despite its other merits, the ‘Kepler high activity + BIS 𝑃 rot rd harm.’ model has the unusual disadvantage of having locked on towhat appears to be a spurious or at least unusually short stellar rota-tion period, unlike the GP model and some of the other parametricactivity models (see Table 11).The MAP values and their corresponding planetary charac-teristics for our two favoured models are summarised in Table 13.The inferred minimum masses for HD 13808 b and HD 13808 care ∼ 𝑀 ⊕ and ∼ 𝑀 ⊕ (model dependent), and their orbitalsemi-major axes around the host star are ∼ .
11 AU and 0 .
26 AU,respectively. As Neptune has a mass of ∼ 𝑀 ⊕ , this leads us tothe conclusion that both planets are likely warm Neptunes.The phase folded RVs for both planets and both models aredisplayed in Fig. 4, and the corner plot corresponding to the poste-riors of both models is presented in Fig. 5. The comparison of the posteriors of both favoured models shows that – despite taking verydifferent approaches to modelling stellar activity – the planetaryparameters do agree within their 1 𝜎 credible intervals. It is alsointeresting to note how narrow (and likely wildly over-optimistic)the posterior distribution for the stellar rotation period inferred bythe harmonic model is.To demonstrate the goodness-of-fit, the RV predictions of ourtwo favoured models in comparison to the data are shown in twoparts of the data set in Fig. 6, one being in the high- and one in thelow-activity subsets of observations. We presented a comprehensive study of the HD 13808 RV dataprovided by the HARPS spectrograph. In particular, we tested andcompared multiple approaches to stellar activity modelling.We found that Bayesian evidence comparison is, by itself, notsufficient for deciding on the number of planets present: when stellaractivity is not modelled adequately, extra ‘planets’ might be favouredto account for residual variability. Only our GP model favoured twoplanets (both of which we believe to be extremely secure detec-tions) over a higher number of planets under all circumstances weconsidered, including modelling only a high-activity subset of ob-servations; the parameters for the additional ‘planets’ suggestedby other models were highly variable and model-sensitive, and/orcorresponded to likely activity variability (cf. periods in Table 9,particularly for the high-activity subset). However, some of our morecomplex stellar activity models such as the multi-harmonic and 𝐹𝐹 (cid:48) stellar activity models did seem to describe the activity variabilitybetter than more simplistic parametric models, based both on the de-crease in residual RV rms scatter and on the decrease in evidence fora three-planet solution vs. the two-planet solution. In short, modelcomparison using Bayesian evidence is only as good as the modelsconsidered.Reassuringly, when modelling the full set of observations, theRV signals of HD 13808 b and HD 13808 c were without exceptiondetected at a > 𝜎 level with every activity model we considered,with the planets’ characteristics remarkably consistent across allmodels. Taking all the observations into account, our two favouredmodels found the two planets HD 13808 b and HD 13808 c to orbittheir host star with periods of ∼ . . ∼ .
11 and 0 .
26 AU, with minimum masses of 11 and 10 𝑀 ⊕ , re-spectively. The planetary parameters inferred from the low-activityobservations alone are also mostly within agreement with thosederived from the full data set.The influence of stellar activity on the planets’ characteristicswas evident, however, when considering the parameters yieldedby the models applied to the high-activity data subset where thelog 𝑅 (cid:48) HK value of the star was > − .
90. For example, even under asingle model (‘Kepler high activity + BIS 𝑃 rot rd harm.’), the MAPRV semi-amplitude of HD 13808 b decreases from 3 .
71 m s − to2 .
19 m s − when moving from the low- to high-activity subset ofobservations, although these two semi-amplitudes are inconsistentat only a ∼ 𝜎 level. Using the GP – the only model reliably todetect HD 13808 c in the high-activity subset of observations – theMAP RV semi-amplitude of HD 13808 c decreases from 1 .
78 m s − to 1 .
04 m s − , although in this case at least the semi-amplitudes doremain consistent within ∼ 𝜎 .To conclude, even though in our case the planetary parametervalues were relatively insensitive to the choice of stellar activitymodel when analysing the full data set (246 measurements), cau- MNRAS000
04 m s − , although in this case at least the semi-amplitudes doremain consistent within ∼ 𝜎 .To conclude, even though in our case the planetary parametervalues were relatively insensitive to the choice of stellar activitymodel when analysing the full data set (246 measurements), cau- MNRAS000 , 1–17 (2020) E. Ahrer et al.
Table 13.
MAP values and ± 𝜎 credible intervals of the planetary parameters of the two planets in the HD 13808 system for the two favoured models:‘Gaussian process’ and ‘Kepler + BIS 𝑃 rot rd harm.’ HD 13808 b HD 13808 cParameter BIS 𝑃 rot rd harm. Gaussian process BIS 𝑃 rot rd harm. Gaussian process 𝑃 (d) 14 . + . − . . ± . . + . − . . + . − . 𝐾 (m s − ) 3 . ± .
24 3 . ± .
22 2 . ± .
25 2 . + . − . 𝑒 . + . − . . + . − . . + . − . . + . − . 𝑀 sin ( 𝑖 ) ( 𝑀 ⊕ ) 11 . ± .
80 11 . + . − . . ± . . + . − . 𝑎 (AU) 0 . ± . . + . − . . ± . . + . − . (a) Kepler + BIS 𝑃 rot rd harm.: 14.2 d (b) Kepler + BIS 𝑃 rot rd harm.: 53.8 d(c) Gaussian process: 14.2 d (d) Gaussian process: 53.8 d Figure 4.
Phase folded RVs of HD 13808 for the planets HD 13808 b (left) and HD 13808 c (right) i.e. the leftover RV when the other planetary signal, thestellar activity model and the polynomial offset is subtracted. The black dots are the averaged RV values computed over bins of 36 degrees in mean longitudeand the error bars correspond to the rms of the mean, while the gray points are the observed RV values with their error. The red line represents the calculatedcurve of the planet with the MAP planetary parameter values of each model. tion would be essential on smaller data sets, or indeed when tryingto characterise weaker planetary signals. In particular, the changesin the inferred planet parameters when considering only the high-activity subset of our data underscores the importance of stellaractivity modelling even for a host star nominally classified as ‘inac-tive’.
ACKNOWLEDGEMENTS
The authors are grateful to the anonymous reviewer for con-structive feedback that helped improve this manuscript, and toXavier Dumusque for a valuable discussion. VMR acknowledges the Royal Astronomical Society and Emmanuel College for finan-cial support, and the Cambridge Service for Data Driven Discovery(CSD3) – operated by the University of Cambridge Research Com-puting Service ( ), provided by Dell EMCand Intel using Tier-2 funding from the Engineering and Physi-cal Sciences Research Council (capital grant EP/P020259/1), andDiRAC funding from the Science and Technology Facilities Coun-cil ( ) – for providing computing resources. A.M.acknowledges support from the senior Kavli Institute Fellowships.
MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 . . . K b (m s − ) σ + R V ( m s − ) K a c t ( m s − ) P r o t ( d ) . . . . e c . . . e b . . . P c ( d ) . . . P b ( d ) . . . . K c ( m s − ) . . . K c (m s − ) .
176 14 . P b (d) . . . P c (d) . . . e b . . . . e c
32 34 36 38 40 P rot (d) K act (m s − ) σ +RV (m s − ) Kepler + BIS P rot rd harm.Gaussian process Figure 5.
Corner plot of the Keplerian semi-amplitudes 𝐾 𝑏 , 𝐾 𝑐 , periods 𝑃 𝑏 , 𝑃 𝑐 and eccentricities 𝑒 𝑏 , 𝑒 𝑐 of HD 13808 b and HD 13808 c, respectively,under our two favoured models: the GP model (blue) and the Kepler + BIS 𝑃 rot rd harmonic model (red). 𝐾 act is the semi-amplitude of activity-dependent RVterms combined in quadrature, for the harmonic model or the GP model, and 𝜎 + RV is the additive noise (‘jitter’) parameter in either model. The dark and lightfilled regions correspond, respectively, to 1 𝜎 (39 . 𝜎 (86 . DATA AVAILABILITY
The data underlying this article will be available via VizieR at CDS.
REFERENCES
Aigrain S., Pont F., Zucker S., 2012, MNRAS, 419, 3147Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R.,2018, AJ, 156, 58 Balan S. T., Lahav O., 2009, MNRAS, 395, 1936Baranne A., et al., 1996, A&AS, 119, 373Boisse I., et al., 2009, A&A, 495, 959Boisse I., Bonfils, X. Santos, N. C. 2012, A&A, 545, A109Bonfils X., et al., 2018, A&A, 613, A25Bradshaw S. J., Hartigan P., 2014, ApJ, 795, 79Chaplin W. J., Cegla H. M., Watson C. A., Davies G. R., Ball W. H., 2019,AJ, 157, 163Collier Cameron A., et al., 2019, MNRAS, 487, 1082Delgado Mena E., et al., 2019, A&A, 624, A78MNRAS000
Aigrain S., Pont F., Zucker S., 2012, MNRAS, 419, 3147Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R.,2018, AJ, 156, 58 Balan S. T., Lahav O., 2009, MNRAS, 395, 1936Baranne A., et al., 1996, A&AS, 119, 373Boisse I., et al., 2009, A&A, 495, 959Boisse I., Bonfils, X. Santos, N. C. 2012, A&A, 545, A109Bonfils X., et al., 2018, A&A, 613, A25Bradshaw S. J., Hartigan P., 2014, ApJ, 795, 79Chaplin W. J., Cegla H. M., Watson C. A., Davies G. R., Ball W. H., 2019,AJ, 157, 163Collier Cameron A., et al., 2019, MNRAS, 487, 1082Delgado Mena E., et al., 2019, A&A, 624, A78MNRAS000 , 1–17 (2020) E. Ahrer et al. (a) Kepler + BIS 𝑃 rot rd prediction(b) Kepler + BIS 𝑃 rot rd residuals(c) Gaussian process model prediction(d) Gaussian process model residuals Figure 6.
Model prediction (red) with posterior predictive uncertainty (shaded blue) and their residuals for ‘Kepler + BIS 𝑃 rot rd harm.’ and ‘Gaussian process’models; RV measurements and corresponding uncertainties appear in black. Two representative subsets of observations are shown, one from the high-activity(left) and low-activity (right) phases of the star. MNRAS , 1–17 (2020) wo Neptune mass planets orbiting HD 13808 Dumusque X., Udry S., Lovis C., Santos N. C., Monteiro M. J. P. F. G.,2011a, A&A, 525, A140Dumusque X., et al., 2011b, A&A, 535, A55Dumusque X., et al., 2012, Nature, 491, 207Dumusque X., et al., 2014, The Astrophysical Journal, 789, 154Faria J., et al., 2016, A&A, 589, A35Faria J. P., et al., 2020, A&A, 635, A13Feroz F., Hobson M. P., 2013, MNRAS, 437, 3540Feroz F., Hobson M., Bridges M., 2009, MNRAS, 398, 1601Feroz F., Balan S. T., Hobson M. P., 2011, MNRAS, 415, 3462Fischer H., 2011, A History of the Central Limit Theorem. Sources andStudies in the History of Mathematics and Physical Sciences, Springer,doi:10.1007/978-0-387-87857-7Foreman-Mackey D., Agol E., Angus R., Ambikasaran S., 2017, The Astro-nomical Journal, 154, 220Fressin F., et al., 2011, The Astrophysical Journal Supplement Series, 197,5Gelman A., Carlin J. B., Stern H. S., Dunson D. B., Vehtari A., Rubin D. B.,2013, Bayesian data analysis. Chapman and Hall/CRCGillon M., et al., 2017, A&A, 601, A117Gladman B., 1993, Icarus, 106, 247Goodman S. N., 1999, Annals of internal medicine, 130, 1005Gray R. O., Corbally C. J., Garrison R. F., McFadden M. T., Bubar E. J., Mc-Gahee C. E., O’Donoghue A. A., Knox E. R., 2006, The AstronomicalJournal, 132, 161Grunblatt S. K., Howard A. W., Haywood R. D., 2015, ApJ, 808, 127Hall R., Thompson S., Handley W., Queloz D., 2018, MNRAS, 479, 2968Handley W., Hobson M., Lasenby A., 2015, MNRAS, 453, 4384Hatzes A., et al., 2010, A&A, 520, A93Haywood R., et al., 2014a, International Journal of Astrobiology, 13, 155Haywood R. D., et al., 2014b, MNRAS, 443, 2517Henry T. J., Soderblom D. R., Donahue R. A., Baliunas S. L., 1996, TheAstronomical Journal, 111, 439Higson E., Handley W., Hobson M., Lasenby A., 2018, Bayesian Anal., 13,873Jeffreys H., 1983, Theory of Probability. International series of mono-graphs on physics, Clarendon Press, https://books.google.co.uk/books?id=EbodAQAAMAAJ
Jurgenson C., Fischer D., McCracken T., Sawyer D., Szymkowiak A., DavisA., Muller G., Santoro F., 2016, Proc. SPIE, 9908, 99086Kharchenko N. V., 2001, Kinematika i Fizika Nebesnykh Tel, 17, 409Lomb N., 1976, Ap&SS, 39, 447Lovis C., Pepe F., 2007, A&A, 468, 1115Lovis C., Mayor M., Pepe F., Queloz D., Udry S., 2008, in Fischer D., RasioF. A., Thorsett S. E., Wolszczan A., eds, Astronomical Society of thePacific Conference Series Vol. 398, Extreme Solar Systems. p. 455Lovis C., et al., 2011, arXiv e-prints, p. arXiv:1107.5325Malavolta L., et al., 2017, The Astronomical Journal, 153, 224Mamajek E. E., Hillenbrand L. A., 2008, The Astrophysical Journal, 687,1264Mayor M., Queloz D., 1995, Nature, 378, 355Mayor M., et al., 2003, The Messenger, 114, 20Mayor M., et al., 2009a, A&A, 493, 639Mayor M., et al., 2009b, A&A, 507, 487Mayor M., et al., 2011, arXiv:1109.2497Mortier A., Faria J. P., Correia C. M., Santerne A., Santos N. C., 2015,A&A, 573, A101Nelson B. E., et al., 2020, AJ, 159, 73Pasquini L., et al., 2008, Proc. SPIE, 7014, 70141IPepe F., Mayor M., Galland F., Naef D., Queloz D., Santos N. C., Udry S.,Burnet M., 2002, A&A, 388, 632Pepe F., et al., 2014, Astronomische Nachrichten, 335, 8Queloz D., et al., 2001, A&A, 379, 279Queloz D., et al., 2009, A&A, 506, 303Rajpaul V. M., 2017, PhD thesis, University of OxfordRajpaul V., Aigrain S., Osborne M., Reece S., Roberts S., 2015, MNRAS,452, 2269Rajpaul V., Aigrain S., Roberts S., 2016, MNRAS, 456, L6 Rajpaul V., Buchhave L. A., Aigrain S., 2017, Monthly Notices of the RoyalAstronomical Society: Letters, 471, L125Rasmussen C. E., Williams C. K. I., 2006, Gaussian Processes for MachineLearning. MIT PressRoberts S., Osborne M., Ebden M., Reece S., Gibson N., Aigrain S., 2013,Philosophical Transactions of the Royal Society A: Mathematical, Phys-ical and Engineering Sciences, 371, 20110550Santos N. C., et al., 2002, A&A, 392, 215Santos N. C., et al., 2013, A&A, 556, A150Santos N. C., et al., 2014, A&A, 566, A35Scargle J., 1982, ApJ, 236, 835Schwab C., et al., 2016, in Proc. SPIE. p. 99087H, doi:10.1117/12.2234411Skilling J., 2006, Bayesian Anal., 1, 833Skrutskie M. F., et al., 2006, The Astronomical Journal, 131, 1163Suárez Mascareño A., et al., 2020, A&A, 639, A77Thompson S. J., et al., 2016, Proc. SPIE, 9908, 99086Weiss L. M., et al., 2016, The Astrophysical Journal, 819, 83Zechmeister M., 2018, A&A, 619, A218Zechmeister M., Kürster M., 2009, A&A, 496, 577This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000