An Apodized Kepler Periodogram for Separating Planetary and Stellar Activity Signals

8 MB·

by  Philip C. Gregory

MNRAS000,1–??(2015) Preprint1February2016 CompiledusingMNRASLATEXstylefilev3.0
An Apodized Kepler Periodogram for Separating
Planetary and Stellar Activity Signals
Philip C. Gregory1(cid:63)
1Physics and Astronomy Department, University of British Columbia, 6224 Agricultural Rd., Vancouver, BC V6T 1Z1, Canada
MNRASacceptance13Jan.2016
6
1
ABSTRACT
0
A new apodized Keplerian (AK) model is proposed for the analysis of precision ra-
2
dial velocity (RV) data to model both planetary and stellar activity (SA) induced
n
RV signals. A symmetrical Gaussian apodization function with unknown width and
a
center can distinguish planetary signals from SA signals on the basis of the span of
J
theapodizationwindow.ThegeneralmodelformapodizedKepleriansignalsincludes
9
a linear regression term between RV and the stellar activity diagnostic log(R(cid:48)hk), as
2
well as an extra Gaussian noise term with unknown standard deviation. The model

parametersareexploredusingaBayesianfusionMCMCcode.Adifferentialversionof
]
R theGeneralizedLomb-Scargleperiodogramthatemploysacontroldiagnosticprovides
S an additional way of distinguishing SA signals and helps guide the choice of new peri-
. ods.ResultsarereportedforarecentinternationalRVblindchallengewhichincluded
h
multiple state of the art simulated data sets supported by a variety of stellar activity
p
diagnostics. In the current implementation, the AK method achieved a reduction in
-
o SA noise by a factor of approximately 6. Final parameter estimates for the planetary
r candidates are derived from fits that include AK signals to model the SA components
t
s andsimpleKeplerianstomodeltheplanetarycandidates.Preliminaryresultsarealso
a reported for AK models augmented by a moving average component that allows for
[
correlations in the residuals.

1
Key words: stars: planetary systems; methods: statistical; methods: data analysis;
v
techniques: radial velocities.
5
0
1
8
0
1 INTRODUCTION and EXPRES that aim to improve RV precision by a fac-
.
1 torofapproximately10overthecurrentbestspectrographs,
At the “Towards Other Earths II” meeting held in Porto
0 HARPS and HARPS-N. Clearly, the success of these devel-
PortugalinSeptember2014,onethemeemergedthatbares
6 opments hinges on our ability to distinguish true planetary
1 directly on the subject of this paper. Intrinsic stellar ac- signals from SA induced signals. At the same time good
: tivity (SA) has become the main limiting factor for planet
v progresshasbeenmadeinsimulatingstellaractivitysignals.
searches in both transit and radial velocity (RV) data. SA
i AtthePortomeetingXavierDumusquechallengedthecom-
X includes stellar oscillations, granulation, spots and facu-
munitytoalargescaleblindtestusingsimulatedRVdataat
lae/plagesandlong-termmagneticactivitycycles(see,e.g.,
r the0.7m/slevelofprecision,tounderstandthelimitations
a Saar & Donahue (1997), Santos et al. (2010), Boisse et al.
ofpresentsolutionstodealwithstellarsignalsandtoselect
(2011); Dumusque et al. (2011a)). Observational strategy
thebestapproach.TheRVchallengedata1 weregenerated
(Dumusqueetal.2011)canmitigateoraverageoutsomeof
using a modification of the SOAP code, which is described
theseactivity-inducedsignal.Thepresenceofactiveregions
in Boisse et al. (2012) and Dumusque et al. (2014).
onthestellarsurfacecanshowperiodicitiesandamplitudes
ThispaperdescribesanewapproachtoRVanalysisus-
similar to those induced by real planetary signals. Indeed,
ingapodizedKeplerianmodelstodistinguishbetweenplan-
these signals can mimic planetary signals (e.g., Figueira et
etary and SA induced RV signals. The methodology is de-
al.(2010),Santosetal.(2014),Haywood(2014),Robertson
veloped in the next section using the initial Challenge test
et al. (2015)).
CurrentlythebestRVprecisionisatthem/slevelbut
new spectrographs are under development like ESPRESSO
1 Details of the RV challenge data sets and results (Du-
musque in preparation 2016) are available from: http://rv-
(cid:63) E-mail:[email protected] challenge.wikispaces.com/
(cid:13)c 2015TheAuthors
2 P. C. Gregory
data set for which the signals present were known in ad- radial velocity is given by
vance. A summary of the method is given in Section 2.7.
This is followed by our analysis and results for the first five (cid:88)m (cid:20) (t −t )2
RV Challenge data sets. y(ti) = V+ Kjexp[− i2τ2ja ]×
j=1 j
(cid:21)
(cos{θ(t +χP )+ω }+e cosω )
i j j j j
+ β rhk(t), (1)
i
2 METHODOLOGY
where rhk ≡ log(R(cid:48)hk). The model also includes an extra
2.1 Apodized Keplerian models Gaussianwhitenoisetermwithunknownstandarddeviation
swhichcanaccountforstellarjitter.Theequationinvolves
ThefirststepinthenewapproachemploysapodizedKeple- 7×m+3 unknown parameters
rian(AK)modelsinwhichthesemi-amplitudeoftheKepler
RV model is multiplied by a symmetrical Gaussian of un- V = a constant velocity.
knownwidthτandwithanunknowncenteroftheapodizing β= the log(R(cid:48)hk) correlation parameter.
windowta.Thesamemodelisusedtofitboththeplanetary Kj= velocity semi-amplitude for jth signal.
and SA induced RV signals. Since a true planetary signal τ = apodization time constant for jth signal.
j
spans the duration of the data the apodization time, τ, will t = center of apodizing window for jth signal.
ja
be large while SA induced signals generally vary on shorter P = the period for the jth signal.
j
timescales.Oneexceptiontothisistherapidlyrotating(P e = the eccentricity for the jth signal.
j
= 0.59 d) active M4 dwarf GJ 1243 which exhibits a star ω = the longitude of periastron for jth signal.
j
spot signal that remains stable over the entire baseline of χ =thefractionoftheperiod,priortothestartofdata
j
theKeplermission(Davenportetal.2015).Asanadditional
taking, that periastron occurred at. For a planetary signal,
checkonthevalidityofplanetarycandidatesignals,weem- χ P =thenumberofdayspriortot =0thatthestarwasat
j j i
ployed a control diagnostic periodogram which is explained periastron, for an orbital period of P days. The time t =0
j i
in Section 2.4.
corresponds to the unweighted average of the observation
InadditiontotheRVmeasurements,eachofthe15chal-
times.
lenge data sets includes simultaneous observation of three s= standard deviation of extra Gaussian noise term.
stellar activity diagnostics. Two of these come from addi-
tional information on the spectral line shape that are ex- The item θ(t +χ P ) = the true anomaly, the angle of the
i j j
tractedfromthecrosscorrelationfunction(CCF),theaver-
star in its orbit relative to periastron at time t.
i
age shape of all spectral lines of the star. These two shape
As described in more detail in Gregory 2011, we em-
parameters are the CCF width (FWHM) and bisector span
ployed a re-parameterization of χ and ω to improve the
(BIS).Thethirddiagnosticlog(R(cid:48)hk)isbasedtheCaIIH&
MarkovchainMonteCarlo(MCMC)convergencespeedmo-
K line flux that is sensitive to active regions on the stellar
tivatedbytheworkofFord(2006).Thetwonewparameters
surface. A preliminary analysis of the first 5 data sets plus are ψ = 2πχ+ω and φ = 2πχ−ω. Parameter ψ is well de-
testdataindicatedastrongcorrelationbetweenRVandthe
termined for all eccentricities. Although φ is not well deter-
log(R(cid:48)hk)diagnosticandslightlyreducedcorrelationwiththe
minedforloweccentricities,itisatleastorthogonaltotheψ
FWHMdiagnostic.Forthedatasetsanalyzedinthisstudy,
parameter.Weuseauniformpriorforψintheinterval0to
thecorrelationwiththeBISdiagnosticwasgreatlyreduced 4π and uniform prior for φ in the interval −2π to +2π. This
and not considered further.
insuresthatapriorthatiswraparoundcontinuousin(χ,ω),
ThegeneralmodelformapodizedKepleriansignalsalso
maps into a wraparound continuous distribution in (ψ,φ).
included a linear regression term 2 between RV and the SA
ToaccountfortheJacobianofthisre-parameterizationitis
diagnosticlog(R(cid:48)hk)aswellasanextraGaussiannoiseterm
necessary to multiply the Bayesian integrals by a factor of
withunknownstandarddeviation s.Thecorrelationtermis (4π)−np, where np= the number of periods in the model.
expected to be particularly useful in removing SA signals
The AK models were explored using an automated fu-
associated with the star’s long-term magnetic cycle (Du-
sion MCMC algorithm (FMCMC), a general purpose tool
musque2011c). It is also expected to partially remove SA
fornonlinearmodelfittingandregressionanalysis(Gregory
associated with the rotation of active regions (spots and
2011, 2013). FMCMC is a special version of the Metropo-
plages).
lis MCMC algorithm that incorporates parallel tempering,
For an m signal apodized Kepler model the predicted
genetic crossover operations, and an automatic simulated
annealing. Each of these features facilitate the detection of
a global minimum in chi-squared in a highly multi-modal
environment. By combining all three, the algorithm greatly
increasestheprobabilityofrealizingthisgoal.TheFMCMC
2 Our treatment of this regression is only approximate because
is controlled by a unique adaptive control system that au-
weonlyhavemeasuredvaluesofthelog(R(cid:48)hk)diagnosticnotthe
tomatesthetuningoftheMCMCproposaldistributionsfor
truevalues.Ingeneral,theeffectofthemeasurementerrorsisto
efficientexplorationofthemodelparameterspaceevenwhen
reduce anycorrelation. Abetter approach wouldbe to treat the
true values of the diagnostic as hidden parameters that can be the parameters are highly correlated. The AK models com-
marginalized in a Bayesian hierarchical analysis. The virtues of binedwiththeFMCMCalgorithmconstituteamulti-signal
thecurrentapproximatetreatmentaresimplicityandspeed. apodized Keplerian periodogram.
MNRAS000,1–??(2015)
Apodized Kepler Periodogram 3
2.2 Priors Lomb-Scargle (GLS) periodograms (Zechmeister & Ku¨rster
2009).Tofirstoderthethreesignalsappearhighlycorrelated
AmoredetaileddescriptionofFMCMCisavailableinChap-
butcloseinspectionoftheperiodogramsindicatedifferences
ter 1 of the ‘Supplement to Bayesian Logical DataAnalysis
between the RV data and the diagnostics in the vicinity of
forthePhysicalSciences,’afreesupplementavailableinthe
theknown16dplanetarysignal.TheGLSperiodogramp(ω)
resources section of the Cambridge University Press web-
measures the relative χ2 -reduction, p(ω) , as a function of
site for my Textbook ‘Bayesian Logical Data Analysis for
frequencyωandisnormalisedtounitybyχ2 (theχ2 forthe
thePhysicalSciences:AComparativeApproachwithMath- 0
weightedmeanofthedata).Thequantity p(ω)isnormalized
ematica Support.’ A Mathematica implementation of fusion
to unity with p(ω)=0 indicating no improvement of the fit
MCMC is also available from the resource section. Chapter
and p(ω) = 1 a“perfect”fit. For all the GLS periodogram
1alsoincludesacomparisonofthreemarginallikelihoodes-
plots the abscissa is expressed as a period where P = 2π.
timatorsusedforBayesianmodelcomparisonandconcludes ω
Frequentist p-values 3 of 0.1, 0.01, 0.001 are shown in the
infavoroftheNestedRestrictedMonteCarlo(NRMC)esti-
upper RH periodogram 4.
matorwhichisisusedinthiswork.Thesupplementincludes
In Fig. 2, the upper left panel shows the RV data af-
a detailed discussion of the priors adopted by this author
ter removing the best linear regression fit with log(R(cid:48)hk) as
for exoplanet RV analysis. Two different frequency priors
(cid:112) (cid:112) (cid:112) theindependentvariable.ThisisreferredtoasRV(rhkcor-
were discussed: (a) p(f|I)=1/ f ×0.5/( f − f ) and (b)
H L rected)inthepanelcaption.Removingthecorrelationwith
p(f|I)=1/f ×1/ln(f /f ), where f and f are the prior up-
H L H L log(R(cid:48)hk)reducedthestandarddeviationoftheRVdatafrom
per and lower frequency bounds, respectively. In this work
8.55to2.72m/s.Wealsotriedremovingalinearregression
(a) was used primarily for parameter estimation purposes
model with log(R(cid:48)hk) and the FWHM as two independent
which helps in the detection of shorter period signals while
regression variables. This resulted in only a very slight re-
thescaleinvariantform(b)wasadoptedformodelcompar-
duction in the standard deviation of the residuals to 2.71
ison purposes which utilized NRMC method for marginal
m/s. For this reason, in the apodized Kepler model fitting
likelihood estimation.
describedbelowweonlyincludeasinglelog(R(cid:48)hk)regression
Aneweccentricitypriorwasutilizedinthisworkinter-
term as indicated in Equation (1).
mediatebetweenKipping’slowperiod(P<382.3d)caseand
hisaverageforall396planets(Kipping2013).Theequation
for the new eccentricity prior is
2.4 Control periodogram
p(e|I)=7.354∗(1−e0.3)1.5. (2)
It proved useful to employ the FWHM diagnostic data as
The two additional parameters for a one signal AK a control after removing the best linear regression fit to
modelaretheapodizationtime,τ,andtheapodizationfunc- log(R(cid:48)hk). This reduced the standard deviation of the raw
tioncentertimet .Weusedascaleinvariantpriorforτand FWHM diagnostic from 26.2 to 5.6 m/s. This control data
a
auniformpriorfort .Theupperlimitonτwassetat4times appearsinthelowerleftpanelofFig.2andisreferredtoas
a
the duration of the data set which corresponds to a a max- FWHM (rhk corrected) in the caption. It provides another
imum variation in the window height of approximately 1% check on validity of planetary candidate signals identified
acrossthedatasetforawindowwithacenteratthemiddle
of the data or by 3% for one centered at either end of the
data. The lower limit was set to the data duration divided 3 Itiscommonpracticetorefertothep-value(e.g.,Zechmeister
by 40 which corresponds to 37 d. Some of the calculations & Ku¨rster (2009) as the false alarm probabilities (FAP) of the
highest peak in the periodogram. This is incorrect, the p-value
were redone using an upper limit of 8 times the duration
saysnothingabouttheprobabilityofthisparticularsignal,rather
and a lower limit of 10 d with no noticeable improvement
it’sthefractionofpossiblehypotheticalsignals,generatedunder
inourabilitytodistinguishSAandplanetarysignals.Much
thenullhypothesisthatthereisnorealsignalpresent,thatwould
more analysis would be needed to select an optimum lower have a p(ω)(cid:62) the actual value observed at the peak of the peri-
limit for τ. The upper and lower limits on the uniform ta odogram.EveryoneofthesehypotheticalsignalshasaFAP=1
priorwere−T to+T relativetothemeandatasampletime, forproducingthep-value.ForanysignaltohaveaFAP(cid:44)1,alter-
where T is the data duration. nativestothenullmustsometimesact.Wecannotcalculatethe
As we will shall see the primary role of the AK models FAPoftheactualpeakwithoutspecifyingthealternativestothe
is to distinguish planetary signal candidates from SA sig- nullhypothesisandhowoftentheyact.Weshouldsticktocalling
nals. Suppose the results indicate that k of the signals are the quantity a p-value and, following Tom Loredo’s suggestion,
planetary and m−k are SA signals. Final model parameter thinkofitasameasureofhowsurprisedwearetoobservesuch
a peak under the null hypothesis. The difficulty in interpreting
estimates and model comparisons are based on subsequent
p-values has been highlighted in many papers (e.g., Berger and
runs using a model of k Keplerians and m−k apodized Ke-
Sellke 1987; Sellke et al. 2001). The focus of these works is that
plerians.
p-values are commonly considered to imply considerably greater
evidenceagainstthenullhypothesisthanisactuallywarranted.
4 Theinterestingregionofp(ω)iswherethefrequentistp-valueis
2.3 Challenge test data set small((cid:28)1).Inthisregionp-value≈M∗(1−z))(N−3)/2(Zechmeister
& Ku¨rster 2009), where z is the corresponding maximum p(ω),
The test data set contained a single planetary signal in a
M=thenumberofindependentfrequenciesandN isthenumber
circular orbit with P = 16.0 d, K = 1.5 m/s and phase = ofdatapoints.Cumming(2004)recommendssetting M=∆f/δf,
0.0, together with a variety of stellar activity signals. Fig. 1 where ∆f is the frequency range examined ≈ fmax, and δf = the
shows a comparison between the raw RV test data and two resolution of the periodogram ≈ 1/T where T is the duration of
diagnostics,theFWHMandlog(R(cid:48)hk),andtheirGeneralized thedata.
MNRAS000,1–??(2015)
4 P. C. Gregory
RVtestdata GeneralizedLomb(cid:45)ScarglePeriodogram
30
0.6
s 20 (cid:76) 0.5
m 10 (cid:72)(cid:144) Ω 0.4 (cid:76)
V p 0.3 (cid:72)
R 0 0.2
0.1 p(cid:45)value(cid:61)0.001
(cid:45)10 p(cid:45)value(cid:61)0.01
0.0 p(cid:45)value(cid:61)0.1
0 200 400 600 800 1000 1200 1400 0.5 1 10 100 1000
Time d Period d
FWHMtest GeneralizedLomb(cid:45)ScarglePeriodogram
100
80 (cid:72) (cid:76) 0.6 (cid:72) (cid:76)
60 0.5
M
FWH 2400 pΩ 00..34 (cid:72)(cid:76)
0 0.2
(cid:45)20 0.1
(cid:45)40 0.0
0 200 400 600 800 100012001400 0.5 1 10 100 1000
Time d Period d
RHKtest GeneralizedLomb(cid:45)ScarglePeriodogram
0.10
(cid:72) (cid:76) (cid:72) (cid:76)
0.6
k (cid:76) 0.5
h 0.05
ogR' 0.00(cid:72) pΩ 00..34 (cid:72)(cid:76)
l 0.2
0.1
(cid:45)0.05 0.0
0 200 400 600 800 1000 1200 1400 0.5 1 10 100 1000
Time d Period d
Figure 1. The left column of panels are plots of the raw RV test data and two diagnostics, the FWHM and log(R(cid:48)hk). The test data
setisknowntohaveasingleplanetarysignalwithaperiodof16d.Foreachplotthemeanvaluehasbeensubtracted.Therighthand
(cid:72) (cid:76) (cid:72) (cid:76)
columnofpanelsshowthecorrespondingGeneralizedLomb-Scargle(GLS)periodograms.ForalltheGLSperiodogramplotstheabscissa
isexpressedasaperiodwhere P= 2π.
ω
withtheapodizedKeplerianmodel.Avalidplanetarysignal thenegativeoftheperiodogramofthecontrolandthelight
is not normally expected to have a significant counterpart graytraceshowsthedifference,theblacktraceplusthegray
intheGLSperiodogramofthecontroldata5.TheGLSpe- trace. The differential GLS is shown for 6 period intervals.
riodograms of the RV (rhk corrected) data and control are The strongest feature occurs at 16 d and on either side are
showninthetwoRHSpanels.TheperiodogramofRV(rhk the twoone year aliases. One day aliases canbe seeninthe
corrected) exhibits a dominant peak with a period of 16 d vicinity of P = 0.94 d and P = 1.06 d. For all three, the
andaforestofsmallerpeaks.Theperiodogramofthecontrol light gray trace and the black trace coincide closely near
exhibitsmultiplelowerlevelpeakswhichwereinterpretedas the positive peaks, consistent with a planetary signal or its
residual SA signals or instrumental systematics. alias.Signalsincommontobothblackanddarkgraytraces
(e.g., as occurs near P=6.3 d) indicates stellar activity. At
eachstageintheanalysistheGLSperiodogram(Fig.2)and
2.5 Differential periodograms thedifferentialGLSperiodogram(Fig.3)ofthecurrentRV
residualsprovidesguidanceonthechoiceofthenextperiod
It also proved useful to construct a differential form of the
to include in the AK model.
GLS periodogram of the RV residuals for selected period
regions as is shown in Fig. 3. The black trace is from the Fig.4showssampleFMCMCresultsforaoneAKsignal
upper right periodogram in Fig. 2. The dark gray trace is model. The lower left panel shows the span of the apodiza-
tion window within the overall data window for each signal
(gray trace for MAP values of τ and t , black for a rep-
a
5 Cuntz, Saar, & Musielak (2000) suggested that there may be resentative set of samples which in this case is completely
an observable interaction between a parent star and a close-in hidden below the gray). The lower and upper boundaries
giant planet, specifically, an external heating of the star’s outer of the apodization window are defined by t −τ and t +τ,
atmosphere. Shkolnik, Walker, & Bohlender (2003) detected the a a
respectively.
synchronous enhancement of CaII H and K emission with the
short-period (3.093 d , Mpsini = 0.84MJ) planetary orbit in HD ThetopleftpanelofFig.5showsthe1signalAKmodel
179949. residuals (i.e., m = 1 in Equation (1)). The top right panel
MNRAS000,1–??(2015)
Apodized Kepler Periodogram 5
TestRV rhkcorrected GeneralizedLomb(cid:45)ScarglePeriodogram
10
0.30
s 5 (cid:76) (cid:72) (cid:76) 0.25
m (cid:144) 0.20
0 (cid:72) Ω 0.15 (cid:76)
V p (cid:72)
R (cid:45)5 0.10
p(cid:45)value(cid:61)0.001
0.05 p(cid:45)value(cid:61)0.01
(cid:45)10 p(cid:45)value(cid:61)0.1
0.00
0 200 400 600 800 1000 1200 1400 0.5 1 10 100 1000
Time d Period d
TestFWHM rhkcorrected GeneralizedLomb(cid:45)ScarglePeriodogram
20 (cid:72) (cid:76) Control 0.08 (cid:72) (cid:76)
(cid:72) (cid:76)
M
10 0.06
H
W Ω (cid:76)
F 0 p 0.04 (cid:72)
(cid:45)10 0.02
0.00
0 200 400 600 800 1000 1200 1400 0.5 1 10 100 1000
Time d Period d
Figure 2.TheleftcolumnofpanelsareplotsofthetestdataRVandFWHMafterremovingthebestlinearregressionfitstolog(R(cid:48)hk)
(rhkcorrected).TherighthandcolumnofpanelsshowthecorrespondingGLSperiodograms.
(cid:72) (cid:76) (cid:72) (cid:76)
Figure 3. A differential GLS periodogram of selected period regions for the test data based on the RV and FWHM control (both rhk
corrected).
MNRAS000,1–??(2015)
6 P. C. Gregory
Figure 4. The upper left panel is a plot of the Log10[Prior × Likelihood] versus iteration for the 1 signal apodized Kepler (AK)
periodogramofthetestdata.TheupperrightshowsLog10[Prior×Likelihood]versusperiodshowingthe1perioddetected.Themiddle
left shows the values of the 1 unknown period parameter versus iteration number. The middle right shows the eccentricity parameter
versusperiodparameter.Thelowerleftshowstheapodizationwindowforeachsignal(graytraceforMAPvaluesofτandta,blackfor
a representative set of samples which is hidden below the gray). The lower right is a plot of the apodization time constant, τ, versus
apodizationwindowcentertime,ta.
shows the corresponding GLS periodogram. The remaining apodization window center time, t . Based on the span of
a
panelsshowthedifferentialperiodogramoftheresidualsfor the apodization window, all of the signals, with the excep-
selected period intervals. The P=16 day signal and its one tion of P=16 d, are consistent with SA. As a check on the
year and daily aliases have been successfully removed. The planetary interpretation of the 16 d signal, Fig. 3 indicates
featureatP=6.3disnowthehighestpeakintheresiduals. that p(ω) for the control trace is within 0.03 of zero for a
periodof16d.Atthispointintheanalysis,thepeakinthe
The analysis proceeds by repeatedly adding an addi- GLSperiodogramofthe5signalAKfitresidualsatP=0.95
tional AK signal to the model with an initial period esti- d was down to p(ω)=0.045. The p-value was >0.01 and no
matebasedonthepreviousresidualsandexaminingthethe further signals were extracted.
FMCMCresultsplustheGLSanddifferentialperiodograms
of the residuals. Note: the the linear log(R(cid:48)hk) correlation is
2.6 Apodized Keplerian plus ordinary Keplerian
not removed prior to each fit but is accomplished by the
models
β× rhk(t) term in Equ. 1 as part of each fit. This process
i
continuesuntilthep-value>0.01forthehighestpeakinthe TheobtainthebestP=16dparameterestimatesthemAK
GLS periodogram of the residuals and the FMCMC results modeltermsinEquation(1)werereplacedbyfourapodized
for the next more complicated model indicates there is no Keplerians plus one ordinary Keplerian to model the 16 d
welldefinedsolutionforanadditionalsignal.Intheinterest planetary signal, and the data refit. The starting FMCMC
of space we jump to our final 5 AK signal model results. parametervaluesweretakenfromthebestfitting5AKsig-
Fig.6 shows sample FMCMC results for a five AK sig- nal model. Fig. 7 shows a sample of resulting 16 d parame-
nalmodel.TheupperleftpanelisaplotoftheLog [Prior× tercorrelationplots,whichinaBayesiananalysisreflectour
10
Likelihood]versusiterationforthe5signalAKperiodogram current state of knowledge for the 16 d planetary signal in
of the test data. The upper right shows Log [Prior × Like- the test data set.
10
lihood] versus period showing the 5 periods detected. The The P and K values have been normalized by divid-
middle left shows the values of the 5 unknown period pa- ing by their true values to convert to ratios. Since the true
rameter versus iteration number. The middle right shows eccentricity = 0 the relevant phase for comparison is our
the eccentricity parameters versus period parameters. The re-parameterized phase term ψ = 2πχ+ω. In our analysis
lower left panel shows the span of the apodization window the phase is referenced to the mean observation time. This
within the overall data window for each signal (gray trace was converted to a phase referenced to the time of the first
forMAPvaluesofτandt ,blackforarepresentativesetof sampleandaddedtoanadditionalπ/2tocomparewiththe
a
samples which is mainly hidden below the gray). The lower XavierDumusque’sdefinitionofphasewhichreferstoasine
and upper boundaries of the jth apodization windows are function while ours refers to a cosine function.
defined by t −τ and t +τ , respectively. The lower right Itisnottoomeaningfultodraw anyconclusionsabout
ja j ja j
panel is a plot of the apodization time constant, τ, versus parameter biases on the basis of this one data set. Later
MNRAS000,1–??(2015)
Apodized Kepler Periodogram 7
Figure 5.ThepanelsshowtheonesignalAKmodelresidualsforthetestdataset(topleft),itsGLSperiodogram(topright),andthe
differentialperiodogramforselectedperiodranges.
in Section 4, we compare the K and P parameter estimates parameterestimatesforcandidateplanetswerebasedonthe
derived from the m AK model fits to the k AK plus (m−k) k AK plus (m−k) Kepler model fits, where k is the number
Kepler model fits for the six data sets analyzed (including of signals attributed to SA.
thetestdata).Nineplanetswerecorrectlyidentifiedinthese
Ideally,wewouldliketoseetheBayesianmarginalsspan
sixdatasets.FromthisitisclearthatAKmodelestimates
the true values of P,K, and phase slightly better than indi-
of the K parameter are significantly biased to larger values
cated by Fig. 7. This might be an indication that the true
compared to their known true values. As a result, our final
scatter is underestimated possibly as a result of correlated
MNRAS000,1–??(2015)
8 P. C. Gregory
Figure 6. The upper left panel is a plot of the Log10[Prior × Likelihood] versus iteration for the 5 signal AK periodogram of the test
data.TheupperrightshowsLog10[Prior×Likelihood]versusperiodshowingthe5periodsdetected.Themiddleleftshowsthevaluesof
the5unknownperiodparameterversusiterationnumber.Themiddlerightshowstheeccentricityparametersversusperiodparameters.
The lower left shows the apodization window for each signal (gray trace for MAP values of τ and ta, black for a representative set of
sampleswhichishiddenbelowthegray).Thelowerrightisaplotoftheapodizationtimeconstant,τ,versusapodizationwindowcenter
time,ta,wherethesymbolsrefertothedifferentsignalperiodsspecifiedinthepanelsabove.
Figure7.ParametercorrelationplotsfortheP=16dplanetarysignalinthetestdatasetbasedonamodelconsistingoffourapodized
KepleriansplusoneordinaryKeplerian.
MNRAS000,1–??(2015)
Apodized Kepler Periodogram 9
3. Remove the best linear regression fit from the raw
FWHM data to obtain the control FWHM (rhk corrected)
data, using the log(R(cid:48)hk) diagnostic as the independent re-
gression variable.
4. ComputetheGLSperiodogramofcontrolFWHM(rhk
corrected).
5. ComputethedifferentialGLSperiodogramforselected
period regions, showing on the same plot, the periodogram
ofRV(rhkcorrected)inblack,theperiodogramofthecon-
trol FWHM (rhk corrected) plotted negatively and shown
in dark gray, and the their sum shown in gray. Planetary
candidate signals for which the gray trace and black trace
are significantly different (i. e., the dark gray trace has a
p(ω) more negative than −0.05) indicates that the signal is
Figure8.Thetoprowshowstheautocorrelationfunctionofthe
possibly SA in origin.
(rhkcorrected)RVtestdataandtherawlog(R(cid:48)hk)diagnostic.The
6. FitthemodeldescribedbyEqu.1totherawRVdata,
leftpanelinthesecondrowshowstheautocorrelationfunctionof
in our case using the FMCMC algorithm. Use the highest
theRVresidualsafterremovingamodelconsistingof4apodized
Kepleriansignalsplus1Keplersignalandthelog(R(cid:48)hk)regression. peakintheGLSperiodogramoftherhkcorrectedRVdata
Thebottompanelshowsaplotofthenumberofsuchsamplepairs astheinitialestimatefortheperiodoftheAKmodelsignal,
asafunctionofthelag. regardlessofwhetherthedifferentialGLSperiodogramindi-
catesitisSAinorigin.WeareusingtheapodizedKeplerians
to model both types of signals.
noise in the residuals. One way of exploring the residuals is
7. ComputetheGLSanddifferentialperiodogramsofthe
to compute the autocorrelation function, ρ(j) equation (3).
fit residuals.
(cid:80)
ρ(j)= (cid:113)(cid:80) overlap[(xi−x(cid:113))(cid:80)(xi+j−x)] , (3) AK8.siRgenpaelawtitthheanlasintittiwaol piteermiosdbeystiimncaltueddinfgroamntahdedliatriognesatl
overlap(xi−x)2× overlap(xi+j−x)2 peak in the GLS periodogram of the previous fit residuals.
where x is the ith residual, j is the lag and x is the mean of 9. Repeatthepreviousstepuntilthep-valueofthehigh-
i
est peak in the residuals is greater than 0.01 and the FM-
the samples in the overlap region.
CMC fit indicates there is no well defined solution for an
ThetoprowofFig.8showstheautocorrelationfunction
ofthe(rhkcorrected)RVtestdataandtherawlog(R(cid:48)hk)di- additional signal.
10. Examine the span of the apodization window for the
agnostic.Note:thehighestpeakintheautocorrelationfunc-
tionforthelog(R(cid:48)hk)diagnosticcorrespondsroughlywiththe MAP parameter values and representative samples for each
star’s true rotation period of 25.05 d (vertical dashed line). signal extracted. The lower and upper boundaries of the jth
The rotation period was not provided during the competi- apodization windows were chosen to be tja−τj and tja+τj,
respectively.
tion.SAsignalsfrequentlyoccurattherotationperiodand
harmonics.Theleftpanelinthesecondrowshowstheauto- 11. If the apodization window spans the data duration
correlation function for the RV residuals for the 4 apodized and, p(ω) for the dark gray trace in the differential peri-
Kepler plus 1 Kepler model including the log(R(cid:48)hk) regres- odogram of the residuals is less negative than −0.05, then
sion. Clearly, after after removing a model consisting of 5 treat the signal as a candidate planetary signal. If the
signals and the log(R(cid:48)hk) regression, the autocorrelation of apodization window spans the data duration and p(ω) for
the residuals is looking much close to white noise but there the dark gray trace is more negative than −0.05, then treat
isstillevidenceforasmallpositiveautocorrelationforalag the signal as a possible candidate planetary signal. If the
ofoneday.Theremainingpanelshowsaplotofthenumber apodization window does not span the data duration clas-
of such sample pairs as a function of the lag. Because the sify it as a SA signal.
dataarenotuniformlysampled,foreachlagallsamplepairs 12. Maintain a list of the order in which the signals were
that differed in time by this lag ±0.1 d were utilized. extractedandthesignaldesignationasPforplanetarycan-
didate,P?forpossibleplanetarycandidate,andSAforstel-
lar activity signal.
2.7 Summary and Discussion of Methodology 13. Suppose there are m signals in total, k of which are
SA signals and (m−k) are planetary candidates or possible
In this section we summarize the steps in the methodology
planetary.FitamodelwithkAKsignalsand(m−k)ordinary
andillustrateitssuccessinmodelingthebothlong-termand
Kepleriansignalstogetherwiththeβ× rhk(t)termofEqu.1
short-term SA induced RV variations. The method consists i
to derive final parameter estimates for the planetary and
of the following steps:
possible planetary candidates.
1. Remove the best linear regression fit from the raw RV 14. ComputeBayesfactorsasafurthertestofthevalidity
data to obtain the RV (rhk corrected) data, using the SA ofanyplanetarycandidate.Thisisdiscussedfurtherstarting
diagnostic log(R(cid:48)hk) as the independent regression variable. in Section 3.2.
2. Compute the GLS periodogram of RV (rhk corrected)
and p-value levels corresponding to 0.1, 0.01, 0.001. The IntheupperleftpanelofFig.9,the(opencircles)rep-
highest peak in this periodogram is used as the initial es- resenttherawRVtestdataminustheMAPfitforthe16d
timate for the period of the AK model signal. planetwhichwedesignateasRV .Thedotsshowhow
1planCorr
MNRAS000,1–??(2015)
10 P. C. Gregory
Figure 9. The open circles in the upper left panel show the raw RV test data minus the MAP fit for the 16 d planet. The dots show
howmuchofthisremaindercanbeaccountedforbyacorrelationwiththestellaractivitydiagnosticlog(R(cid:48)hk).Thestarredpointsshow
thedifferencebetweentheopencirclesandthedots.Theotherpanelsshowexpandedversions,oneforeachobservingseason.
muchoftheRV canbeaccountedforbyacorrelation
1planCorr Table 1. Order of signals extracted from RV 1 together with
with the SA diagnostic log(R(cid:48)hk), based on a weighted least approximate parameter values and designation of signal type as
squares linear regression fit using log(R(cid:48)hk) as the indepen- planetary(P)orstellaractivity(SA).
dentregressionvariable.Thestarredpointsshowthediffer-
ence between the open circles and the dots. Removing the Period K ecc τ ta Signal
correlationwithlog(R(cid:48)hk)reducedthestandarddeviationof (d) (m/s) (d) (d) type
RV from8.38to2.29m/s.Theotherpanelsshowex-
1planCorr 33.35 2.5 0.32 1356 -133 P
pandedversions,oneforeachobservingseason.Clearly,the
log(R(cid:48)hk) regression fit tracks the long-term SA very well.
9.89 1.8 0.0 1274 -341 P
Closeinspectionoftheseasonalpanelsindicatestheregres-
sionfitonlyprovidesafirstordertrackingoftheshortterm 23.4 1.6 0.24 1747 -1204 P
variationsinRVarisingfromindividualactiveregions(spots
and plages) . 12.5 2.5 0.37 557 410 SA
10.6 1.0 0.0 494 -199 SA
8.38 2.1 0.39 38 -93 SA
IntheupperleftpanelofFig.10,the(opencircles)rep-
resent the raw RV data. The dots show the MAP fit of the
full model consisting of four apodized Keplerians plus one
ordinaryKeplerian(tomodelthe16dplanetarysignal),to-
getherwiththelinearregressiontermusinglog(R(cid:48)hk)asthe
3 RESULTS
independentregressionvariable.Thestarredpointsshowthe
differencebetweentheopencirclesandthedots.Employing The procedures outlined in Section 2 were applied to the
the full model reduced the standard deviation of from 8.55 first five RV challenge data sets and the results are given
to 1.35 m/s. The other panels show expanded versions, one below.Table1showstheorderinwhichthesignalsemerged
for each observing season. Clearly, the full model provides insuccessiveAKmodelstogetherwiththesignaldesignation
a much better fit to the RV data including the short term asplanetary(P)orstellaractivity(SA)inthefinalcolumn.
variations arising from individual active regions.
MNRAS000,1–??(2015)