Orbital Characterization of Multi-Object Exoplanetary Systems with Radial Velocity Observation

MISSING IMAGE

Material Information

Title:
Orbital Characterization of Multi-Object Exoplanetary Systems with Radial Velocity Observation
Physical Description:
1 online resource (154 p.)
Language:
english
Creator:
Guo, Peng-Cheng
Publisher:
University of Florida
Place of Publication:
Gainesville, Fla.
Publication Date:

Thesis/Dissertation Information

Degree:
Doctorate ( Ph.D.)
Degree Grantor:
University of Florida
Degree Disciplines:
Astronomy
Committee Chair:
Ford, Eric B
Committee Members:
Matcheva, Katia Ivandva
Ge, Jian
Dermott, Stanley F
Ghosh, Malay

Subjects

Subjects / Keywords:
bayesian -- evidence -- exoplanet -- keplerian -- mcmc -- model -- nbody -- orbit -- rv
Astronomy -- Dissertations, Academic -- UF
Genre:
Astronomy thesis, Ph.D.
Electronic Thesis or Dissertation
bibliography   ( marcgt )
theses   ( marcgt )
government publication (state, provincial, terriorial, dependent)   ( marcgt )

Notes

Abstract:
Exoplanets are the planets around stars other than the sun. The detection and characterization of exoplanetary systems nowadays is one of the fields of the greatest interest in astronomy and planetary science. As various exoplanets detection programs have been keeping delivering observations of potential exoplanets, their database is growing rapidly in terms of the size and the detection precision. It is crucial to make the best use of these data, not only for discovering the exoplanets, but also for determining the most appropriate model for describing the data set and measuring the model parameters (e.g. masses and orbital parameters) more accurately. One of the most productive detection methods is the Radial Velocity(RV) method. It measures the Doppler shift of a stellar spectrum caused by the gravity perturbation of the orbiting planet(s). Dispersed Fixed-delay Interferometry(DFDI) is a technique to economically implement the RV detection. Chapter 2 describes my involvement in two instrumentation projects based on this technique and focus on a method of increasing the short-term stability of the instruments. When interpreting the reduced RV data, people usually use a Keplerian model to describe the orbits of exoplanets. However, if there is more than one planet in the system, the Keplerian model neglects the interaction between planets. Theoretically, an n-body simulation should be used for describing the system in terms of accuracy. However, the N-body model is much more computationally expensive. In an attempt to determine when a full N-body model is important, I compare the Keplerian model with N-body model for characterizing the orbits of a set of simulated interacting two planet systems, which is described in Chapter 3. Given a RV data, Markov Chain Monte Carlo provides a convenient tool to sample from the posterior distribution of the model parameters. One challenge that remains is computing the Marginal Likelihood of the Posterior or the Bayesian Evidence, as it requires the calculation of a high-dimensional integral. A simple algorithm for calculating the Bayesian Evidence is introduced in Chapter 4. Chapter 5 presents several applications of the algorithm introduced in Chapter 4. And the summery of this dissertation research is given in the end.
General Note:
In the series University of Florida Digital Collections.
General Note:
Includes vita.
Bibliography:
Includes bibliographical references.
Source of Description:
Description based on online resource; title from PDF title page.
Source of Description:
This bibliographic record is available under the Creative Commons CC0 public domain dedication. The University of Florida Libraries, as creator of this bibliographic record, has waived all rights to it worldwide under copyright law, including all related and neighboring rights, to the extent allowed by law.
Thesis:
Thesis (Ph.D.)--University of Florida, 2012.
Local:
Adviser: Ford, Eric B.
Statement of Responsibility:
by Peng-Cheng Guo.

Record Information

Source Institution:
UFRGP
Rights Management:
Applicable rights reserved.
Classification:
lcc - LD1780 2012
System ID:
UFE0044871:00001


This item has the following downloads:


Full Text

PAGE 1

ORBITALCHARACTERIZATIONOFMULTI-OBJECTEXOPLANETARYSYSTEMSWITHRADIALVELOCITYOBSERVATIONByPENG-CHENGGUOADISSERTATIONPRESENTEDTOTHEGRADUATESCHOOLOFTHEUNIVERSITYOFFLORIDAINPARTIALFULFILLMENTOFTHEREQUIREMENTSFORTHEDEGREEOFDOCTOROFPHILOSOPHYUNIVERSITYOFFLORIDA2012

PAGE 2

c2012Peng-ChengGuo 2

PAGE 3

Apieceofhumbleunderstandingofthenature,tomyparents,Mr.Guo,Shoutao,andMrs.Liu,Suqin 3

PAGE 4

ACKNOWLEDGMENTS UponthenishofthisPhDdissertation,Ialmostreachedtheendofmyuntypicalacademiccareer.However,Ibelievethisisnottheterminationofmyresearchlife:exploringanewdomain,learningnewknowledge,solvingnewproblems,andcontributingtothedatabaseofthehuman'sunderstandingoftheworld.Thespecialthanksgotomyphysicsteachers,Mr.Zhang,inmymiddleschool,andMr.Zhao,myphysicsteacherinhighschool,fortheirintroductionandguidanceofmeintothephysicalscience.MyPhDadvisor,Dr.EricB.Ford,hasplayedthemostimportantroleinmakingmeaqualieddoctor.Withhisabundantacademicknowledgeandresearchexperience,hetaughtmemanyveryusefulmethodstoconductresearch.Iespeciallythankhimforhispatientsupportandencouragement,andinspiringmetogeneratetheresearchideasandconcentrateonthetopics.IthankDr.JianGe,mypreviousadvisor,amemberofmyacademiccommitteeaswell,forgivingmetheopportunitiestoparticipateinseveralcutting-edgeastronomicalinstrumentationprojects.Meanwhile,hisstrictworkingstylealsotrainedmetobecomemoreprofessionalandresponsive.IalsoespeciallythankDr.StanF.Dermott,Dr.MalayGoshandDr.KatiaIvandvaMatcheva,whoservedasmyacademiccommitteemembers,fortheirprecioustimeandsuggestions,inbothacademicproblemsandpersonalmotivation.Theirsupportsarehighlyappreciated.IparticularlythankMs.LinWang,whohasbeenselesslysupportingmeinmanyaspectsofmylifeduringthepreparationofthisdissertation.Meanwhile,withoutthehelpfrommyfriendsandfellowsintheDepartmentofAstronomy,UniversityofFlorida,especiallyCatherineCassidy,ourofcesecretary,IamafraidIwasnotabletoachievethegoalwiththePhDdegree.Itisanandunforgettableexperiencetospendmylifeasagraduatestudentwiththem. 4

PAGE 5

TABLEOFCONTENTS page ACKNOWLEDGMENTS .................................. 4 LISTOFTABLES ...................................... 8 LISTOFFIGURES ..................................... 9 ABSTRACT ......................................... 13 CHAPTER 1INTRODUCTION ................................... 15 1.1DetectionMethods ............................... 17 1.2RVDetection .................................. 19 1.3DataanalysisofRVobservation ....................... 21 1.3.1KeplerianModelv.s.N-bodyModel .................. 21 1.3.2UncertaintyoftheModelParametersandModelSelection ..... 21 1.4OrganizationofThisDissertation ....................... 22 2THECONTROLSYSTEMOFAMULTI-OBJECTRADIALVELOCITYTRACKER 24 2.1BackgroundandMotivation .......................... 24 2.2PrincipleofDFDI ................................ 25 2.2.1InstrumentationBasedontheDFDITechniques ........... 29 2.3TheExoplanetTracker ............................. 29 2.3.1PhaseLocker .............................. 30 2.3.2DisadvantagesofthePhaseLocker .................. 34 2.3.3OtherInstrumentationControlandUserInterfaceUpgrading .... 35 2.3.4CurrentStatusofET .......................... 36 2.4MARVELSPilotProgram ........................... 36 2.4.1TheTwoPointsPhaseLocker ..................... 36 2.4.2UpgradeofMARVELS ......................... 38 2.5Summary .................................... 42 3KEPLERIANMODELSVERSUSN-BODYMODELSONRADIALVELOCITYOBSERVATIONSOFRESONANTPLANETARYSYSTEMS ........... 43 3.1BackgroundandMotivation .......................... 43 3.2Methodology .................................. 44 3.2.0.1GenerationofResonantPlanetarySystem ........ 45 3.2.0.2GenerationofRadialVelocityDataSets .......... 48 3.2.1Multi-planetKeplerianFitting ..................... 50 3.3ResultsandAnalysis:ComparisonBetweenKeplerianandN-bodyModels 53 3.3.1GoodnessofFit ............................. 53 3.3.2ParameterCharacterization ...................... 57 5

PAGE 6

3.3.2.1Pi,o .............................. 59 3.3.2.2ei,o ............................... 60 3.3.2.3Ki.o .............................. 60 3.4SummaryandDiscussion ........................... 69 4BAYESIANMODELSELECTIONONEXOPLANETSRVDETECTION ..... 77 4.1Motivation .................................... 77 4.2Background ................................... 78 4.2.1Outline .................................. 79 4.3TheBayesianModelSelectionFramework .................. 80 4.3.1BayesianModelSelectionFrameworkinRVExoplanetDetection 82 4.3.2BackgroundofMarkovChainMonteCarlo,Metropolis-HastingsAlgorithmandGibbsSampling .................... 85 4.4CalculatingEvidenceBasedonMCMCSamples .............. 86 4.4.1TheFractionEstimator ......................... 87 4.4.2ErrorAnalysisofFractionEstimator .................. 88 4.4.2.1VarianceofFn(x) ...................... 89 4.4.2.2VarianceofGn(x) ...................... 90 4.4.2.3CombineGn(x)andFn(x) .................. 92 4.4.3FractionEstimatorwithMCMCSamples ............... 92 4.4.4Multi-dimensionalModelSelection .................. 93 4.4.5TestwithMCMCSamplesfromMulti-variantGaussian ....... 94 4.4.5.1GeneratingMCMCsamplesfrommulti-variantGaussian 94 4.4.5.2EvidencecalculationwithMCMCsamplesfrommulti-variantGaussian ........................... 95 4.4.5.3TherelationbetweenNMCMCandNMCI .......... 97 4.4.5.4Principlecomponentsanalysis ............... 102 4.5EvidenceCalculationBasedonMCMCSamplesfromRVData ...... 105 4.5.1HD88133 ................................ 105 4.5.2EvidenceCalculationwithFractionEstimator ............ 108 4.5.2.1Findingmodeofposteriordistribution ........... 108 4.5.2.2FractionEstimatorwithHD88133RVData ........ 110 4.5.3ComparisonwithOtherEstimators .................. 110 4.5.4ChoiceoftheSizeofVn ........................ 114 4.6SummaryandDiscussion ........................... 120 5APPLICATIONS ................................... 123 5.1BackgroundandMotivation .......................... 123 5.2ModelSelectionwithSimulatedData ..................... 125 5.2.1OnePlanetModel ........................... 126 5.2.2TwoPlanetModel ............................ 127 5.2.3ThreePlanetModel ........................... 132 5.2.4DiscussionandModelSelection .................... 134 5.3EvidenceCalculationwithRealObservations ................ 139 6

PAGE 7

5.3.1AndwithThreePlanetsandOneVelocityReference ....... 139 5.3.261VirwithThreePlanetsandTwoVelocityReferences ....... 140 5.3.3ArawithFourPlanetandThreeVelocityReferences ....... 140 5.3.4Discussion ................................ 141 5.4MARVELSReferenceStars .......................... 144 5.5BayesianModelSelectionforInstrumentAnalysis ............. 148 5.6Summary .................................... 150 REFERENCES ....................................... 152 BIOGRAPHICALSKETCH ................................ 154 7

PAGE 8

LISTOFTABLES Table page 3-1MERCURYIntegrationParameters ......................... 48 3-2Selectedknownplanetarysystemaround2:1MMR ............... 75 4-1SAMSIExoplanetWorkingGroupReferencePriors ................ 107 5-1OrbitalparametersofthetwoplanetsystemforBayesianselection ....... 126 5-2ThreeselectedplanetarysystemsforexaminingtheFractionEstimator .... 139 8

PAGE 9

LISTOFFIGURES Figure page 2-1Dispersedinterferometerschematic ........................ 26 2-2Simulatedinterferometerpatterns ......................... 28 2-3TheuserinterfaceofthePhaseLocker ...................... 31 2-4BeamsplitteropticallayoutofMARVELSpilotdesign ............... 37 2-5Thetwopointphaselocker ............................. 39 2-6ThecontroldiagramofMARVELSpilotprogram ................. 40 2-7StabilityofMARVELSpilotinstrument ....................... 41 3-1Samplesofevolvingsystemswithdifferentmassesandmassratio ....... 47 3-2Samplesofevolvingsystemsindifferentevolvingstages ............. 49 3-3ThedistributionofPi,Po,eiandeo ......................... 51 3-4AnexamplesystemwithKepleriants ....................... 54 3-5DistributionofRSS(ResidualSumofSquare)ofN-bodytsandKepleriants 56 3-6ContoursofthepvalueofF-testbetweenthegoodness-of-toftheKeplerianmodelandtheN-bodymodel ............................ 58 3-7RelativemeasurementbiasofPi .......................... 61 3-8RelativemeasurementerrorofPi .......................... 62 3-9RelativemeasurementbiasofPo .......................... 63 3-10RelativemeasurementerrorofPo ......................... 64 3-11Relativemeasurementbiasofei .......................... 65 3-12Relativemeasurementerrorofei .......................... 66 3-13Relativemeasurementbiasofeo .......................... 67 3-14Relativemeasurementerrorofeo .......................... 68 3-15RelativemeasurementbiasofKi .......................... 69 3-16RelativemeasurementerrorofKi .......................... 70 3-17RelativemeasurementbiasofKo .......................... 71 3-18RelativemeasurementerrorofKo ......................... 72 9

PAGE 10

4-1ExamplesofmarginalposteriorprobabilitydistributionsofrandomlyselecteddimensionsfromasequenceofMCMCsamples ................. 96 4-2EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma7-variantNormaldistribution. ........................ 98 4-3EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma14-variantNormaldistribution ........................ 99 4-4EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma21-variantNormaldistribution. ....................... 100 4-5EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma28-variantNormaldistribution. ....................... 101 4-6Choiceofthesub-domaininoriginalparameterspace .............. 104 4-7Choiceofthesub-domaininoriginalparameterspaceshowninthetransformedparameterspaceusingPCAtechniques ...................... 105 4-8Choiceofthesub-domaininthetransformedparameterspaceusingPCAtechniques ...................................... 106 4-9Choiceofthesub-domaininthetransformedparameterspaceusingPCAtechniquesshownintheoriginalparameterspace ................ 108 4-10RVobservationsofHD88133 ............................ 109 4-11Transforming!formodeselection ......................... 111 4-12Transformed!formodeselection ......................... 112 4-13TransformingM0formodeselection ........................ 113 4-14TransformedM0formodeselection ......................... 114 4-15EvidenceestimationofoneplanetmodelforHD88133 .............. 115 4-16EvidenceestimationuncertaintycalculatedusingFractionEstimatorofoneplanetmodelofHD88133 .............................. 116 4-17Comparisonofseveralestimatorsofthemarginalposteriorprobabilityforsingleplanetmodelsandthe17RVobservationsofHD88133 ............. 117 4-18Comparisonofseveralestimatorsofthemarginalposteriorprobabilityforsingleplanetmodelsandthe17RVobservationsofHD88133. ............. 118 4-19RelativeestimationstandarddeviationcalculatedwiththeFractionEstimatorwithoutPCA ..................................... 119 10

PAGE 11

4-20RelativeestimationstandarddeviationcalculatedwiththeFractionEstimatorusingPCA ...................................... 120 5-1EvidenceofoneplanetmodelforaN-bodysimulatedtwo-planetsystem .... 128 5-2FractionofMCMCsamplesusedforevidencecalculationofoneplanetmodelforaN-bodysimulatedtwo-planetsystem ..................... 129 5-3RelativestandarddeviationoftheevidencecalculationofoneplanetmodelforaN-bodysimulatedtwo-planetsystem ..................... 130 5-4PosteriorPDFsoftwoplanetmodelparameters .................. 131 5-5PosteriorPDFsoftwoplanetmodelparameterstransformedintoPCAspace 132 5-6EvidenceofthetwoplanetmodelofaN-bodysimulatedtwoplanetsystem .. 133 5-7FractionofMCMCsamplesforthetwoplanetmodelevidenceofaN-bodysimulatedtwoplanetsystem ............................ 134 5-8RelativestandarddeviationsoftheevidenceestimationofthetwoplanetmodelforaN-bodysimulatedtwoplanetsystem ..................... 135 5-9EvidenceofthreeplanetmodelofaN-bodysimulatedtwoplanetsystem ... 136 5-10FractionofMCMCsamplesusedforevidenceestimationofthreeplanetmodelofaN-bodysimulatedtwoplanetsystem ..................... 137 5-11RelativestandarddeviationoftheevidenceestimationofthreeplanetmodelforaN-bodysimulatedtwoplanetsystem ..................... 138 5-12EvidenceofthreeplanetmodelofAnd ...................... 141 5-13FractionofMCMCsamplesincludedinthesub-domainofthreeplanetmodelforAnd ....................................... 142 5-14RelativestandarddeviationoftheevidenceestimationofthreeplanetmodelforAnd ....................................... 143 5-15Evidenceofthreeplanetmodelwithtwovelocityreferencesfor61Vir ..... 144 5-16FractionofMCMCsamplesinsub-domainfor61Vir ............... 145 5-17Relativestandarddeviationoftheevidenceestimationofthreeplanetmodelwithtwovelocityreferencesfor61Vir ....................... 146 5-18EvidenceoffourplanetmodelofafourplanetsystemAra ........... 147 5-19FractionofMCMCsamplesincludedinthesub-domainforAra ........ 148 11

PAGE 12

5-20RelativestandarddeviationoftheevidenceestimationoffourplanetmodelforAra ........................................ 149 12

PAGE 13

AbstractofDissertationPresentedtotheGraduateSchooloftheUniversityofFloridainPartialFulllmentoftheRequirementsfortheDegreeofDoctorofPhilosophyORBITALCHARACTERIZATIONOFMULTI-OBJECTEXOPLANETARYSYSTEMSWITHRADIALVELOCITYOBSERVATIONByPeng-ChengGuoDecember2012Chair:EricB.FordMajor:Astronomy Exoplanetsaretheplanetsaroundstarsotherthanthesun.Thedetectionandcharacterizationofexoplanetarysystemsnowadaysisoneoftheeldsofthegreatestinterestinastronomyandplanetaryscience.Asvariousexoplanetsdetectionprogramshavebeenkeepingdeliveringobservationsofpotentialexoplanets,theirdatabaseisgrowingrapidlyintermsofthesizeandthedetectionprecision.Itiscrucialtomakethebestuseofthesedata,notonlyfordiscoveringtheexoplanets,butalsofordeterminingthemostappropriatemodelfordescribingthedatasetandmeasuringthemodelparameters(e.g.massesandorbitalparameters)moreaccurately. OneofthemostproductivedetectionmethodsistheRadialVelocity(RV)method.ItmeasurestheDopplershiftofastellarspectrumcausedbythegravityperturbationoftheorbitingplanet(s).DispersedFixed-delayInterferometry(DFDI)isatechniquetoeconomicallyimplementtheRVdetection.Chapter2describesmyinvolvementintwoinstrumentationprojectsbasedonthistechniqueandfocusonamethodofincreasingtheshort-termstabilityoftheinstruments. WheninterpretingthereducedRVdata,peopleusuallyuseaKeplerianmodeltodescribetheorbitsofexoplanets.However,ifthereismorethanoneplanetinthesystem,theKeplerianmodelneglectstheinteractionbetweenplanets.Theoretically,ann-bodysimulationshouldbeusedfordescribingthesystemintermsofaccuracy.However,theN-bodymodelismuchmorecomputationallyexpensive.Inanattemptto 13

PAGE 14

determinewhenafullN-bodymodelisimportant,IcomparetheKeplerianmodelwithN-bodymodelforcharacterizingtheorbitsofasetofsimulatedinteractingtwoplanetsystems,whichisdescribedinChapter3. GivenaRVdata,MarkovChainMonteCarloprovidesaconvenienttooltosamplefromtheposteriordistributionofthemodelparameters.OnechallengethatremainsiscomputingtheMarginalLikelihoodofthePosteriorortheBayesianEvidence,asitrequiresthecalculationofahigh-dimensionalintegral.AsimplealgorithmforcalculatingtheBayesianEvidenceisintroducedinChapter4. Chapter5presentsseveralapplicationsofthealgorithmintroducedinChapter4.Andthesummeryofthisdissertationresearchisgivenintheend. 14

PAGE 15

CHAPTER1INTRODUCTION Inthepasttwodecades,thestudyoftheplanetsaroundstarsotherthanthesun(Extra-solarPlanetorExoplanet)hasemergedtobecomeoneofthehottestdomainsinAstronomyandPlanetaryScience.Weareabletoobservethesignalfromtheexoplanets,eitherindirectlyordirectly.Thankstotheadvanceddevelopmentofenergytechnology,materialengineering,electronicengineering,andcomputingandinformationtechnology,whichhavebeenprovidingtoolsforastronomerstoprobedeeperuniverse,detectweakersignals,andarchiveandanalyzeincreasinglylargerandlargerdatabase.SincethedetectionoftherstexoplanetsPSR1257+12[ 43 ]orbitingapulsar,thenumberoftheknownexoplanetshasreached777planetsin623planetarysystems,105ofwhichhostmorethanoneplanet(multipleplanetsystem)1.Thenumberofknownexoplanetsisstillingrapidlyincreasing.Meanwhile,astheinformationaboutthesedetectedplanetarysystemsaccumulates,theuncertaintyofmodelparametersdecreases. Theexoplanetarysystems,likeoursolarsystem,usuallycontainplanetsorbitingacentralstar(orstellarremnant,andtherealsoexistsomesystemsthataroundbinarystarsorevenpossiblymorecomplicatedsystem).Intermsofsize,mass,orbrightness(reectedorininfrared),theplanetsaresmallcomparedwiththeirhoststars.Therefore,itisdifculttodetecttheplanetdirectly.Becauseaplanetarysystemisboundedbygravity,aplanetorbitingthestarwillcauseperturbationinthelightofthestar. MostmethodsofdetectingexoplanetsarebasedonanalyzingtheelectromagneticradiationsignalobservedfromthestarslikeRadialVelocitymethods,Astrometry 1AsofAug,01,2012accordingtoTheExtrasolarPlanetsEncyclopediahttp://exoplanet.eu 15

PAGE 16

methods,andTransitingmethods,orfromthebackgroundlightsourcelikeMicrolensing,whilesomemethodstrytoobservetheplanetdirectlybyexcludingthelightfromthecentralstarduringtheobservation.Eachmethodhasitsownadvantagesanddisadvantages,whicharesummarizedbelow.Forsomesystems,combiningmultipledetectionmethodscanimprovetheinformationaboutthemassesandorbitalparametersofthewholeplanetarysystem. Inordertodescribeaplanetarysystem,usuallytwomainkindsofcharacteristicsaboutaplanetshouldbetakenintoconsiderationintheobservation.First,thedynamiccharacteristicsdescribetheplanets'motion,includingtheorbitingperiod,major-axisoftheorbit,andeccentricityoftheorbit,etc.Alsoplanetaxialtiltanditsrotationarealsoveryimportantintheplanetarydynamics,butaregenerallyinaccessibleforexoplanet.Second,thephysicalcharacteristicsoftheplanet,includingplanetmass,itsinternalstructuredifferentiation,atmospherecharacteristics,andmagnetospherecharacteristic.Duetotheobservationallimitation,astronomer'sworksillfocusonthemassandorbitalcharacteristics.Theorbitalinformationiscrucialforunderstandingtheformationandevolutionoftheplanetarysystems.Accuratecharacterizationoftheorbitscanprovideconstrainsfortheplanetformationandevolutionmodels.Asthetotalnumberofdetectedplanetarysystemsincreases,betterstatisticsandstrongerconclusioncanbedrawn. Inatwobodyproblem,whereboththecentralstarandtheorbitingplanetareregardedasrigidbody,theKepler'slawsofplanetarymotionpreciselydescribestherelativemovementofthetwobodies.Eachplanet'sorbitisellipticalandisdescribedbyasetofparameters: P,theperiodoftheplanetorbitingthecentralstar,orthea,thesemi-majoraxisoftheorbit,wherePandaarerelatedbythethirdKepler'slaw; e,theeccentricityoftheellipticalorbit; 16

PAGE 17

i,theinclinationofaplanet'sorbitalplane,theanglebetweentheplanetaryplaneandthereferenceplane,whichisusuallyperpendiculartothelineofobserver'ssight; ,thelongitudeoftheascendingnode,whichistheanglebetweenthereferenceplane's0longitudeandtheplanet'sascendingnode,whereascending(descending)nodeisthepointatwhichaplanetcrossesabove(below)itsreferenceplane; !,theargumentofperiastron,theanglebetweenaplanet'sascendingnodeanditsperiastron(theclosestapproachtoitsstar); M0,theorbitalphaseoftheplanetatcertainepoch.Sometimesanotherrelatedparameterisused,forexample,timeofperipassagetp,thetrueanomaly,eccentricanomalyEormeananomalyM; Theratioofthemassoftheplanettothatofthestar.Sincethemassofthestarcanbeesimatedseparately,themassoftheplanetcanbeestimatedaswell.InRadialVelocity(RV)detection,therelevantparameterisK,thesemi-amplitudeofthevelocityofthestarwhichisperturbedbythemotionoftheplanet. InRVdetection,themostcommonlyusedparametersareP,e,K,!,M0(ortp),becauseRVobservationsarenotsensitivetoandthereisadegeneracybetweenmassandi.Furthermore,RVobservationsarealsoaffectedbythewholestar-planetsystemoranylongperiodplanet,eitherofwhichleadstoalinearvelocityconstant.Iftheobservationdataarecombinedfrommultipleobservatories,thereshouldbeanotherconstantvelocityparameterCforeachinstrument. 1.1DetectionMethods Dependingonthephysicsoftheexoplanetarysystem,differentdetectingmethodshaveyieldedsuccess.ThemostestablisheddetectingmethodsareRadialVelocity(RV)methods,astrometry,transitmethods,gravitationalmicrolensing,directimaging,andpulsartiming. RVmethodsarebasedonthephysicsthatwhenaplanetorbitsaroundastar,thestarwillalsomovearoundthecenterofmassofthestar-planetsysteminresponseoftheplanetgravity.Theprojectionofthiswobbleinthelineofobserver'ssightwillcausetheshiftofthestellarspectralduetoDopplereffect,i.e.blueshiftwhenstarmovingtowardstheobserverandredshiftwhenstarmovingawayfromtheobserver. 17

PAGE 18

Theplanet'smassisrelatedtoamplitudeoftheRVwhileotherorbitinformationcanbedeterminedbytheperiodicalchangeoftheRV.Itiseasiertodetectmassiveplanetandshortperiodsplanetthanthosewithsmallmassesandlongorbitalperiods.Furthermore,sinceRVonlymeasuretheinformationalonglineofobserver'ssight,itcannotgiveinformationabouttheinclinationoftheplanetaryplanettotheskyplane,thereforecanonlyyieldtheminimummassoftheplanet(i.e.Mpsini).TheRVmethodisoneofthemostproductivemethodssofar:outofthetotalnumberoftheplanetswehaveknown,475planetsin368planetarysystemsincluding78multipleplanetsystemsweredetectedbythismethod2.FamousRVdetectionprojectsincludeHARPS(HighAccuracyRadialVelocityPlanetSearcher),CaliforniaandCarnegiePlanetSearch,etc. Transitsoccurwhenaplanetcrossesinfrontofthehoststarinthelineofobserver'ssight(i.e.orbitisviewededgeon),causingthebrightnessofthestartodecreaseasmallamountthatdependsontherelativesizesofthestarandtheplanet.Thelowprobabilitythattheplanetshappentobeinthesightofobservationcanbecompensatedbyalargernumberofobservationtargets.TransitsearchprojectslikeCOROTandKeplerhavedetectedmanyexoplanetssystems.Transitobservationsarecapableofmeasuringthediameterofthetransitingplanetaswellasitsatmosphericinformation.Furthermore,theTransittimingvariationmethod(TTV)andtransitdurationvariationmethod(TDV),whicharedenedasthevariationsinthetimingordurationofthetransitofaplanetthathasbeendetectedbythetransitmethod,canprovideanextremelysensitivemethodcapableofdetectingadditionalplanetsinthesystem. Thereareseveralotherpotentiallyeffectiveexoplanetsdetectionmethods.Eachmethodhasitownadvantageanddisadvantageanddifferentsensitivityfordetectingplanetsdependingontheplanetmassandtheorbitalsize.Themostinterestingdomainwheretheexoplanetsexististhehabitablezone(HZ)(moreaccurately,circumstellar 2http://exoplanet.eu 18

PAGE 19

habitablezoneorCHZ),whichmeansthezonewherethedistanceoftheplanetfromthehostingstarandthemassoftheplanetaresuitablefortheexistenceofliquidwateronsuchplanet,whichinferstheconditionsfavorableforlifeonEarth.Lunine[ 31 ]gaveagoodexamplecomparingseveraldifferentexoplanetsearchmethodsinastandardizeddepthofsearchdisplaydevelopedbytheExoplanetTaskForce3.Itcanillustratethediscoverypotentialofplanet-searchsurveysofdifferentsizesandtechniques,overthemassesMandthethermal-equivalentorbitalradiusRER=p Lofexoplanetsbeinglookedfor,whereListheluminosityofthestar.Accordingto[ 31 ],thetransitstechniques,e.g.,theKeplermission,isthemostpromisingtodiscovertheHZexoplanets. AlthoughtheRVmethodsarenotthemostidealforthedetectionofHZplanet,themanyRVprojectswillcontinueondeliveringincreasinglyhighqualitydataasthetechniquesofobserving,calibrationanddataanalyzingadvance.SinceRVsearchstartedtheearliesttoaccumulatedata,theyarecurrentlythemostusefulforlongperiodscandidates.AndwiththedevelopmentofIRRVtechniques,searchingofplanetsaroundlowmassstaralsomakesthediscoveryofHZplanetpossible.Therefore,intheresearchpresentedinthisdissertation,IwillfocusonRVobservationsonly. 1.2RVDetection TheutilizationoftheDopplershifttodetecttheradialvelocityofcelestialbodyisnotnew.Butitwasonlyinthe1990sthattheprecisionofRVmeasurementsmettherequirementfordetectingexoplanets.JupitercausestheSun'svelocitytochangeby12.5ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1overaperiodof12years,Saturninducesa2.7ms)]TJ /F5 7.97 Tf 6.58 0 Td[(1signalwitha30yearperiod,andtheEarthcausesonly0.1ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1variationsoveraperiodof1year.Fortypicalgas-giantplanetstheRVvelocityvariationsareintherangeof1100ms)]TJ /F5 7.97 Tf 6.58 0 Td[(1. 3http://www.nsf.gov/mps/ast/exoptf.jsp 19

PAGE 20

Givenabove,long-termobservationsbystableinstrumentswithaveryhighresolutionarerequired.AtypicalexoplanetRVdetectioninstrumentisamidium/highresolutionspectrograph.Resolutionisrelatedtoresolvingpower: R= =c v (1) whereisthewavelengthofthelightandcandvarethespeedoflightandtheradialvelocity,respectively.Forawavelengthreference,astronomercommonlyuseaniodineabsorptioncelloremissionlampslikeThorium.Byaveragingtheshiftsinmultiplestellarspectrallines,highprecisionofRVmeasurementcanbeachieved.Toincreasetheresolvingpowerofthespectroscopy,insteadofusingsinglestandarddiffractiongrating,highprecisionspectroscopyuseechellegratingforadditionaldiffraction.Forexample,HighAccuracyRadialVelocityPlanetSearcher(HARPS)4usesechellespectroscopytoachieveameasurementaccuracyof0.97ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1. ThereisanothertechniquecalledDispersedFixed-delayInterferometry(DFDI)[ 16 17 ]orExternallyDispersedInterferometry(EDI)[ 7 ],whichsuperposesasetofinterferencefringepatternoverthestellarspectrabyplacinganMichelsoninterferometerinfrontofthegratingintheopticalset.Itcanbeappliedtomediumresolutionspectroscopysinceitcanincreasetheeffectiveresolutionoftheinstrument.Usingthistechnique,theExoplanetTracker(ET)[ 16 ]alreadyfoundaplanetET-1[ 18 ].BesidesupgradingengineeringworkinET,Ihaveparticipatedininstrumentationworkofamulti-objectversionofDFDIinstrument,calledMulti-objectAPORadialVelocityExoplanetLargeAreaSurvey(MARVELS).Inthisdissertation,IwilldescribeamethodsthatincreasetheshorttermstabilityofDFDI. 4http://www.eso.org/sci/facilities/lasilla/instruments/harps/index.html 20

PAGE 21

1.3DataanalysisofRVobservation Thegoalofdataanalysisistounderstandthephysicsmodelsthatproducethedatawehaveobserved.AfterthereductionfromtherawRVdata,i.e.backgroundremoving,ateldcalibration,wavelengthcalibration,barycentricvelocitycalibration,etc,basicallyKeplerianmodelsareusedtotthedata,andtheparametersetthatyieldsthebestgoodnessoft(e.g.,2statistic)givesthebestestimateoftheorbit.Certainlytherearemorecomplicatedprocedurestomakesurethevalidityofthedetection,whicharedescribedinChapter3and4. 1.3.1KeplerianModelv.s.N-bodyModel BasedonKeplerianModel,theRVoftheplanetshoststaris,v(t)=XKp[cos(!p+p)+epcos(!p)]+Cj (1) wherepdenotespthplanet,Kisthesemiamplitudeofthereexvelocity,Cjisthejthreference,andisthetrueanomaly,relatedwithMviaKeplerequation. TheKeplerianmodelisatwo-bodymodelanddoesnotconsidertheinteractionbetweenplanetsiftherearemorethanoneinasystem.Inthatcase,anN-bodymodelshouldbeusedandanN-bodyintegrationshouldbeconductedtosimulatethemulti-objectplanetarysystem.Thetwomodelsareexpectedtogivedifferentorbitalparameterestimationsonlyiftheplanet-planetinteractionissignicant.Inmanyscenarios,theinteractionsbetweenplanetmightbenegligibleiftheyareminorcomparedtothestellar-planetinteraction,forexample,duetothelowmassordistancebetween.Therefore,itisusefultocomparethesuitabilityoftwomodelsforanumberofdifferentsimulatedplanetarysystemsandobservationalprogramstodeterminedeterminewhenafulln-bodymodelisrequired. 1.3.2UncertaintyoftheModelParametersandModelSelection OnetraditionalfrequentistwaytodeterminetheuncertaintyofthettedmodelparametersistheBootstrap:repeatedlyttingthemanysyntheticdatasetstoyielda 21

PAGE 22

numberofsetsofdifferentbestttedparameters.Instead,theMarkovChainMonteCarlomethodprovidesaBayesianwayofcalculatingtheposteriordistributionofthemodelparameters(seeChapter4andreferencesthereinfordetails).TheBayesianmethodsaremoremathematicallyrigorousandcangivetherealposteriordistributionoftheparameters. AsanimportantpartofBayesianinference,themodelselectioninvolveswiththecalculationofthemarginallikelihooddistributionorBayesianevidence.Inlowdimensionthisisalreadycomputationallypractical,butinthecaseoftheorbitofexoplanets,thepotentiallyhighdimensionality(e.g.7foroneplanetmodel,or14fortwoplanetmodel)makesthecalculationoftheintegralextremelycomputationallyexpensive.Manyattemptshavebeenmadetosolvetheproblem.InChapter4,IintroduceanalgorithmtocalculatetheBayesianevidence. 1.4OrganizationofThisDissertation Inthisdissertation,IdescribeanddiscussthreeprojectsthatarerelevanttotheorbitalcharacterizationoftheexoplanetarysystemsusingRVobservations,andthechaptersareorganizedasfollows. InChapter2,acoupleofinstrumentationprojectsaredescribed.Besidescertainlabbuildup,theworkinvolvedtheinstrumentationcontrolutilizingLabView5fortwoexoplanetrelevantdetectioninstruments,i.e.,ETandMARVELS(seeChapter2fordetails).Andamethodofreducingtheobservingnoisebyimprovingtheinstrument'sshorttermopticalstabilityisdiscussed. InChapter3,Ipresentaprojecttoexaminethedifferencesinplanetmassesandorbitalparametersderivedusingtwomodels,theKeplerianmodelandtheN-bodymodel.Keplerianmodelisbeautifullyconciseandcheapintermsofcomputability,sinceitonlyconsiderstheinteractionbetweentheplanetandthehoststar.Moreover, 5Agraphicalcodingprogram:www.ni.com. 22

PAGE 23

itneglectstheinteractionbetweentheplanetsifthesystemownsmorethanoneplanet.Meanwhile,theN-bodymodelincludestheinteractionsbetweeneachpairofthebodiesinthesystemItismoreaccuratethanKeplerianmodel,butithasmuchhighercomputationaldemands.Aplanetdetectionofasimulatedtwoplanetsystemisstudiedinchapter3. InChapter4,Iintroduceanalgorithmtocalculatethemarginallikelihood,orBayesianevidenceofacertainmodel,basedontheMCMCsamplesfromRVobservations.Thiscalculationisusefulwhendecidingwhichmodeltochoose,i.e.,oneplanet,twoplanets,etc,ornoplanet,giventhedatathathavebeenobserved.Iusebothsimulateddataandobserveddatatotestthealgorithmsandcomparewithotheralgorithms. InChapter5,IappliedthealgorithmsdevelopedinChapterIVtoannouncedexoplanetdata,includingsystemswithmultipleplanets.Summaryanddiscussionofmyresearcharegivenintheendofthischapter. 23

PAGE 24

CHAPTER2THECONTROLSYSTEMOFAMULTI-OBJECTRADIALVELOCITYTRACKER 2.1BackgroundandMotivation Sincethediscoveryoftherstexoplanet,theRVdetectionshavebecomethemostmatureandproductivemethodofexoplanetarysystemdetectionandobservationsfortheorbitalparametermeasurements.Theorbitalparametersoftheplanetarysystemsprovidetheconstrainsforthemodelsoftheformationandevolutionoftheplanetarysystems.Inordertoobtainthecomprehensiveinformationaboutthestatisticsoftheorbitalparametersoftheexoplanetarysystems,weneednotonlylargerbutmorehomogenoussamplesofobservingobjects.Thisrequiresanewtechniquestoconductthesurveyofexoplanetarysystemsmoreefciently. Thedispersedxed-delayinterferometer(DFDI)isatechniquecombiningMichelsoninterferometerinserieswithamediumresolutionspectrograph,whichcanbeusedforhigh-precisionradialvelocitysurveysofexsoplanets[ 7 ]andhasthecapabilityofperformingmulti-objectsurveying. TheimplementationofDFDIforsingle-objecthadbeenrealizedbyExoplanetTracker(ET),whichwasinstalledatKittPeakNationalObservatoryandhadobservedin2.1mtelescopeandthe0.9mcoudeFeedtelescope.Itservedtoprovetheconceptwithasuccessfulconrmationoftheknownplanet,51Pegb[ 39 ],andlaterdiscoveredanewplanet,ahotJupiterHD102195b(ET-1)[ 18 ]. Themulti-objectversionofDFDIinstrumentforexoplanetdetectionwasbuiltasW.M.KeckExoplanetTracker(KeckET)atApachePointObservatoryduringthepilotstage,apartofSloanDigitalSkySurveyIII(SDSSIII)andlater(since2008)wasnamedMulti-objectApachePointObservatoryRadialVelocityExoplanetLarge-areaSurvey(MARVELS).Ithasalsodemonstrateddetectionofknownplanets,andhasbeenusedstudybrowndwarfandbinaries[ 12 ][ 29 ]. 24

PAGE 25

Thedevelopmentofnewdatareductionpipelinetoturntherawdataoftheseinstrumentsintomeaningfulresultsefcientlyandreliablyrequiresacarefulconsiderationofprecisionlevelsandsourcesoferror,whichinturnhasrequiredestablishingadetailedunderstandingofthephysicalprinciplesbehindtheinstruments(referthedissertationofvanEykenand[ 40 ]).However,theinstrumentationengineeringwilldeterminetheprecisionlevelsandsourcesoferrorbeyondthetheoreticaldesign.Inthefollowingsections,IdescribetheengineeringworkIparticipatedforbothETandKeckETandfocusonmethodstoincreasetheprecisionlevelsandtoreducethesourcesoferror. IparticipatedthetwoinstrumentationprojectsbasedonDFDItechniques.Particularly,IwasinchargeofthealgorithmsandimplementstostabilizetheMichelsoninterferometer(PhaseLocker)intheDFDIinstruments,andthedesignandprogrammingoftheuserinterfaceusedforcontrollingtheinstrumentduringthecalibrationandobservations.InChapter2,IintroducetheprinciplesoftheDFDItechniquesandthetwoinstrumentsforexoplanetdetections,andfocusonthePhaseLockerandtheuserinterface. 2.2PrincipleofDFDI Fixed-delayInterferometermeansaninterferometerthatoperatesinanarrowrange(m)ofdelayscomparedwithinterferometersinFouriertransformspectrographs,whichscanoverlargedelayrange.DispersedFixed-delayinterferometersmeanstheinterferometerscoupledwithspectrographs.Figure 2-1 givestheillustrativelayoutofanDFDIdevice[ 40 ].AMichelsoninterferometerisplacedinfrontofadispersedevice.OneofthetwoarmsoftheMichelsoninterferometerisxedbyapplyingamirrorononesideofthebeamsplitter,whileanotherischangeablewithamirror(mirror2asinFigure 2-1 )mountedinsuchawaythatitcanbetip-tiltadjusted.Givenacollimatedbeamofmonochromaticinputlightandthemirror2tiltalongtheaxisinthepaper,theoutputlightfromtheinterferometerwilldisplayafringepatterninthedirectionasshown. 25

PAGE 26

Figure2-1. Dispersedinterferometerschematic.ypositioncorrespondstopositionintheslitdirection(directedoutofthepageintheinterferometerschematic),andindicateswavelengthinthedispersiondirection.AtrightweshowA)Outputfrominterferometeralonewithmonochromaticlightinput,andmirror2untilted.B)Thesamewithmirror2tiltedalongtheaxisintheplaneofthepage,asshown.C)Imageondetectorwithmonochromaticlightatveryhighresolution.Onefringedemissionlineisseen.D)Detectorimagewithwhitelightinput.E)Imagewithstellarspectruminput.F)AsforEbutatlowresolution.ReprintedbypermissionfromvanEyken,J.C.,2007,PhDdissertation(Figure2-1),andvanEyken,J.C.,2010,ApJS,189,156(Figure1). 26

PAGE 27

Foraninterferometerwithanopticalpathdifferenced,misthefringeorder,andthemonochromaticwavelengthofinterest:m=d, (2) som=d 2, (2) whichimpliesthatasmallchangeinthewavelengthduetoaDopplershiftwillalsocausethefringeorderchangeifthedremainsthesame.Consideringnon-relativisticDopplershiftcausedbyaradialvelocityv, =v c (2) Equation 2 becomes,m=d v c (2)v=c 2d (2) where=2misthephaseoftheinterferenceorder. Equation 2 infersthattheradialvelocityvcausingDopplershiftcanbemeasuredfromthephaseshift(),whichissinusoidalpatternandcanbeeasilymeasured. Iftheinputlightcontainslargewavelengthrange,therewillbealmostnointerferencefringesvisibleintheoutputpattern,sincedifferentwavelengthswillhavedifferentphaseshiftsoverlappedwitheachother.Thespectrographisusedtodispersethelightoverthedetector,allowingthemeasurementoffringeoverasignicantlymorewidewavelengthrange.Realstellarspectrahavethousandsofstellarlines,andinordertoyieldthebestvelocityprecisionrequiredforexoplanetdetection,thelargestpossiblenumberoflinesshouldbeused. 27

PAGE 28

Figure 2-2 showsschematicsimulatedinterferometercomb,thedispersedinterferometerfringepattern,whichisexpectedfromthedispersionbythemedium-resolutionspectrograph.Intheselected'window',orregionofinterest,thefringesareapproximatelyparallel. Figure2-2. Simulatedinterferometerpatterns.Schematicshowofthedispersedinterferometerfringepattern,assuminginniteresolution.Intheselected'window',orregionofinterest,thefringesareapproximatelyparallel.ReprintedbypermissionfromvanEyken,J.C.,2007,PhDdissertation(Figure2-4),andvanEyken,J.C.,2010,ApJS,189,156(Figure3). Iftheinterferometerdelay(d)isknown,thephasecanbemeasuredinasingleexposureusingatwodimensionalimageandvcanbecalculatedaccordingtoEquation 2 iftheinterferometerisstable,thatis,theinterferometerdelay(d)willnotchangeovertimeateverypositionontheinterferometerplane.However,inreality,dneedcarefullydesignandmeasurement,andtheinterferometer'sstabilityisinuenced 28

PAGE 29

bythetemperature,airpressure,andambientvibration.Foraccurateddesignandmeasurement,onecanreferto[ 32 ]and[ 41 ].Iftheinterferometerisnotstable,thefringeshiftwillwashouttheinterferencepattern,washoutthefringes,reducingthemeasuredfringevisibilityandthevelocityprecision.Thefringevisibilityisdeneasfollows:=Imax)]TJ /F3 11.955 Tf 11.96 0 Td[(Imin Imax+Imin (2) whereImaxandIminarethemaximumandminimumintensitiesofthesinusoidalfringepatternalongtheslitdirection,respectively. 2.2.1InstrumentationBasedontheDFDITechniques IhaveparticipatedintheinstrumentationworkfortwoversionsofDFDIinstrumentsforexoplanetdetection.TherstoneistheExoplanetTracker.Myinstrumentationworkincludes1)Improvingtheinstrumentcontrolsystem,includingthetemperaturemonitor,theiodinecellmotorcontrol,andthephotonmetercontrol.2)Revisingandimprovingthelockingsystemtostabilizetheinterferometer.3)UpgradingtheuserinterfacefortheExoplanetTracker,includingcombiningtheelectronicobservationlogging.ThesecondisMARVELS 2.4 .Myinstrumentationworkinclude1)Improvingtheinstrumentcontrolsystem,includingthetemperaturemonitor,theiodinecellmotorcontrol,andthephotonmetercontrol.2)ImplementationofthelockingsystemforlargeversionofthatofExoplanetTracker.3)DevelopmentoftheberfeedingsystemforMARVELS. Inthefollowingsections,Iwilldescribeeachinstrumentationworkandfocusonthemethodswetestedtoincreasethestabilityoftheinterferometer. 2.3TheExoplanetTracker TheExoplanetTracker(ET),duringitdutybetween2002and2008,waslocatedinthebaseofKittPeakNationalObservatory(KPNO)2.1mdome.TheETteammadethebestefforttomakesurethespectrographstable,airinsulated,builtonvibrationisolatedopticalbenchanditstemperatureonlyrespondsslowlytotheoutsidetemperature.And 29

PAGE 30

atechniquecalledPhaseLockingwasusedtolocktheinterferometerfromphasedrifting. 2.3.1PhaseLocker KPNOETusedapiezoelectrictransducer(PZT)tokeeptheinterferometerfromdriftingduringanexposureandkeepthedelaystable.Piezoelectricmaterialscangeneratemechanicalstressinresponsetoappliedelectricvoltage,enablingtheinterferometerdelaytobemaintainedeveninthepresenceofvibrationandtemperaturedrifts.Oneoftheinterferometermirrors(mirror2inFigure 2-1 )wasmountedonthePZT,whichiscontrolledusingaLabView1program.Asingle-modeberfedwithastabilizedHe-Nelaser(632.8nm)isusedtoilluminateasmallportionofthebeamsplitter.ThefringesformedarerecordedbyavideocameraandareusedastheinputtotheLabViewprogram.ThepositionofthesefringesisthenlockedbyapplyingtheappropriatevoltagestothePZTasshowninFigure 2-3 .Thismethodcaneffectivelystabilizetheinterferometerfringepatternovertheperiodofanobservingrun(7-12days).ToobtainatsandcalibrationlampspectrathePZTvoltageisvariedveryquickly,leadingtothePZTmirrorvibratingandwashingouttheinterferencefringes.ThePZTisturnedoffaftereveryobservingrun,re-initializedanddialedbacktoitspreviouspositionatthestartofanewobservingrun. Figure 2-3 isasnapshotoftheLabviewprogramcontrollingthePZT.Ontheleftisthefringepatternofthereferencelaserbeam,registeredattheverybeginningoftheobservation,andontherightisthefringepatternofthereferencelaserbeammonitoredduringtheobservation.Asshownasthepartintheframe,aregionofinterestischosentodiscardunevenilluminatedpartofthefringefrombythelaserreferencebeam.VoltagestothreePZTtipsarecalculatedbasedonthedifferencebetweentwofringepatterns,whichareevaluatedbythreeparameters,followingthecalculations 1www.ni.com 30

PAGE 31

Figure2-3. TheuserinterfaceofthePhaseLocker.Left:fringepatternofreferencelaserbeam,registeredatthebeginningoftheobservation.Right:fringepatternofreferencelaserbeam,monitoredduringtheobservation.VoltagestothreePZTtipsarecalculatedbasedonthedifferencebetweentwofringepatterns.ThePhase(position)ofthefringes,Densityofthefringes,andTiltofthefringesaretheinformationdescribingthefringes. 31

PAGE 32

methodsbelow.Allcalculationsarebasedontheirimagesafterthenormalizationoftheilluminationprole,whichisrecordedduringtheattakingandvoltagerampingstage. Phase:representsthepositionofthefringesinROI.Averagingcrossthefringedirection(alongtheslitdirection),andttedwithsinusoidalcurve,thephaseatthetopofROIis. Density:representsthenumberofthefringesinROI.Sameprocedurewiththecalculationof,buttheperiodofthettedsinusoidalcurveis. Tilt:representstheanglesofthefringestothehorizonalline.Picktwosymmetricalverticalregionswithinaframe,onefromlefthalfandtheotherfromtherighthalf,eachaveragedcrossthefringedirectionandttedwithsinusoidalcurves,andusecrosscorrelationbetweentwottedcurves,whichyieldsthedifferencebetweentheleftandrighthalfinpixels,. Thesamecalculationsaremadeforbothregisteredpatternssavedattheverybeginningofobservationsandthelivelyrecordedfringepattern.Thedifferencesbetweentwoimagesthencanbeparameterizedas,,and. ThePZThasthreetipstowhichthemirror2wasmounted.Eachtipcanbestressedaccordingtotheinputvoltage(biasvoltage)controlledbythePhaseLockerLabviewprogram,denotedwithV1,V2andV3,andthedevicealsoprovideacontrolvoltagecalledpiston,whichisactuallyapplyingthesamevoltageonthreebiasesatthesametime.Sothemirror'sposturehas3degreesoffreedom.However,evenifwehadfullknowledgeaboutthemechanicalpositionsofthemirrorrelativetothePZTandtheopticalpathlengthalongthevideocameraplane,itisstillverydifculttoknowhowthefringebehavesduetothechangesofV1,V2andV3.Thatis,weneedthecombinationofthevaluesofV1,V2andV3ifwewouldliketoonlydecreaseonevalueof,andwithoutaffectingtheothertwo.Withoutaccuratecombinationparameters,thephaselockingmightundergoneanunexpectedwobbleanduponsuddenshift,itisdifculttoquicklylockthephasebackbecauseofthedegeneracyofthefringeparameterintermsofthebiasvoltages. 32

PAGE 33

ALabviewprogramisusedtocalibratethebiasvoltagescombinationforeachparameter.ToanalyticallyderivetherelationsbetweenViandfringeparameters,wehavetohavetheaccurateinformationoftheopticaldesignoftheinstrument,themechanicaldesignoftheinstrumentespeciallyhowthemirrorismountedonPZT,andthevoltageresponseofeachprobeonthePZT,whichisunnecessarilycomplicated.Sincewhenthephaseislocked,PZTisonlyactiveinaverysmallrange,itisconvenienttoapproximatealocallylinearrelationbetweenViandfringeparameters.Inaveryshorttime(3seconds),averylowvoltage(0.001V)isappliedtoeachbias,andthechangeofthefringeparameterswillhavechangesateachtime:V1V1 =1 +1 +1 (2)V2V2 =2 +2 +2 (2)V3V3 =3 +3 +3 (2) where denotetheunitelementinbiasorfringeparameterspace,Videnotethebiasvoltageappliedeachtime,whichcanbeequaltoV=0.001V,andidenotethephasechangegiventheithbiasvoltage.Sowehave,V266664V1 V2 V3 377775=266664111222333377775266664 377775 (2)266664 377775=266664111222333377775)]TJ /F5 7.97 Tf 6.58 0 Td[(1V266664V1 V2 V3 377775 (2) 33

PAGE 34

Thecalculationcanbedonebytheprograminseconds(settobeprocessedinsecondsincaseoflowresponseofthecontrolsystem),andthebiascombinationforcertainfringeparametersareobtained.Differentgainisappliedonphase,densityandtilttocontrolthespeedofphaselocking,whenthephaseisclosedtolockedposition,thegainscanbeturnedtoverysmallnumbertopreventthefringeoscillationduetooverlocking. Notethereisdegeneracyincalculatingthe.Whenthefringeshiftinbigamplitude,integerwavesshiftwillberegardedasthesamevaluein.Duringtheobservation,thephasedriftwillbecontinuouslymonitoredbythephaselockeranditkeepthetrackofthephasechangesothatitbreaksthedegeneracy.However,atbeginningofeveryobservationrun,aphasewrapmeasurementoverthewholedetectorbytakingtungstenlampwithiodineabsorptionlines. 2.3.2DisadvantagesofthePhaseLocker Thephaselockerdescribedinlastsectionrunsaclosed-loopcontrollooponinterferometer.Itsometimesmakesobservationsdifcult.Themajorpartofthedifcultyliesinthecomplicationofthephaselocking,andtherelativeinstabilityofWindowsystem.Unfortunately,LabviewonlysupportWindowsandpartiallyMacOS,mainlybecausetheOSsystemotherthanWindowshavepoorcompatibilityforhardware,butLabviewheavilyusetheirextensioncardpluggedinPCforinstrumentcommunicate.Intermsofalongdurationsurveyforexoplanets,failureofthePZTitselfcanleadtotheinabilitytocoherentlylinkradialvelocities.Forthemulti-objectversionofET,KeckETorMARVELSpilot,thephaselockerbecomesincreasinglymoredifcultsinceitusemuchlargerbeamsplitterstoaccommodatemorebers,areferencelaserbeamonlysheddingononespotcannotfullyrepresentthefringeinformationoverthewholebeamsplitter. 34

PAGE 35

2.3.3OtherInstrumentationControlandUserInterfaceUpgrading Theobservationprocedurewasdetailedin[ 18 ].Thetungstenlampatswereusedtocorrectpixel-to-pixelvariationsandilluminationvariations,andthenon-fringingTh-Aremissionlampspectrawereusedtodeterminetheslantofthespectraintheslitdirection.Aphotonmultipliertube(PMT)wasalsousedtomonitoranduxfedintotheberfromthetelescopebymeasuringthescatterlightatthegrating.Aniodineabsorptioncellwasusedformonitorthedriftoftheinstrumentaswellasanwavelengthcombwithitsabsorptionlines.Duringtheobservation,besidescontrollingthetelescopeandthephaselocker,theobserverneedcontroltwolamps,onemotortomovetheiodinecell,monitorthetemperatureandux,takeexposureswithdifferentcombinationsoftarget,lamps,iodineinornot,etc.,andwriteobservingnotes.Manyexposurestookjustminutesandinanobservingnighttherecouldbehundredsofexposure.Itwasverycommonforaobservertomakemistakesinsuchscenario,andittooklongtimetotrainnewpeopletolearntheobservationprotocol. Inordertodecreasetheoperationerror,IinitiatedaupgradeofETcontrolprogram.Thewholeprogramcombinedthefunctionofphaselocking,temperaturemonitorandlogging,CCDshutterinterruptedPMTandlogging(todeterminethequalityofexposureanditscentertime),iodinecellsoftwarecontrol,andobservationlogginginasingleuserinterface.Duetotheolddevicesetupofthelampscontrol,theTungstenandTh-Arlampsstillneedtobeturnedon/offviahardware.ButtheETcontrolprogramprovidedtheremindericonsandcheckboxincasetheobserverdidnotcorrectlyoperate.Foraplannedobservationnight,atargetlistneedtobeimported,andtheobservedcanjustselectthetargetandexposuretype.Thecertaincontrolactionwouldbetakenbythecontrolprogram,whichmaketheobservationnightmucheasier. Unfortunately,thisprogramonlyworkedforlessthanyearbeforeitwasreplacedbyaUnix-basedcontrolsystem,whichisabyproductoftheupgradeofMARVELS.And 35

PAGE 36

thePZTphaselockerwasalsoreplacedbymonolithicinterferometerinstrument[ 32 ],ExtremelyHighPrecisionExtrasolarPlanetTrackerInstruments(EXPERT). 2.3.4CurrentStatusofET Currently,EThasbeenreplacedbyEXPERT.EXPERTisacombinationofathermallycompensatedmonolithicMichelsoninterferometerandacross-dispersedechellespectrographforextremelyhighprecisionDopplermeasurementsfornearbybrightstars(e.g.,1ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1foraV=8solartypestarin15minexposure).EXPERTalsohasadirectcross-dispersedechellespectroscopymodefedwith50micronbers.IthasspectralresolutionofR=27,000andasimultaneouswavelengthcoverageof390-1000nm.AndithasbeenworkingonKPNO2.1msince2008. 2.4MARVELSPilotProgram MultiobjectAPORadialVelocityExoplanetLarge-areaSurvey(MARVELS)istherstfull-scale,multi-objectexoplanettrackerinstrumentusingtheSloan2.5mtelescopeatApachePoint,aimingtobecapableofawide-eldradialvelocitysurveyforexoplanetsaroundTYCHO-2starsinthemagnituderangeV7V12.ThepilotprogramwasfromDecember2006toMay2008,alsocalledKeckETatthemoment.ItspilotdesignisbasedontheideaofsingleobjectET,butitcansimultaneouslyobservenolessthan60objects,asshowninFigure 2-4 2.4.1TheTwoPointsPhaseLocker TheprincipleoftheactivephaselockerisbasicallythesameasthatofET.Butduetotheneedoffeedingmoreobjects,thebeamsplitterandthemirroraremuchlargerthanthoseinET.Usingonereferencelaserbeamforphaselockingisnotsufcient.IfthereferencelaserbeamislocatedatbottomofthebeamsplitterandthephasewaslockedbycontrollingthemirrormountedonthePZT,thenitwasmeasuredthat,sometimesthefringesonthetopofthecouldhavelargeshiftsovertheobservationtime.Phasewrapsmeasurementareabletomeasuretheoverallphaseshift,however,thephaselockerusingthebottomreferencebeamfringeisindifferenttotheshiftonthe 36

PAGE 37

Figure2-4. BeamsplitteropticallayoutofMARVELSpilotdesign.ReprintedbypermissionfromvanEyken,J.C.,2007,PhDdissertation(Figure3-2). top.Figure 2-2 illustratesthatthefringesonthetopishigherorderedthanthoseinthebottom.Thephaseshiftduetotheinstabilityoftheinterferometerislargeronthetopthaninthebottom.Thesecondreferencebeamandvideocameraareusedtomonitorthefringepatternonthetopofthebeamsplitter. TheLabviewprogramwasalsomodiedasinFigure 2-5 .Followingthemethodsofparameterizingthefringepattern,wecanalsomeasurethephase,densityandtiltofupperandlowerimagecomparedwiththesavedimage.Butoutofthesesixparameters,onlythreeareindependentsincethedegreeoffreedomisequalto3for 37

PAGE 38

describingtheposition/postureofthemirrorordescribingthefringes.Soweonlyneedtopickthreeofsome.Sincetheparameterphaseismoredirectlyrelatedwithournalphasemeasurementforradialvelocity,twophaseparametersarechosen.Andsincethefringepatternissensitiveonthetop,itisreasonabletopickanotherparameteronthetop.Inourcase,itistilt. ThecalibrationforbiascombinationandphaselockingprocessaresimilartothatlistedintheKPNOET,exceptthatinsteadofusingphase,densityandtiltfromonesingleframe,forMARVELSpilot,thephaselockerusetheparametersphasefromupperandlowerimage,andtiltfromlowerimage.TheobservationprocedureisalsosimilartothatofET.AndthecontroldiagramofMARVELSpilotprogramareshowninFigure 2-6 .Althoughthetwopointphaselockerareprovedtobeeffectivetokeeptheinterferometerstable.TheMARVELSpilotprogramservedasanexploreranddemonstrationofmulti-objectplanetsurveyinstrumentanditexperiencedseveralengineeringupgrade,additionallytheobservingenvironmentwerenotspecicallydesigned,itsdatasufferfrompooruncertainties. Figure 2-7 [ 19 ]showsthemeasurementofthesolarRVbyoneselectedberwiththeTwoPointPhaseLockerequippedMARVELSpilotinstrument.TheuncertaintyoftheRVmeasurementsare9.8ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1,whilethephotonlimitederrorsare6.6ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1.Consideringthatotherpossiblenoisesourcesarelikelytoexistinthepilotinstrument,thisindicatestheTwoPointPhaseLockercansuccessfullystabilizethephase. 2.4.2UpgradeofMARVELS Basedontheconsiderationofthelongtermstability,theMARVELSteamabandonedWindowsOSandturnedtoUnix.Asaresult,theoriginalcontrolprogramwasberecoded.BecauseLabviewanditshardwarearenotcompatiblewithUnix,alloftheinstrumentcontrolwereredesignedandimplementedusingJava.Furthermore,thesuccessofthedemonstrationofmonolithicinterferometerofsingleobjectET[ 32 ] 38

PAGE 39

Figure2-5. Thetwopointphaselocker.Thetwopointsphaselockerhaveamorecomplicateduserinterface.Itmonitortwofringepatterngeneratedbytworeferencebeams,oneonthetopthebeamsplitterandtheotheroneonthebottom.Thefringesparameterarethephaseandtiltofupperframeandthephaseoflowerframe. 39

PAGE 40

Figure2-6. ThecontroldiagramofMARVELSpilotprogram.ET1andET4arethenameoftwoWindows-basedcomputersthatcontroldifferentcomponentofMARVELSinstrument. 40

PAGE 41

Figure2-7. ThestabilityofMARVELSpilotinstrument.Theinstrumentstabilitywasmeasuredbyobservingthesolarspectrumwiththeopticalberspointingtothecloudlesssky,andusethedatareductionpipelinetomeasuretheRVofthesun.Topandbottompanelsshowthemeasurementwithoneselectedber.TheRVmeasurementfrombothbercanobtaintheuncertaintyabout9.8ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1,whilethephotonlimitingerroris6.6ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1foreachber.ReprintedbypermissionfromGe,J.2006,AAS,20919102(Page16,TopFigure). 41

PAGE 42

encouragedtheMARVELSteamswitchedintomonolithicinterferometerintheofcialprogram,whichhavereducedthedifcultyoftheinstrumentcontroldevelopment. 2.5Summary Theengineeringworkdescribedabovearefocusedontheinstrumentcontrol,especiallyonthephaselockermethodsinbothKPNOETandMARVELSpilotproject.MARVELShadbeenrunningforabouttwoyears,anditsRVprecisioncanbearound10ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1for8magnitudestars.Thestabilityoftheinstrumentandthedatapipelinehaveimprovedcomparedwiththepilotinstrument,whichmakesitpossibletousethedataforscience.Actually,Themulti-objectradialvelocitysurveydataofMARVELShasbeenabletoproducescienticresulttostudybrowndwarf[ 29 ]andeclipsingbinaries[ 12 ],etc.Butduetothechallengeinbothsoftwareandhardware,MARVELScurrentlyisnotabletoreachtheRVprecisionforexoplanetdetection.UnderstandingtheoriginofthenoiseinthisinstrumentisalsoveryhelpfulforfutureimprovementoftheDFDItechniques.IneverysurveyplateofMARVELS,usuallythereisatleastonebrightstarknowntohostagiantexoplanetandshowsalargeamplitudeRVcurve.ThesestarswithacceptableprecisionarethereferencestarsusedtocomparewithMARVELStargets.InlaterChapters,IwilluseBayesianmethodstoanalyzeselectedreferencestarsusingMARVELSRVdata. 42

PAGE 43

CHAPTER3KEPLERIANMODELSVERSUSN-BODYMODELSONRADIALVELOCITYOBSERVATIONSOFRESONANTPLANETARYSYSTEMS 3.1BackgroundandMotivation Mostexoplanetsystemswerediscoveredwithradialvelocity(RV)techniquesandcharacterizedusingKeplerianorbitalparametersincludingperiod,P,eccentricity,e,inclinationoftheorbitalplanetotheobservedsky,i,argumentofperiastron,!,longitudeofascendingnodes,,andmeananomaly,M.Keplerianmodelneglectsmutualinteractionsbetweenplanets,soitisonlysuitableforsingleplanetorweaklyinteractingmulti-planetsystem,whenalloftheorbitalelementsexceptMareconstants.Inmulti-planetsystemsthemutualinteractionbetweenplanetsingeneralisnegligibleintermsoftimescaleoramplitude,Keplerianmodelcanbeagoodapproximation.Thereforethereexmotionofthecentralstarcausedbytheplanetsisthesuperpositionsofthereexmotioncausedbyeachsingleplanet'sKeplerianmovement,i.e.,v(t)=XpKp[cos(!p+p)+epcos(!p)], (3) wherepdenotewhichplanet. Whentheplanet-planetinteractioninmulti-planetsystemmodelbecomestrongenoughthattheKeplerianmodelwillfailtopreciselydescribetheorbits,andthemotionofthestardivergesfromthatdescribedwithKeplerianmodel.AtypicalexampleisGJ876withtwoplanetin2:1meanmotionresonance(MMR)[ 33 ]:N-bodytsarerequiredinordertoaccuratelycharacterizetheorbitalparameters.Therefore,itisusefultoquantifytheboundaryofthesuitabilityofKeplerianmodelfordescribingthemultipleexoplanetsystems,andwhenaN-bodymodelisrequired.AsystematicandstatisticalstudycanhelpinterprettheRVdataanddesignmoreefcientobservingstrategies. InChapter3,IpresentaprojectexaminingtheperformanceoftheKeplerianmodelfordescribingmulti-planetsystembycomparisonbetweentheKeplerianmodelandtheN-bodymodelwhentheyareusedtocharacterizetheorbitsofasametwoplanet 43

PAGE 44

systemin2:1MMR.TheMMRcaptureistheresultfromplanet-planetinteractionduringtheevolutionoftheplanetarysystem.Itispredictedtobeacommonoutcomeofdiskmigrationandisobservedinsomeexoplanetsystems,especiallyforthelowerorderMMR(e.g.2:1MMR).resonantangle TheparametercharacterizationisbasedonamockRVdatasetfromasimulatedsystemwithN-bodyintegrations.ByvaryingtheN-bodysimulationparameters,includingmasses,massratiobetweeninnerandouterplanets,anddifferentevolvingstagesduringtheMMR,Iinvestigatedthegoodnessoftsoftwomodels,thebiasanderrorofmeasurementofthettedparameters,indifferentscenario.ThereasonofchoosingRVobservationonlyisthatcurrentlytheRVobservationscovermostknownexoplanetsandsinceitstartedtheearliesttoobservetheexoplanetssystemcomparedwithotherdetectionmethods,itistheonewhichcanprovidethemostinformativedataforthelongperiodexoplanetcandidatesandconstraintheorbitalparametersformulti-planetsystems. Chapter3isorganizedasbelow:InSection 3.2 ,Iintroducethemethodtosimulateatwo-planetsystemindifferentevolvingstagesof2:1MMR,thegenerationofmockRVobservations,andorbitalparametercharacterizationwiththedual-Kepleriantting.InSection 3.3 ,IcomparetheKepleriantswithN-bodyts,andthettedparameterswithtwomodels.InSection 3.4 isthediscussionandconclusion. 3.2Methodology OurrststeptoinvestigatetheorbitalcharacterizationusingKeplerianmodelwithRVobservationdatafrominteractingplanetarysystemsistogenerateatwo-planetsystemusingN-bodysimulation.Byvaryingthesystemparametersincludingtheevolutionstagesaroundthe2:1MMR,themassesandmassratioofthetwoplanets,360systemsaregenerated.Thesecondstep,applyarealobservationtimeseriestotheN-bodyintegrationoutputstomockaRVobservationsforeachofthe360planetarysystems.Inthenalstep,themockRVobservationsarettedwithKeplerianmodel,and 44

PAGE 45

theresultsarecomparedwiththemeasurementfromtheN-bodyintegrationsoutputs,whichareregardedasthetruevalueoftheorbitalparameters.Thedetailofeachstepisintroducedinthefollowingsections. 3.2.0.1GenerationofResonantPlanetarySystem Tosimplifytheproblem,westartthesimulationfromacoplanar,relativelynon-inclinedtwo-planetsystem,i.e.,inclinationsforbothplanetarei=0deg),oninitiallycircularorbits(ei,o=0,hereinafter,i,oareusedtodenotetheinnerandouterplanet,respectively)aroundasolar-masscentralstar(M?=1M).Allotherorbitalparameters,i.e,thelongitudeofascendingnode,,theargumentofperiastron,!andthemeanAnomalyMi,oatstartingepocharesettobe0sincetheRVobservationisnotsensitivetothem.ThemassratiobetweentheinnerandouterplanetsvariedMi:Mo=3:1,2:1,1:1,1:2,and1:3inunitJupitermass(MJ),andthemassofthelessmassiveplanetrangedfrom0.1to1.2MJ. WechoseMERCURY6,ageneral-purposesoftwarepackagetoperformN-bodyintegrationstosimulateresonatethetwo-planetsystems.ItincludesmultipleN-bodyalgorithmswhichisoptimizedforplanetarysystem[ 3 ].ThegeneralBulirsch-Stoer(BS)algorithmisused,whichisthemostaccurate.MERCURYallowsuser'sdenedforcestobeappliedonthesimulatedbodies,whichenabletosimulatetheplanetinwardsmigrationcausedbyplanet-diskinteraction.TheintegrationparameterssettingislistedinTable 3-1 Theinitialsemi-majoraxisofouterplanetissettobeao=0.83AUandtheintegrationstartswithforcedinwardsmigrationontheouterplanet.Theinitialsemi-majoraxisoftheinnerissettobeai=0.5AUandthereisnoforcedmigrationappliedonit.DuringtheN-bodyintegration,astheouterplanetmigratescloseto0.794AUwithmoderateenoughmigrationrate(=)]TJ /F6 11.955 Tf 9.3 0 Td[(10)]TJ /F5 7.97 Tf 6.59 0 Td[(5AU=Year),the2:1MMRcaptureoccursandbothplanetsarelockedintheresonanceandstarttomigrateinwardstogetherwitha 45

PAGE 46

migrationratelowerthanthepresetvalue.Theeccentricityofbothplanets,eiandeoareexcited. Themigrationwillgoonuntiltheintegrationsnishinthepresettime(2107days5.5104years).ThreesamplesofsuchlongN-bodyintegrationsareillustratedinFigure 3-1 :eachpanelshowsatypicalevolvingsystemwithdifferentmassesandmassratioofthetwoplanets.Inthetoppanelofeachplot,besidesthecurvesofsemi-majoraxis,aiandao,acurveofperiodsratioP=Po=Pi,derivedfromao=ai,decreasestoandlocksaround2asthesystemevolves.Thelibrationoftheperiodratioiscausedbythexedmigrationratewhileaiandaodecrease(themigrationratearesupposedtoproportionaltothesemi-majoraxis,butforonlythegenerationofresonantsystemsinsteadofthestudyoftheplanetaryevolution,thissimplicationisacceptable).Theresonanceangles,sin(i,o)inthebottompanelineachplotchangefromradonvaluestolibratingvalue,indicatingthesystemisnotonlywithperiodscommensurability,butindeedinMMR.Theincreaseoftheeccentricitiesofbothplanetsdependsonthemassesandmassratio,butusuallytheinnerplanet'seccentricityrisesfasterthantheouter,exceptwhenthemassratioMi:Mo=3:1. Sixsnapshotsaretakenatdifferentevolvingstagesastheplanetsmigrateinwardsandtheireccentricitiesareexited.Thechosencriteriaofthesixstagesaredescribedasfollow. Immediatelybeforetheresonantcaptureoccurs,denotedwithbr.Giventhemigrationrateof)]TJ /F6 11.955 Tf 9.3 0 Td[(10)]TJ /F5 7.97 Tf 6.58 0 Td[(5AU=Yearandstartingsemi-majoraxisofinnerplanet0.83AU,wechosetheorbitaloutputat1106days,whentheouterplanetmigratedtoao=0.803AU. Immediatelyaftertheresonancecaptureoccurs,denotedwithar.WetaketheMERCURYintegrationoutputat1.66days.Themigrationratevarywiththedifferentmassesandmassratiooftwoplanets,resultinginthedifferenteccentricitiesatthisstages. Deepintheresonance,whentheeccentricityoftheplanetwithhighereccentricityareexitedandreachedthevaluesof0.2,0.3,0.4,0.5,denotedwithe1,e2,e3,e4, 46

PAGE 47

Figure3-1. Samplesofevolvingsystemswithdifferentmassesandmassratio.Themassesandmassratioarelabeledineachpanel.Black:parametersofinnerplanet.Red:parametersofouterplanet.Topplot:thesemi-majoraxisandtheperiodratio(Po=Piindashblackline).Middle:theeccentricitiesandthesinevalueofresonanceangles. 47

PAGE 48

respectively.ExceptwhenmassratioMi:Mo=3:1,theinnerplanethasahighereccentricityanditseccentricityisusedfortheeccentricitychoosingcriteria. Withthesixsnapshots,theorbitaloutputareextractedastheinputforafollowedshortN-bodyintegration.Theshortintegrationtimeisabout30years,withhigherintegrationaccuracyandoutputfrequency.Theforcedmigrationsareturnedoffinordertosimulateastableresonantplanetarysystem.Somesystems,especiallythosewithhighmassesandhigheccentricities,arenotdynamicallystableinlongrun.5systemsareproventohavecloseencounterswith10000yearsN-bodyintegration.Informatof(Mi,Mo,evolvingstage)andinunitofMJ,theyare(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4).Thecloseencounteredsystemareexcludedfromnalresults.Table 3-1 liststheinteractionparametersfortheshortintegrationandcloseencounteredintegrations. Table3-1. MERCURYIntegrationParameters.BSalgorithmisused.Ejectiondistance=100AUandtheradiusofcentralbody=0.005AU ParameterLongShortCloseEncounterCheck ForcedMigrationYesNoNoDuration(days)21071.11044.01106OutputInterval(days)40120IntegrationTimeStep(days)80.24AccuracyParameter10)]TJ /F5 7.97 Tf 6.59 0 Td[(1210)]TJ /F5 7.97 Tf 6.59 0 Td[(1610)]TJ /F5 7.97 Tf 6.59 0 Td[(14 Figure 3-2 illustratesthreesamplesoftheoutputsfromtheshort(around30yrs)integrationsforthreeevolvingstagesofasystemwithMi:Mo=0.5:1.5MJ,usingthesamenotationoflinestyleandcolorwithFigure 3-1 .Asthesystemevolved,ai,odecreasewhileei,oincrease,andthesini,ochangefromrandomvaluetoperiodicallibration. 3.2.0.2GenerationofRadialVelocityDataSets Basedontheoutputofthe30yearshortN-bodyintegrationsataoutputrateofonedatapointperday,thereexmotionofthecentralstararecalculatedbasedonthevelocitiesofinnerandouterplanetsinBarry-centriccoordinates.Weapplythecoordinatesrotationforthereexmotion,bysettingi,o=0sinceitwillnotchange 48

PAGE 49

Figure3-2. Samplesofevolvingsystemsindifferentevolvingstages.WithMi:Mo=0.5:1.5[MJ]indifferentevolvingstages,immediatelybeforetheresonance(br),inMMRwithei=0.2(e1),andinMMRwithei=0.5(e4).Black:parametersofinnerplanet.Red:parametersofouterplanet.Topplot:thesemi-majoraxisandtheperiodratio(Po=Piindashblackline).Middle:theeccentricitiesandthesinevalueofresonanceangles. 49

PAGE 50

thepatternofRVdata,theinclinationofthesystemii,o=90(edge-on),andthevelocityofthecenterofthemassoftheplanetarysystemsas0.ThenweinterpolatethetransformedRVvelocitiestoarealobservationtimeswith150datapointsover15years.Consideringatypicalvalueofground-basedRVobservationuncertainty,wethenaddaGaussiannoiseof=2ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1tothepureRVdata.The15yrsobservationscovermorethan10orbitalperiodsevenforthesystemwiththelongestperiod,andtheobservationtimesareevenlydistributed. NotewegenerateRVdatausingBarrycentriccoordinatesfromtheoutputsofN-bodyintegrations,butuseJacobbicoordinatestoextractthevaluesoftheorbitalparameters.Inperturbedplanetarysystem,theorbitalparametersareosculatingratherthanapparent.BychoosingtheJacobbicoordinates,thevariationsoftheorbitalparametersarelessthanthoseinothercoordinates[ 5 ].Inthecaseofoursimulatedsystem,Jacobbicoordinatesregardthecentralstarastheoriginoftheinnerplanetforthemeasurement,andthecenterofmassesofinnerplanetandthecentralstarastheoriginofouterplanetfortheorbitalmeasurement.IfwesayKeplerianmodelcanappropriatelycharacterizetheresonantsystems,thebestcoordinatechoiceshouldbeJacobbi.SoweaveragedtheJacobbicoordinatesoutput(theoneperdaydatabeforeinterpolationtotherealobservationtimes)overtheobservedtimespanforthetruevalueoforbitalparameters,andfollowedtheformulainJacobbicoordinatestocalculateindirectparameterslikeKi,oandPi,o,tobeusedforinitialguessofKepleriantandttingresultscomparisons. Wedisplaythedistributionoforbitalparameters,Pi,oandei,o,inFigure 3-3 ,forallshortiteratedsystemsindifferentevolvingstages.ThevaluesarebasedonaveragedJacobbiorbitalparameters. 3.2.1Multi-planetKeplerianFitting WeuseRVLINby[ 44 ],amultiplanetKeplerianttingpackagetondthebest-tteddual-KeplerianparametersofeverysimulatedsystemwithmockRVdata.Theresultcan 50

PAGE 51

Figure3-3. ThedistributionofPi,Po,eiandeo.Jacobbiorbitalparametersatdifferentstatesofthemigration,fromthelefttotheright,ar,e2,ande4arelledhistograms,andbr,e1ande3areopenhistograms. 51

PAGE 52

becomparedwiththetruevalueoftheorbitalparametersmeasuredduringtheN-bodyintegrations.BasedonKeplerianModel,theRVoftheplanethostingstarisv(t)=Xp=i,oKp[cos(!p+p)+epcos(!p)]+Cj (3) wherepdenotesi,oforinnerorouterplanet,Kisthesemiamplitudeofthereexvelocity,Cjisthereferencevelocities,andisthetrueanomaly,relatedwithMviaKeplerequation.Becauseinthedatasimulationdoesnotinvolveanychangingreference,Cjissupposedtobe0. Thettingpackageisbasedonlocalminimumsearchingalgorithm,soasetofgridsofinitialguessesofplanetorbitalparametersareusedtondtheglobalminimuminthegoodness-of-t(2statistics)space.AssumingthecorrespondingbestttedorbitalvalueshouldbeclosetothetruevaluemeasuredfromtheN-bodyintegrations,thetruevaluesareusedtobethecentersoftheinitialguessesgrids.Thetsgivingtheminimum2amongthegriddedinitialguessesgivesthebestttedparametersforthemockRVdataset.AndforeachRVdataset,wescramblethenoise,andrepeatthegriddedtting30times.Thenthestatisticsofthets(mediananduncertainty)ofeverysystemareobtainedwithBootstrap(N=30). Theinitialmodelparameterguessesincludeperiod(Pi,o),eccentricity(ei,o),timeofpassageofperiastron(tp,i,o),amplitudeofRV(Ki,o),andargumentofperiastron(!i,o),aswellasthevelocityreference()andalongtermslope(_v).Thevalueofand_varesettobe0.WendthequalityofthettingisnotsensitivetotheinitialguessesofKi,oand!i,o.Wealsondthattheoverallminimumisnotsensitivetotheinitialguessforperiodsaslongastheirtruevaluesareincludedinthegrids.Inordertoreducethecomputation,P,Kand!areallsetasthetruevalue.Weonlyuseei,oandtp,i,oaroundtherealvalueforgriddedinitialguesses.Foreachplanet,wetried8differente's,whicharecenteredwiththetruevalueesim,startingfrom0,anduseagridspacingofaquarterofthetruevalueofei,o.Similarlyourgridsinclude8differentti,o,witha 52

PAGE 53

spacingofaquarteroftherealvaluePi,o.ForeverysetofRVdata,thereare84=4096initialguessesgridpointstherefore4096ts,andtheonewiththeoverallminimum2isregardedastheglobalminimumt.AndwithBootstrap,N=30globalbesttsyield30bestttedvaluesforeachorbitalparameters.Foreveryorbitalparameter,weusethemedianfromthe30tsasthereportedtsforthegivenplanetarysystemand34%quantilefromthemedianastheerrorofthets. Figure 3-4 givesanexampleofpureRVdatageneratedandinterpolatedwithrealobservationtimesinthetoppanel,andanoisedRVdatawithdual-Keplerianttinginthebottompanel.ThesystemhasMi:Mo=0.5:1.5MJandinthestageofe4,indicatingthesystemisindeepMMR.ThereisanobvioussignatureofMMRfromthepureRVcurve.ThestarsaretheinterpolatedRVdata,expanding15yrswith150points.Inthebottompanel,thestarsaretheRVdataaddedwithGaussiannoise=2ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1;thesolidlineisthedual-Keplerianttedcurve. 3.3ResultsandAnalysis:ComparisonBetweenKeplerianandN-bodyModels WecancomparehowwelltheKeplerianortheN-bodymodeldescribetheorbitalparametersforeachsimulatedsystembasedonRVobservations.Weusetwomethods:1)Comparethegoodness-of-toftheobservationwiththetwomodels,and2)comparethebestttedorbitalparametersforthetwomodels. 3.3.1GoodnessofFit Tocomparethegoodness-of-toftheRVdataforeachsystemwithtwomodels,wecalculateFstatisticbetweentheKepleriantsandtheN-bodyts.TheFstatisticisF=RSSK=Df,K RSSN=Df,N, (3) whereRSSKandRSSNaretheresidualsumofsquareofKeplerianmodelandN-bodymodel,andDf,KandDf,Narethedegreeoffreedomforeachmodel. TheRSSNoftheN-bodymodelforeachRVdatasetisthesumofsquareoftheaddednoisetothedataset.AndthedegreeoffreedomfortheN-bodymodel 53

PAGE 54

Figure3-4. AnexamplesystemwithKepleriants.Mi:Mo=0.5:1.5[MJ]ine4.Toppanel:IdealRVdatageneratedfromthesimulatedsystem.SolidlinepureRVfromshortN-bodyintegrationaftercoordinatesrotation;starsinterpolatedRVdatatorealobservationtimes(150observationsover15yrs).Midpanel:MocknoiseRVdata.StarsnoisedRVobservationswith=2ms)]TJ /F5 7.97 Tf 6.58 0 Td[(1;solidlinedual-Kepleriantwithmedianglobalminimum2basedonBootstrap(N=30).Bottompanel:FitresidualandarticialnoiseaddedformocknoisedRV.OpencirclesresidualsofN-bodyt,whicharethenoisesaddedtogeneratethemockRVdata;starsresidualsofdual-Kepleriants.Forthissystem,theresidualoftheKeplerianmodelgiveslowqualityts. 54

PAGE 55

Df,N=150)]TJ /F6 11.955 Tf 11.95 0 Td[(1.ForthetwoplanetKeplerianmodel,weuse2tocalculateRSSK,RSSK=22, (3) whereisthesimulateduncertaintyoftheRVobservationand=2.0(m=s). SinceduringthetwoplanetKeplerianttingfor150datapoints,wesetfree5p+1=11parameters,sothedegreeoffreedomfortwoplanetKeplerianmodelDf,K=150)]TJ /F6 11.955 Tf 12.54 0 Td[(11=139.Therefore,theFstatisticsinEquation 3 followstheFdistribution,F(150)]TJ /F5 7.97 Tf 6.59 0 Td[(11,150)]TJ /F5 7.97 Tf 6.58 0 Td[(1). InFigure 3-5 ,weplotthedistributionsofresidualsumofsquare(RSS)oftheN-bodymodelts(toppanel)andthatoftheKeplerianmodelt(bottompanel)forRVobservationsofallofthesimulatessystems.ThedistributionoftheN-bodyRSSappearsnearlyGaussianwithameanof600,whichisexpectedsincethereare150datapointsandtheuncertaintyofthesimulatedobservationnoise=2ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1.TheRSSofthetforallsystemswithKeplerianmodelisscaledbydividingtheratioofthedegreeoffreedombetweentwomodels:Df,K Df,N=139 150)]TJ /F5 7.97 Tf 6.59 0 Td[(1.ThedistributionoftheRSSwithKeplerianmodelspanfromaconcentrationaround600toover1000(notincludedintheplots).Thisindicatesthatintermsofthegoodness-of-t,theKeplerianmodelresultsinasimilarqualitytastheN-bodymodelinsomecases,butinothercases,thetwomodesgiveverydifferentqualityts.SinceweregardthetoftheN-bodymodelistherealtofthedata,theKeplerianmodelfailstosuccessfulttheRVdataforsomesystems. Figure 3-6 showsthecomparisonbetweentwomodelsindifferentscenariowiththecontoursofthepvaluesofFstatisticintroducedinEquation 3 .Differentpanelsrepresentdifferentmassratio.AlongthexdirectionthetotalmassoftwoplanetincreasesandyrepresentsthestateofthesystemasitmigratesfrombeforetheMMRcapture(bottom)todeepinMMRwithlargeeccentricity(top).Smallerpvalue(darker 55

PAGE 56

Figure3-5. DistributionofRSS(ResidualSumofSquare)ofN-bodytsandKepleriants.Toppanel:distributionoftheRSSoftheN-bodyts,i.e.,thedistributionofthenoiseaddedtotheidealRVdatatogeneratethemockRVdata.Bottompanel:distributionoftheRSSoftheKepleriant,scaledbytheratioofthedegreeoffreedomofthetwomodels,i.e.,RSSK150)]TJ /F5 7.97 Tf 6.58 0 Td[(1 150)]TJ /F5 7.97 Tf 6.59 0 Td[(11. 56

PAGE 57

color)impliesthatforthecorrespondingsystemitismorelikelytorejecttheKeplerianmodelovertheN-bodymodelasanappropriatemodel. ThereissharpboundarytodeterminewhethertheKeplerianmodelshouldberejected.Ifwechoosea5%signicancelevelfortheFtest(correspondingapvalueof0.05),avalueofthetotalmassoftwoplanetaround2.5MJcanberegardedasathreshold.Forsystemwithtotalmassoverthisvalue,theKeplerianmodelshouldberejectedbasedontheF-test.Whentheouterplanetismoremassive,thatis,forthesystemswiththemassratioof1:2and1:3,thisthresholdispushedtosmallervaluesofthetotalmass.Andthereisatrendbetweentheevolvingstagesandthetotalmassthreshold. WecansummarizethatforatwoplanetsystemevolvingintoMMR,whentheouterplanetismoremassive,therejectionthresholdoftheKeplerianmodelisafunctionoftotalmassandevolvingstages,andwhentheinnerplanetismoremassive,therejectionthresholdoftheKeplerianmodelisafunctionofonlythetotalmass.Itisbecausethatwhentheouterplanetislessmassive,itsreexmotiononthecentralstarissmallerthanthatoftheinnerandmoremassiveplanet,sotheRVsignalofthecentralstarisdominatedbythechangecausedbytheinnerplanet.EvenindeepMMR,thesignalcausedbytheplanet-planetinteractionissmallenoughthattheRVcanstillbeapproximatedwithKeplerianmodel.However,whenthetotalmassofthetwoplanetincreasesoverthethreshold,theplanet-planetinteractionsarenolongernegligible.Ontheotherhand,whentheouterplanetissignicantlymoremassive,thereexmotionofthecentralstarcausedbytheinnerplanetdoesnotdominate,sotheeffectoftheplanet-planetinteractionsarerelatedwithtotalplanetmassandthestateofthesystemsmigratingfrombeforeMMRcapturetodeepinMMR. 3.3.2ParameterCharacterization InthissectionweconsidertheparametercharacterizationusingtheKeplerianmodel.First,wedenetwofunctions,thebiasandtheerrorofthemeasurementsofthe 57

PAGE 58

Figure3-6. ContoursofthepvalueofF-testbetweenthegoodness-of-toftheKeplerianmodelandtheN-bodymodel.Eachpanelrepresentsdifferentmassratio.xcorrespondsthetotalmassoftheinnerandouterplanets.ycorrespondstheevolvingstages,asdeeperinMMRwithincreasingy.Note5datapointsareexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4).Theregionswithdarkercolorindicateslowerpvalues,thereformorelikelytorejecttheKeplerianmodel. 58

PAGE 59

orbitalparametersPi,Po,ei,eo,Ki,Ko.Forexample,forattedvalueoftheperiodoftheinnerplanetPiforasimulatedsystem,therelativemeasurementbiasofPiis,RelativeBias(Pi)=Pi,t)]TJ /F3 11.955 Tf 11.95 0 Td[(Pi,N)]TJ /F4 7.97 Tf 6.58 0 Td[(b Pi,N)]TJ /F4 7.97 Tf 6.59 0 Td[(b, (3) wherePi,tisthebestttedvalueofinnerplanetperiodwithKeplerianmodel(medianoftheBootstrapvalues,wherefortheBootstrap),andPi,N)]TJ /F4 7.97 Tf 6.59 0 Td[(bisthetruevalueofinnerplanetperiodsmeasuredduringtheN-bodysimulation. AndtherelativemeasurementuncertaintyofPiis,RelativeError(Pi)=^Pi,t Pi,N)]TJ /F4 7.97 Tf 6.58 0 Td[(b (3) where^Pi,tisthe34%quantilefromthemedianofthettedvalueoftheinnerplanetperiodbasedonN=30syntheticdatasetgeneratedfortheBootstrap.Itisequivalentto1ifthedistributionBootstrapoftheparameterisGaussian. 3.3.2.1Pi,o Figure 3-7 and 3-8 displaytherelativebiasanderrorofthemeasurementoftheinnerplanetperiod.Therelativemeasurementerrorissmall10)]TJ /F5 7.97 Tf 6.59 0 Td[(5)]TJ /F6 11.955 Tf 12.08 0 Td[(10)]TJ /F5 7.97 Tf 6.58 0 Td[(3.Asthetotalplanetmassincreases,therelativeerrordecreases.TherelativemeasurementbiasofPishowsdifferentpatterndependingonwhichplanetismoremassive.Iftheinnerplanetismoremassive,Piisunderestimated,andthemagnitudeoftherelativebiasincreasesasthetotalmassofthesystemincreases.Iftheouterplanetismoremassive,thenPiistypicallyunderestimatedasthereisaweektrendwithincreasingbiasasthetotalplanetmassincreases,butthereismuchmorescattercomparedtosystemswithamoremassiveinnerplanet. TherelativemeasurementerrorofPoisslightlargerthanthatofPi,asinFigure 3-10 .Thisisexpectedsincetheouterplanethasalowersignaltonoiseratio.ComparedwithPi,Potendstobemoreunderestimated,asshowninFigure 3-9 ,withsomevalues 59

PAGE 60

evenreaches)]TJ /F6 11.955 Tf 9.3 0 Td[(1%.Themagnitudeoftherelativebiasalsoincreasesasthetotalplanetmassincreases. Overall,themeasurementofperiodsofthetwoplanetsareunderestimatedintheorderof0.1percent,althoughtherelativeerrorofthemeasurementisverysmall(0.01)]TJ /F6 11.955 Tf 12.27 0 Td[(0.1%).Theunderestimationincreasesasthetotalmassincreases.Thiscanbecausedbytheorbitsprecessioncausedbytheplanet-planetinteractions.Notethatin'br'stagethemeasurementoftheperiodsaremorebiasedthaninotherevolvingstageseventhoughwefailtorejecttheKeplerianmodelbasedontheF-test.Thisismightbeduetotheloworbitaleccentricitiesofthetwoplanetsinthisstagesresultinginmorerapidprecession. 3.3.2.2ei,o Figure 3-11 3-12 3-13 and 3-14 showtherelativemeasurementbiasanderrorofeiandeo,respectively. Therelativemeasurementerrorsofei,oclearlydecreaseswithincreasingthetotalplanetmassandmoremigrationfollowingMMRcapture.Itisunderstandablesinceitcoincideswithincreasingsignaltonoiseratio. Themeasurementbiasesofei,oarelargeforbrstage,mainlyduetothenear0valueoftheorbitaleccentricities.Whenthetotalmassislow,themeasurementbiasesofei,oarelarge.Asthetotalmassincrease,iftheinnerplanetismoremassive,themeasurementofei,obecomemoreaccurate,butiftheouterplanetismoremassive,themeasurementofei,ocannotgiveaccurateestimation. 3.3.2.3Ki.o Figure 3-15 3-16 3-17 and 3-18 showthemeasurementbiasanderrorofKiandKo,respectively. BothKiandKoareoverestimated.Thebiasisbiggerforthemoremassiveplanet,whichisreasonablesincemoremassiveplanethasstrongereffectonthereexmotionofthecenterstar.Asthetotalmassincreases,themeasurementsofbothKiandKo 60

PAGE 61

Figure3-7. RelativemeasurementbiasofPi.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 61

PAGE 62

Figure3-8. RelativemeasurementerrorofPi.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 62

PAGE 63

Figure3-9. RelativemeasurementbiasofPo.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 63

PAGE 64

Figure3-10. RelativemeasurementerrorofPo.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 64

PAGE 65

Figure3-11. Relativemeasurementbiasofei.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 65

PAGE 66

Figure3-12. Relativemeasurementerrorofei.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 66

PAGE 67

Figure3-13. Relativemeasurementbiasofeo.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 67

PAGE 68

Figure3-14. Relativemeasurementerrorofeo.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 68

PAGE 69

Figure3-15. RelativemeasurementbiasofKi.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). becomemoreandmoreaccurate.NotethisisintheregionwheretheKeplerianmodelisrejected. 3.4SummaryandDiscussion InthelastSection 3.3 ,onlythemeasurementsofPi,Po,ei,eo,Ki,andKoarepresented.Theotherttedparameters,tp,i,oand!i,oarenotshown,becausewendthatthesevaluesaredifculttoyieldaccurateorprecisevalueswhenusingthe 69

PAGE 70

Figure3-16. RelativemeasurementerrorofKi.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 70

PAGE 71

Figure3-17. RelativemeasurementbiasofKo.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 71

PAGE 72

Figure3-18. RelativemeasurementerrorofKo.Differentpanelrepresentsdifferentmassratio.Ineachpanel,thexaxisisthetotalmassofthetwoplanetsandthedifferentevolvingstageisdrawnwithdifferentlinestyle.Note5datapointsshouldbeexcluded:(3.3,1.1,e3),(3.3,1.1,e4),(3.6,1.2,e3),(3.6,1.2,e4),and(1.2,3.6,e4). 72

PAGE 73

Keplerianmodel.Themostimportantreasonisthatfortheinteractingtwoplanetsystem,theprecessingorbitsmakeitverydifculttodeterminetheorbitalpositionparameters.ItmightbehelpfultousetheprecessingKeplerianmodel,thatis,v(t)=XKp[cos(!p+p+_!pt)+epcos(!p)+_!pt]+Cj, (3) where_!pistheprecessingrateoftheptheplanet. Thisresearchfocusesontheorbitalcharacterizationofthemultipleplanetsystemsindifferentcasesinsteadofstudyingtheformationoftheresonantsystems,whichwouldrequiremorecomplicatedconsideration.Wechosetoconsiderthesystemswithacentralstarofonesolarmass,withtwointeractingplanetevolvingin2:1MMR,whichhaveatotalmassoftwoplanetsrangingfrom0.2to4.8MJ,andthesemi-majoraxiswithin0.8AU.Thesimulationoftheinwardsmigrationoftheplanetsystemisimplementedbyapplyingtheforcedmigrationonlyontheouterplanetleadingtotheresonancecapturewiththeinnerplanet.Thesesimulationsuseasimplieddescriptionformigrationmechanismonlyforthepurposeoftherapidlyproducingoftheresonantsystemstobeusedasinitialconditionsfororbitaltting.Themigrationrateappliedtotheouterplanetwasadjustedsothatmigrationsimulationiscomputationallypredictedyetnosofastthatthesimulatedsystemswereviolentlyinstable.Everysystemwastestedforlong-termstabilityaftertheforcedmigrationwasturnedoff,butbeforeitwasusedasinitialconditionsforshortN-bodyintegrationthatgeneratedthemockRVobservations. TheresultspresentedinSection 3.3 showthatforourchoiceofparameters,theKeplerianmodelshouldberejectedifthetotalplanetmassexceeds2.5MJ,ifweuseasignicancelevelof0.05.Whentheouterplanetismoremassive,thethresholdofrejectingtheKeplerianmodelisalsoafunctionoftheextentofmigrationfollowingtheresonantcapture.Ourtestsimulatedanobservationuncertaintyof2m=s.Thisimpliesthatwithsuchobservationprecision,RVobservationsareoftenunabletodetectthedifferencebetweentheKeplerianmodelandtheN-bodymodel,atthesignicancelevel 73

PAGE 74

of0.05.Itishelpfultoapplytheresultonknownmulti-planetarysystemstodetermineiftheN-bodymodelisrequiredforthestudyoftheorbitalcharacterization. InTable 3-2 ,Ilistedsomesystemsthathavebeenconrmedwithtwoormorethantwoplanets(e.g.,GJ876)butonlythetwoplanetsinorcloseto2:1periodcommensurabilityarelisted.Basedontheinformationofthesystemandtheresultfromthisresearch,everysystemrecommendstheN-bodymodeltostudytheorbitalparameters.Fromthe126multipleplanetsystems,Iexcludedthelongperiodplanetsorplanetwithsignicantsmallmasscomparedtootherplanets.AsimpliedbytheresultinSection 3.3 ,forsystemswithamoremassiveinnerplanet,thetotalplanetmassisthemostimportantfactortodeterminethenecessityofusingN-bodymodel,forsystemswithmoremassiveouterplanet,boththetotalplanetmassandtheextentofmigrationfollowingtheresonancecapture(canbeindicatedpossiblebytheeccentricities)areimportant. GJ876isasystemwellstudiedfornoneligibleplanet-planetinteraction,thereforerequiresN-bodymodelstudy.NSVS1425(AB)b,c,givenitstotalplanetmass>10MJwithamoremassiveouterplanetanditscentralstarwithhalfofthesolarmass,theN-bodymodelisrecommendedfortheorbitalstudy,eventhoughttheperiodsarequitelong.HD82943b,careasystemsimilartothesimulatedsystemsinourresearch,exceptthattheperiodsofthetwoplanetdoubled.Butduetothemoderateeccentricitiesofthetwoplanets,theN-bodymodelisstillrecommended.Kepler-25b,candKepler-31b,caresystemsjustdiscoveredbyKeplerusingtransittechniques.Althoughwehavenoinformationabouttheeccentricities,thetotalplanetmassislargeenoughthatwehavetouseN-bodymodeltostudytheseplanets. Overall,themeasurementofPi,oareunderestimated,whileei,oandKi,oareoverestimated.Whenthetotalmassincreases,therelativeerrorofmeasurementsalwaysdecreases,andtherelativemeasurementbiasesofPi,oincreases,whilethoseofei,oandKi,odecrease.Inotherwords,intheregionwheretheKeplerianmodelisnot 74

PAGE 75

Table3-2. Selectedknownplanetarysystemsaround2:1MMRa.Somesystemshavebeenconrmedwithmorethantwoplanets(e.g.,GJ876),buthereonlythetwoplanetsin2:1MMRarelisted.Thenameoftheplanetisdenoteasstarname+letter,e.g.,GJ876b.iandoareusedtodenotethepropertiesoftheplanetastheinnerortheouterplanet.Basedontheinformationofthesystemandtheresultfromtheresearch,everysystemisrecommendedtousetheN-bodymodeltostudytheorbitalparameters. SystemM(M)Mi(MJ)ai(AU)Pi(day)eiMo(MJ)ao(AU)Po(day)eoN-body? GJ876b,c0.3342.280.20861.120.030.7140.1330.10.256YesNSVS1425(AB)b,c0.5282.81.91276.08.02.92506.00.52YesHD82943b,c1.181.751.19441.20.2192.010.746219.00.359YesKepler-25b,c1.2212.70.0686.23854.160.1112.7204YesKepler-31b,c1.216.80.1620.86134.70.2642.6318Yes ahttp://exoplanet.eu 75

PAGE 76

rejected,themeasurementsofei,oandKi,oarebiased,whileintheregionwheretheKeplerianmodelisrejected,themeasurementsofei,oandKi,otrendtobeunbiased,whichcouldbeduetotheincreasedsignaltonoiseratio.Furthermore,whentheinnerplanetismoremassive,themeasurementsoftheorbitalparametersaremoreaccurate. OurmockRVobservationscovermorethan20periodsfortheplanetsintheN-bodymodelsimulatedsystem,andwith150datapoints,itissufcienttoconstraintheKeplerianorbitalparameters.Ourstudyimpliesthatforverywellobservedinteractingmulti-planetsystems,e.g.in2:1MMR,theKeplerianisnotabletosuccessfullytsomesystemoraccuratelycharacterizetheorbitalparameters.Forothersystems,thesignaltonoiseratioisnothighenoughtotellthedifferencebetweentheN-bodymodeltandtheKeplerianmodelt.Further,themeasurementoftheorbitalparametersexcludingperiodcouldbebiasedwithlargeuncertaintyifweusetheKeplerianmodelforinteractingsystems.Therefore,iftheobservationprecisionwereimproved,itisusefultoconsiderN-bodymodeltsforsystemswithinteractinggiantplanets.Whentheinnerplanetislessmassivethantheouterplanet,usingtheN-bodymodelisevenmoreimportant.Forsomeorbitalparameters(i.e.,tp,!,orM0)intheinteractingplanetarysystem,itishelpfultoconsideratleastthemodiedKeplerianmodel(i.e.,precessingKeplerianmodel)inordertoobtainaccurateestimation. 76

PAGE 77

CHAPTER4BAYESIANMODELSELECTIONONEXOPLANETSRVDETECTION 4.1Motivation Modelselectionbasedongivenobservationdataisacommonprocessinthestudyofthemostsuitablephysicalmodeldescribingtheobservedphenomenon,whicharesupportedbytheobservations.Particularly,fortheexoplanetdetectionandcharacterization,themodelselectioncananswerquestionslike:howmanyplanetsareorbitingtheobservedstar?Ifthenewdetectedplanetonlyaresultfromtheobservationnoise?Bayesianmodelselection,onthecontrarytotheFrequentistmodelselection,usesonlytheobserveddataandmorerigorousmathematicalbase.ItinvolvescalculatingtheBayesfactor(seeSection 4.3 )orcomparingtheevidenceofeachmodel. InChapter3,IdiscussedthemodelselectionbetweenKeplerianmodelandN-bodymodel,andtheirsuitabilityindifferentscenario.InChapter4,Iwilldiscussthemodelselectionbetweenthemodelswithdifferentnumberofplanets,ontheassumptionthatKeplerianmodelissufcientlyaccuratetodescribetheplanetarysystemobserved.Thatis,forhowmanyplanetsistheresufcientevidencebasedontheRVobservationstojustifyaplanetdiscovery. ThecalculationoftheBayesfactororthemodelevidencesishighlycomputationallyexpensiveduetothehighdimensionmodelparameter.Ford[ 13 ]andFord[ 14 ]havedevelopedamethodusingMarchovChainMonteCarlo(MCMC)togeneratesamplesfromtheposteriorprobabilitydistributionfunction(PDF)fortheorbitalparametersoftheexoplanetarysystems,andtheseMCMCsamplesweresuccessfullyusedtoestimatetheuncertaintyoftheorbitalparameterestimationsaccuratelyandefciently.GiventhecomputationcostofthegenerationofMCMCsamples,itisreasonabletousethemformodelselection.BasedontheMCMCsamples,thecalculationofthemodelevidencehasbeenreducedtoaproblemofcalculatinganintegralofamulti-dimensional 77

PAGE 78

functiongivenasetofsamplesdrawnfromthisfunction.Intheexoplanetmodel,thehighdimensionintegrationmakesthecalculationextremelydifcult.Thereareseveralpreviousstudies(seeSection 4.4 )onthisproblem,mostofwhichassumetheposteriorPDFcanbeapproximatebymulti-variantGaussian.AlternativealgorithmsofcalculatingthemodelevidencewithoutassumingtheshapeoftheposteriorPDFbasedontheMCMCsamplesthereforearedesired.InChapter5,IconstructtheFractionEstimator,whichdoesnotrelyonanyapproximationoftheshapeoftheposteriorPDF,anditsestimationuncertaintyiscontrollable.Besidesthemathematicalconstructionoftheestimator,IapplytheFractionEstimatoronrealRVdatasetofaoneplanetsystem,whichwaspreviouslystudiedfortestingothermodelevidenceestimators. 4.2Background InChapter3,Iusedthevalueofeveryttedparameterthatwasestimatedusingapointestimatebasedonthemaximumlikelihoodsolution,i.e.,thelargestvalueofgoodness-of-ts(minimumof2).AndIusedre-samplingmethods,i.e.bootstrap,tocalculatetheuncertaintyofthettedorbitalparametersbasedontheRVobservation.Thisisatraditionalfrequentistwayofstatisticaldataanalysis.FortheproblemofexoplanetdetectionwithRVdata,apaperreportinganexoplanetdiscoverytypicallyfollowsthefollowingprocedure:1)First,aperiodogramanalysisdeterminesthemostpromisingorbitalperiodanddeterminesafalsealarmprobability(FAP)[ 24 ].Thisperiodisusedtoobtainthebest-torbitalsolutionsvialocalminimizationinstep2.TheFAPisthetypeIerrorgiventheNullHypothesis(noplanetpresents)andAlternativeHypothesis(oneplanetpresentwithcircularorbit).AsmallerFAPmeansgreatersignicanceofrejectingtheNullHypothesis. 2)Usealocalminimizationalgorithm(e.g.L-M)toobtainthebesttKeplerianorbitalmodelbasedoncalculatingthegoodness-of-tstatistic.2=Xk(v~(tk,jk))]TJ /F3 11.955 Tf 11.95 0 Td[(vobs(tk,jk))2=2k. (4) 78

PAGE 79

3)Userettingtechniques,e.g.bootstrap,toobtaintheuncertaintyofttedparameters. 4)Usethescrambledvelocityresidualapproach[ 34 ]tomeasuretheFAPallowingfornon-circularKeplerianmodelbytestingthehypothesisthatthenoiseinthedatamighthaveresultedin2aslongasobservedandledtoafalsedetection. 5)Usetheresidualfromthetofstep2andrepeatfromstep1tocheckifthereisadditionalplanetdetectedbasedonthedata. However,thereareseveraldisadvantagesoffrequentistmethodssummarizedabove.First,thecalculationoftheuncertaintyofthettedparametersandtheprobabilityoffalsedetectionarebasedonthesimulateddataratherthantherealdata.Thislacksarigorousmathematicalbasis.Second,theuncertaintyofthettedparametersarebasedontheuncertaintyofthere-samplingofpointestimatorsandthiscanunder-estimatetherangeoftheuncertaintyofthettedparameters.Third,themodelselectionusingtheresultoffrequentistmethodsisalsobasedonacomparisonbetweentwomaximumlikelihoodvalues,whichmightgivebiasedinformationespeciallywhenanadditionalplanetismarginallydetected. Onthecontrary,Bayesianinferenceprovidesrigorousmathematicalbasistoanalyzethestatisticsofdifferentmodelsusingposteriorprobabilityonlyonbasedontherealobservationdata,andtheresultcanbeusedtomakechoiceofmodelsthemoststronglysupportedbythedata.BayesianmodelsselectionisnotnewinmanydisciplinesoreveninAstronomy,butforexoplanetdetection,duetothepotentialhighdimensionproblemsandthepropertiesofthedata,itrequiresparticulartreatment. 4.2.1Outline Inthischapter,IpresentamethodtoconductBayesianmodelselection.IfocusondevelopmentofanmethodcalculatingtheMarginalLikelihoodofthePosterior,orBayesianEvidence(seeSection 4.3 fordetails).Thetopicsareorganizedasfollows.InSection 4.3 ,IintroducetheBayesianmodelselectionframework,andtheconstructionof 79

PAGE 80

MarkovChainMonteCarlo(MCMC)tosamplefromtheposteriorprobabilitydistributionforRVmodelparameters[ 14 ].InSection4.3,IshowthedevelopmentofanestimatornamedFractionEstimatorforEvidence.AndinSection4.4,ItesttheFractionEstimatorusingrealRVdata. 4.3TheBayesianModelSelectionFramework RecallthederivationofBayes'theorem:giventwoeventsAandB,thejointprobabilityofAandBequalstotheprobabilityofAgiventhatBoccurstimestheprobabilityofBoccurringonitsown,ortotheprobabilityofBgiventhatAoccurstimestheprobabilityofAoccurringonitsown: p(AjB)p(B)=p(A,B)=p(B,A)=p(BjA)p(A), orifweconsiderbothAandBareconditionaloninformationI, p(AjB,I)p(BjI)=p(A,BjI)=p(B,AjI)=p(BjA,I)p(AjI),(4) ThenwehaveBayes'Theorem:p(BjA,I)=p(AjB,I)p(BjII) p(AjI) (4) Andthedenominator,theprobabilityofA,isthesumofthejointprobabilityofAandBoverallofthepossibleoutcomesofB,whichisalsocalledthemarginalprobabilityofA, p(AjI)=XBp(A,BjI).(4) IntheBayesianinterpretation,Bayes'theoremexpresseshowasubjectivedegreeofbeliefshouldrationallychangetoaccountforevidence.IfwereplaceAwith~D,thedata 80

PAGE 81

weobservedinanexperiment,IwithM,themodelweusetoexplainthedata,andBwith~f1,2...g,theparametersinmodelM,( 4 )canbewrittenas,p(~j~D,M)=p(~Dj~,M)p(~jM) p(~DjM), (4) wherep(~j~D,M)isthePosteriorProbabilityoftheparameter~ofmodelMgiventheobservation~D,p(~Dj~,M)istheLikelihoodofthedata~Dgiventheparameterof~inthemodelofM,p(~jM)isthePriorProbabilityoftheparameterinmodelM,whichexpressesourknowledgeofthemodelparameterbeforewehavetheobservationof~D,andp(~DjM)istheMarginalLikelihoodorIntegratedLikelihoodoftheobservation~DortheModelEvidence,p(~DjM)=Zp(~Dj~,M)p(~jM)d~ (4) whereRdimdenotestheparameterspacewithdimensionofdim.Marginallikelihoodisanameintermsofhowitisconstructed,evidenceisanameintermsofitsfunction.Hereafter,IwillusethetermEvidence. ConductingBayesianmodelselectioninvolveswithcalculatingtheBayesianfactor,whichis,Bi,j=p(Mij~D) p(Mjj~D) (4) wherep(Mij~D)andp(Mjj~D)aretheposteriorprobabilityofmodelMiandMjgiventheobservationdata~D,respectively.AndBi,j=p(Mi) p(Mj)p(~DjMi) p(~DjMj). (4) Becauseweusuallylacktheknowledgeofthepriorprobabilityofparticularmodel,itisreasonabletoregardsthepriorprobabilityratiofactorp(Mi) p(Mj)asunity.Soweconsiderthe 81

PAGE 82

evidenceratioEi,j=p(~DjMi) p(~DjMj), (4)Ei,j=Rp(~Dj~,Mi)p(~jMi)d~) Rp(~Dj~,Mj)p(~jMj)d~). (4) AvalueofEi,j>1meansthatMiismorestronglysupportedbythedataunderconsiderationthanMj.Notethatclassicalhypothesistestinggivesonehypothesis(ormodel)preferredstatus(the'nullhypothesis'),andonlyconsidersevidenceagainstit.Jeffreys[ 25 ]gaveascalesubjectiveforinterpretationofBi,j:ifBi,jislessthan1:1,thenMjissupported;ifBi,jisbetween3:1to100:1,thenthestrengthofevidenceforMiisstrong;ifBi,jisgreaterthan100:1,thentheevidenceforMiisdecisive. Inastronomy,moreextremevaluesaretypicallyrequired.Theproblemhasbeenreducedtocalculationofevidenceofeachmodel,amulti-dimensionintegralp(~DjM)(Equation 4 ).However,inmanycases,directcalculationofsuchanintegralisnoteasy. 4.3.1BayesianModelSelectionFrameworkinRVExoplanetDetection NowwerestatetheprobleminspecicexoplanetdetectionwithRVobservation,followingtheschemein[ 14 ].WeconsiderasetofRVdata~vobservedattimes~twithobservationuncertainty~.BecauseintheRVmeasurementarebasedontheaverageofalargenumberofspectrallines,theRVuncertaintythereforecanberegardedasgaussian,andthelikelihoodofobservingsuchdatagivencertainmodelMwithparameters~is,L(~jM)p(~vj~,M)=Yk1 p 2(2k)exp)]TJ /F6 11.955 Tf 10.49 8.36 Td[((v~(tk,jk))]TJ /F3 11.955 Tf 11.96 0 Td[(vk)2 2(2k), (4) whereklabelsthekthobservation,vkistheRVmeasurementofthestarrelativetojkthvelocityreference,v~(tk,jk)isthevelocitygeneratedwithgivenmodelM(Keplerian 82

PAGE 83

modelwithpplanets)attimetkonjthobservatory,and2kistheuncertaintyofvk:2k=2k,obs+2+, (4) where2k,obsistheuncertaintyofkthobservationmeasuredfromthespectraonlyand2+=2jitter+2other, (4) where2jitteristhestarjitternoiseand2otherrepresentsanyothernoisesources,e.g.,noisefromunseenadditionalplanetaroundthestarobserved. BasedonBayesiantheoremframework,theposteriorprobabilityofmodelparameter~p(~j~v,M)=p(~jM)p(~vj~,M) Rp(~jM)p(~vj~,M)d~=p(~jM)L(~) p(~vjM). (4) Thepriorprobabilityp(~jM)ofparametersgivenmodelMcanbeobtainedfromourpreviousknowledgeabouttheexoplanetarysystem.Theposterioristheprobabilitydistributionof~ofmodelMgiventheRVobservation~v.Andsinceforeveryvalueof~,the 4 holds,andtheposteriorsatisesZp(~j~v,M)d~=1, (4) theposteriorprobabilityisthenormalizedfunctionofp(~jM)p(~vj~,M)andthenormalizationfactorisjusttheevidence p(~vjM)=Zp(~jM)p(~vjM,~)d~.(4) 83

PAGE 84

Inpractice,thecalculationoftheevidenceisparticularlydifcult.RecallforexoplanetRVdata,basedonKeplerianmodelwithpplanets~(t)=XpKp[cos(!p+p)+epcos(!p)]+Cj (4) wherepdenotespthplanet,Kpisthesemiamplitudeofthereexvelocitycausedbypthplanet,isthetrueanomaly,eistheeccentricity,!istheargumentofperiastron,andCjisthereferencevelocityconstantatjthobservatory.InChapter1,weintroducedthatthefreeparametertodescribeaKeplerianorbitofplanetareperiodP,e,K,!,phaseormeananomalyatcertainepochM0,inclinationiandthelongitudeofascendingnode.Thereare7moreindependentparametersforeachadditionalplanet.RVcannotmeasureiand,whichreducesthenumberofparameterforeachplanetto5,butsometimeswewouldliketoobtaintheposterioroftheCj,+,oralineartrendconstantdvdt,etc.ThereforethedimensionofparameterspaceofKeplerianplanetsystemisapproximatelydim7NpforageneralmodelwithNpplanetsanddim5NpforaRVnon-interactingmodelwithNpplanets.Thehighdimensionoftheparameterspacemakesthecalculationoftheevidenceextremelydifcult.Whenweareonlyinterestedinthedistributionanduncertaintyoftheparametersincertainmodelgiventheobservation,thecalculationoftheevidencecanbeunnecessary.MarkovChainMonteCarlo(MCMC)hasbeendevelopedinrecentyearstosamplefromtheposteriorprobabilityfunction,whichwillbebrieyintroducedinthefollows. ThechallengesforBayesianmethodsliesincalculatingsuchposteriordistributionswhichincludesanintegralofhighdimensionsofvariables.Inagivenmodelwithxednumberofplanets,theMarkovChainMonteCarlo(MCMC)algorithms[ 14 ]canefcientlyproducesampleswithafrequencyproportionaltotheposteriorprobabilityoftheparametervector.Thenthemarginalposteriordistributionofoneparametercanbecalculatedbymarginalizingovertheotherparameters.However,ifwewouldliketocomparedifferentmodelstodeterminehowmanyplanetareintheobservedsystem 84

PAGE 85

basedongivenobservationdata,thenweneedaproportionalmeasurementofthemarginallikelihoodindifferentmodels.ThischapterdescribestheconstructionandtestingofsuchalgorithmstoestimatethemarginallikelihoodusedinBayesianmodelselectionformultipleplanetRVdetection. 4.3.2BackgroundofMarkovChainMonteCarlo,Metropolis-HastingsAlgorithmandGibbsSampling MarkovchainMonteCarlo(MCMC)methodsaredesignedtofacilitatesamplingfromcertaintargetdensitybygeneratingaMarkovchainwhichhasthetargetdensityasitsinvariantdistribution.Ford[ 13 ]hasdevelopedanalgorithmtogenerateaMarkovchainofstates(i.e.,setofparametervalues)thataresampledfromtheposteriorprobabilityfunction 4 ,usingMetropolis-Hastingsalgorithm(M-H)andGibbssampling. M-H,developedby[ 36 ]andgeneralizedby[ 23 ],isaMCMCmethodtoobtainasequenceofrandomsamples(Markovchain)fromaprobabilitydistributionwhosedirectsamplingisdifcult,whichisjustinthecaseoftheposteriorprobabilityfunctionofmodelparametersgivenRVdata.AGibbssamplerspeciesthatonlyassubsetof~ischangedateachstepoftheMarkovchain,thereforeenablethesamplinginamulti-dimensionalparameterspace.Thecomputationefciencyhasbeensignicantlyimprovedby[ 14 ].Onecanrefer[ 13 ]and[ 14 ]forthedetailsabouttheconstructionandtechnicalconsiderationoftheMCMCforRVdata,andthereferencethereinforthetheoryabouttheMCMCwithM-HandGibbssampler.Ijustbrieysummarizebelow. AtcertainstateoftheMarkovchainwithvalueof(~i),atrialstate(~0)isproposedbydrawingfromacandidatetransitionprobabilitydistributionfunction(CTPDF)q(~0j~i,M),whichisaconditionalprobabilityfunctionbasedoncurrentstate(~i)andthemodelM. 85

PAGE 86

TheacceptanceprobabilitycanbeconstructedbytheM-Halgorithm[ 36 ],(~0j~i,M)=minfp(~0j~D,M)q(~ij~0,M) p(~ij~D,M)q(~0j~i,M),1g=minfp(~0jM)p(~Dj~0,M)q(~ij~0,M) p(~i0jM)p(~Dj~i,M)q(~0j~i,M),1g (4) Wherep(~jM)istheposteriorofmodelparameter,p(~jM)isthepriorand(~Dj~,M)isthelikelihood.Thetrialstate(~0)wasacceptedwithaprobability(~0j~i,M).Ifitisaccepted,thenthenextstateoftheMarkovchain(~i+1)=(~0),otherwise,(~i+1)=(~i).TheGibbssamplerupdateonlyonedimensionoftheparametereachstep.NoteinEquation 4 theevidenceinthenumeratoranddenominatorcanceled,thustheevidenceneednottobecalculatedtosamplefromtheposterior. AfteralongtermrunofMCMC,calledtheburninstage,itcanbeshownthatthechainwillconvergetothetargetprobabilityp(~jM).ThechoiceofCTPDFdeterminestheefciencyoftheMCMC.Thereareseveralalternativesforitsuitabledifferentscenarios,whicharediscussedin[ 14 ].InasequenceofMCMCsamples,eachstatecontainsasetofvaluesoftheparameterthatrepresentonesampleformtheposteriordistribution,soitiseasytoobtainthemarginalizedposteriorprobabilityforeachmodelparameter.Byplottingahistogramoveronedimension,themarginalposteriorcanbeshown,andbyamulti-dimensionalhistogramorcontourplot,wecanvisualizethepartiallymarginalposteriorofcertaindimensionsoftheparametervector.Variousstatisticscanbecomputedfromthesesamples.InSection 4.4 ,Iusethealgorithmsfrom[ 14 ]togenerateMCMCsamplesfortestingalgorithmstoevaluatetheevidence. 4.4CalculatingEvidenceBasedonMCMCSamples ThischapterconsidershowtocalculatetheevidencebasedonanMCMCsamplesowecanperformBayesianmodelselection.Severalpreviousalgorithmshavebeenappliedtothisproblem.InFord&Gregory[ 15 ]severalestimatorswereevaluated,usingrealRVdatafromasingleplanetHD88133,includingtheimportancesampling 86

PAGE 87

estimators,theratioestimator,andparalleltempering(seethepaperfordetaileddenitionoftheestimators),buttheseneedmoretest,especiallyformulti-planetsystem.Bovyetal.[ 2 ]developedExtremeDe-convolution,whichuseamixtureofGaussiantosimulatetheposterior;Loredoetal.[ 30 ]useimportancesampling;Weinberg[ 42 ]usedtheweightedharmonicmeanestimator;andtheMultinestSamplingby[ 9 ]useLebesgueintegralsamplinginpriorspacetocalculatetheevidence. IfwealreadyhaveasequenceofMCMCsamplesthathavefullyexploredtheparametersspaceandconvergedintotheposteriorprobabilityofthemodelparameters,thenthecalculationoftheBayesianEvidencehasbeenreducedtoaproblemofcalculatingthenormalizationfactorofafunctiongivenagroupofsamplesdrawnfromthisfunction.Inthefollowing,IproposeasimplealgorithmtocalculatetheBayesianevidence,andshowthatthisalgorithmhascontrollableuncertainty. 4.4.1TheFractionEstimator Aswedescribeabove,theproblemcanbefurthersimpliedandstatedasfollows, Considerthefunctionsg(x)andf(x),whereisproportionaltog(x)byanormalizationfactorE,E=g(x) f(x) (4) whereEZg(x)dx, (4) andZf(x)dx=1 EZg(x)dx=1 (4) Letfxig,i=1,2...NbeasequencewithlengthofNsampledfromf(x)org(x),i.e.,xif(x).UseVN=todenotethefulldomainthatcontainsallNsamplepoints.Ifa 87

PAGE 88

sub-domainVnVNispickedcontainingn
PAGE 89

Accordingtotheerrorpropagationtheory,thevarianceoftheFractionEstimatorintroducedin( 4 )is,Var(E) E2Var(Gn(x)) Gn(x)2+Var(Fn(x)) Fn(x)2)]TJ /F6 11.955 Tf 11.96 0 Td[(2G,F(Gn(x))(Fn(x)) Gn(x)Fn(x), (4) whereisthestandarddeviationandisthecorrelationcoefcient.SincewecalculateGn(x)andFn(x)independently,G,F=0.Therefore,Var(E) E2Var(Gn(x)) Gn(x)2+Var(Fn(x)) Fn(x)2, (4) 4.4.2.1VarianceofFn(x) ThecalculationofFn(x)isstraightforward.IfwespecifyaVn,thenweonlyneedtocountthenumberofsamplesxithatfallinVn,andcomputetheration N.Sincexifollowthedistributionf(x),assumingthesamplesdrawnfromthetargetdistributionareindependentfromeachother,asinglerandomsamplefromVNisdrawnfromVnattheprobabilitypVn=Zx2Vnf(x)dx. (4) Ifeachsampleisindependent,thenthetotalnumbern,countedinVnoutofN,followstheBinomialdistribution,i.e.,nB(N,pVn).Sotheestimationofnis,E(n)=NpVn, (4) anditsvariancewouldbe,Var(n)=NpVn(1)]TJ /F3 11.955 Tf 11.96 0 Td[(pVn). (4) Son=NpVnp NpVn(1)]TJ /F3 11.955 Tf 11.95 0 Td[(pVn) (4) 89

PAGE 90

SinceFn(x)=n Nwehave,E(Fn(x))=NpVn N=pVn, (4) andVar(Fn(x))=Var(n) N2=pVn(1)]TJ /F3 11.955 Tf 11.96 0 Td[(pVn) N. (4) SoVar(Fn(x)) Fn(x)2=1 N(1 pVn)]TJ /F6 11.955 Tf 11.95 0 Td[(1) (4) Equation 4 suggeststhattherearetwowaystodecreasetherelativeerroroftheestimationforthefractionn=N,orthatofthedenominatorin( 4 ):(i)increasingthenumberofthetotalsamplepointsN;(ii)increasingpVn,theprobabilitythatarandompointsfallingintoVn,whichmeansselectingalargeVnorselectingasub-domainwithlargef(x). 4.4.2.2VarianceofGn(x) Theformererroritemintherightsideof 4 couldbedeterminedbythealgorithmscalculatingGn(x),i.e.,theintegralofg(x)overVn.Sincewearefacingthepotentialhighdimensionalintegralproblem,pseudo-randomnumbernumericalintegrationtechniquesuchastheMonteCarloIntegration(MCI)canbeused.Bygeneratingrandompointsovertheindependentvariabledomain,MCIestimatestheintegralofafunctionoverthisdomain.Inourcase,whenNMCI,thenumberofrandompointsdrawnbyMCIoverVnislargeenough,wehave[ 37 ],Zg(x)dVVnhg(x)iVnVns hg(x)2iVn)-221(hg(x)i2Vn NMCI, (4) whereVnalsomeansthevolumeofthissub-domain,andhg(x)iVn1 NMCINMCIXi=1g(xi), (4) 90

PAGE 91

andhg(x)2iVn1 NMCINMCIXi=1g(xi)2, (4) thatis,Var(g(x))Vn=hg(x)2iVn)-222(hg(x)i2Vn1 NMCINMCIXi=1g(xi)2)]TJ /F6 11.955 Tf 11.96 0 Td[((1 NMCINMCIXi=1g(xi))2, (4) whichmeans,E(Gn(x))=Vnhg(x)iVn, (4) Var(Gn(x))=V2nVar(g(x))Vn NMCI, (4) andGn(x)=ZVng(x)dxVnhg(x)iVnVns Var(g(x))Vn NMCI. (4) SotherelativevariancecausedbyGn(x)becomes,Var(Gn(x)) Gn(x)2=Var(g(x))Vn hg(x)Vni2NMCI. (4) Equation 4 infersthattherelativeerrorfromGn(x)canbereducedby1)choosingasub-domainVnwithlowervariance,2)choosingVnwithhigheraveragevalueofg(x),and3)increasingNMCIsincethevarianceofGn(x)decreasesinthescaleof1=NMCI. 91

PAGE 92

4.4.2.3CombineGn(x)andFn(x) Baseonthediscussionabove,therelativevarianceoftheFractionEstimatorasshowninEquation 4 Var(E) E2=Var(g(x))Vn hg(x)i2VnNMCI+1 N(1 pVn)]TJ /F6 11.955 Tf 11.95 0 Td[(1), (4) OntherightsideofEquation 4 ,bothitemsarerelatedwiththechoiceofVn.Inordertominimizedthevalueofthewholeequation,therstitemrequireslargevalueofhg(x)iVnandsmallVar(g(x))Vn,andtheseconditemrequiresaslargevalueofpVnaspossible.AlargevalueofpVncanberealizedbychoosingthedomainwhereg(x)orf(x)islarge,orincreasethesizeofVn.ThemaximumvalueofpVnis1,however,itrequiresusingthewholedomainVN,butthelagerVnwillincreasethevarianceduetotherstitemontherightsideofEquation 4 .IfVnistoobig,itcausesVar(g(x))Vntobelarger,soaverylargeNMCIwouldbeneededtosatisfytheprecisionrequirement.Basedontheseconsideration,theoptimizedchoiceofVnistheasub-domainaroundtheglobalmodeoff(x).BesidesthelocationofVn,thesizeofVnneedstobecarefullychosen,too.ThiswillbediscussedintheSection 4.4.3 Furthermore,bothitemsdecreasewith1=NMCIor1=N,soincreasingthesamplesizeNandMCIdrawsNMCIreducestherelativevariance.ConsiderthesamplesfortheevidencecalculationarebasedonMCMC(replaceNwithNMCMCinEquation 4 ),whereeachstepisnotindependent,soincreasingNMCIiscomputationallycheaperthanincreasingNMCMC,becauseeachstepofthecalculationofg(x)inMCIrequiresonlyonemodelevaluation,whileeacheffectivelyindependentstateofMCMCrequiresseveralmodelevaluation. 4.4.3FractionEstimatorwithMCMCSamples IntheSection 4.4.1 ,weconstructedtheFractionEstimatorunderasimpliedsituation.WemadeanassumptionthatxiareIndependentandidenticallydistributedandfollowsg(x)oritsnormalizedfunctionf(x).Thisassumption,however,isan 92

PAGE 93

approximation,sinceinarandomwalkMarkovchain,eachstateisproposedbythepreviousstatewiththeCTPDF,resultinginsamplesthatareauto-correlated.Var(Fn(x))mayunderestimatethevarianceoftheestimator.Ifthechain'slengthislongenough,andweuseeverymthstatefromthenalMCMCsamplefortheintegrationoff(x)andmislargerthantheautocorrelationlengthoftheMarkovChain,thentheoverallautocorrelationbetweenstatescouldbeneglectedandtheaboveapproximationforVar(Fn(x))wouldbesafe. 4.4.4Multi-dimensionalModelSelection GivenasequenceofMCMCsamplesthathavefullyexploredthedimdimensionfunctionf(~),weuse~=(1,2,..d,...dim)toreplacexanddtodenotewhichparameter.TocalculatepVninEquation 4 ,itisconvenienttochooseahypercubeVn,whichineachdimension,d,ldd,u,whered,landd,uarethelowerandupperboundofd,andVn=Ydim(d,u)]TJ /F7 11.955 Tf 11.95 0 Td[(d,l)=Ydimd (4) SinceusingtheMCMCsamplestoobtainthemarginalposteriorineachdimensionisaseasyastoplotthehistogramofthedthdimensionofthesamples,wecanuselocatedtocoverthemodeindthdimension.Andwecandene,pndp(nd)=p(dd)=nd N (4) 93

PAGE 94

wherenddenotetheeventthattherearendoutofthetotalNsamplessatisfyingtheconditiondd.SopVnisthejointprobabilityofnd,pVn=p(n1,n2,...ndim)=p(n1,n2,...ndim)]TJ /F5 7.97 Tf 6.58 0 Td[(1jndim)p(ndim)...=p(n1jn2,...ndim)p(n2jn3,...ndim)...p(ndim)]TJ /F5 7.97 Tf 6.59 0 Td[(1jndim)p(ndim) (4) Inthesimplestcase,thejointprobabilityofp(n1,n2)=p(n1jn2)p(n2):ifn1andn2areindependent,thenp(n1,n2)=p(n1)p(n2).Ifn1andn2arecorrelated,thenp(n1jn2)p(n1).Andinextremecase,p(n1jn2)1,when1and2arehighlycorrelated.So 4 becomes,pVnp(n1)p(n2)...p(ndim), (4) wheretheequalityholdsifandonlyifdareindependentfromeachother.SowehavefoundthelowerboundofpVn,thustheupperboundof,1 N(1 pVn)]TJ /F6 11.955 Tf 11.96 0 Td[(1)1 N(1 Qdimp(nd))]TJ /F6 11.955 Tf 11.95 0 Td[(1) (4) 4.4.5TestwithMCMCSamplesfromMulti-variantGaussian DuringthegenerationofMCMC,theidealCTPDFisthetargetprobabilityitself,whichisimpossibleinpracticesincethetargetprobabilityisunknown.However,inordertotesttheFractionRatioEstimator,wecanconstructMCMCusingasimpleanalyticprobabilitydistribution,e.g.,Multi-variantGaussian,asthetargetprobabilityandtheCTPDF,andwatchtheperformanceofthisestimatorwiththelengthoftheMarkovchain,thedimensionandthechoiceofthefractionvolume. 4.4.5.1GeneratingMCMCsamplesfrommulti-variantGaussian InpracticefortheMCMCsamplesgeneratedbasedontherealdata,thetargetdistribution,whichislikelihoodtimespriori,isasmoothfunctionovertheparameter 94

PAGE 95

space,soitisreasonabletoapproximatethefunctionaroundamodewithamulti-variantGaussian.Inordertosimplifythecalculation,Ichooseamulti-variantGaussiandistributionwith0meanandidentitycovariancematrix,i.e.,multi-variantnormaldistributionasthetargetdistributionandsamplefromitusingMCMCalgorithms.ThesamedistributionisalsousedastheCTPDFforefcientsampling.p(~)N(,) (4) whereisaconstantdimdimensionvectorwithvalueof(0,0,0,...0),andisadimdimdimensionidentitymatrix. UsingM-HalgorithmsandGibbssampling,fourMCMCchainsaregenerated.Eachofthemhas7,14,21or28dimensionsintheparameterspace,whichisequaltothenumberoforbitalparametersforgeneralmodelwith1,2,3or4planets.Theintegralofthemulti-variantnormaldistributionscanbeanalyticallyevaluated,soistheevidence.Therefore,wecanuseittocomparewiththevaluecalculatedusingtheFractionEstimator. Figure 4-1 showsafewexamplesofthemarginalposteriorprobabilitydistributionsthemulti-variantnormaldistributionparameters.TheyarethehistogramplotsoftheMCMCsamplesoftherandomlyselectedparameters.EachpanelshowsaGaussiandistributionwithzeromeanandunitstandarddeviation,whichimpliestheMCMChasconvergedtothetargetdistributions. 4.4.5.2EvidencecalculationwithMCMCsamplesfrommulti-variantGaussian TheevidencesoffourchainsarecalculatedusingtheFractionEstimatorasafunctionofthesizeofthesub-domainVn,asshowninFigure 4-2 4-3 4-4 and 4-5 .Thesub-domainforFractionEstimatorischosenasahypercubewithaidenticalrangeineverydimensionanditisasymmetricrangearoundthemodewhichisalsothemeanofthemarginaldistribution,i.e.,0.xgivesthehalfwidthoftherangeineachdimensioninunitofthestandarddeviationwhileygivethevalueoftheevidence.In 95

PAGE 96

Figure4-1. ExamplesofmarginalposteriorprobabilitydistributionsofrandomlyselecteddimensionsfromasequenceofMCMCsamples.TheMCMChasconvergedintothetargetprobabilitydistribution,withisamulti-variantGaussianwithzeromeanvectorandidentityvariancematrix. theplots,thetruevalueoftheevidenceis1.ThesolidcurveistheevidenceestimatedbytheratiooftheanalyticvalueoftheintegrationovertheMCMCsamplefractioninthesub-domain.ThedashcurveistheevidenceestimatedbytheratiooftheMonteCarlointegrationovertheMCMCsamplefractioninthesub-domain,i.e.,bytheFractionEstimator.Forclearcomparisons,allchainsaretruncatedintothesamelength(NMCMC=180000)despitethedifferentdimensions.AndthedrawsofMonteCarlointegrationareproportionaltothedimensionsoftheparameterspace,i.e.,NMCI=70000,140000,210000,280000,respectively.Thedeviationofthesolidcurvefrom1impliestheestimationerrorcausedbythecalculationofthefractionoftheMCMC 96

PAGE 97

samplesinthesub-domainVn,Fn,andthedeviationofthedashcurvefromthesolidcurveimpliestheestimationerrorfromtheMonteCarlointegration. Whenthesizesub-domainVnistoosmall0.1aroundthemode,oursearchfailstondanyMCMCsamplesinVn.Asthedimensionincreases,withxedlengthoftheMCMCsequence,itrequireslargerrangeineachdimensiontondsamplesinVn.Itisareasonableresultsincethesizetheofparameterspaceincreasegenerallywiththeorderofdim.Ifwechoosea2rangearoundthemodeineachdimensionandassume3canoverallofthesamples,itgivesVn VN=2 3dim,whichmeansthefractionofthesub-volumedecreaseintheorderofdim. InordertousesufcientnumberofsamplestocalculateFnaccurately,onewayistoincreasethesizeofVn,astheFigure 4-2 4-3 4-4 and 4-5 showthatthesolidcurveconvergeto1asthesizeincreases.However,itconvergesslowerinhigherdimensionthanitdoesinlowerdimensionforxed.SoanotherwayistoincreasethelengthoftheMCMCsequenceasthenumberofthemodelparametersincreases. AlthoughaccordingtothecalculationofGninEquation 4 ,thereisnoimplicittermrelatedwiththedimension,thevarianceofthetargetdistributiong(x)doesincreasewithincreasingdimensionofthemodelparameterspace.ThedashcurveinFigure 4-2 4-3 4-4 and 4-5 indicatestheaccuracyoftheMonteCarlointegration:ifitdivergesfromthesolidcurve,thenitimpliesthevalueofGncalculatedbyMCIisdifferentfromanalyticvalue.For7-dimension,theMCIworksaswellastheanalyticalevaluationofGnwith70000drawsover0.3)]TJ /F6 11.955 Tf 12.3 0 Td[(3.Asthedimensionincreases,theMCIstartstogiveincorrectestimationofGnforx>2.Notethenumberofthedrawsisalreadyincreasedproportionallywiththedimensions. 4.4.5.3TherelationbetweenNMCMCandNMCI AlthoughwecouldincreasetheaccuracyoftheevaluationofGnusingMCIbyincreasingNMCI,alargevalueofNMCIisalsocomputationallyexpensive.However,in 97

PAGE 98

Figure4-2. EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma7-variantNormaldistribution.x-axis:thehalfwidthoftherangeineachdimensioninunitofthestandarddeviation;therangeineachdimensionissymmetricaroundthemode,whichis0.y-axisgivesthevalueoftheevidence.ThetotalnumberofMCMCsamples(NMCMC)andthedrawsoftheMonteCarlointegration(NMCI)aredenotedinthelegend.Thetruevalueoftheevidenceis1.Thesolidcurve:theevidenceestimatedbytheratiooftheanalyticvalueoftheintegrationovertheMCMCsamplefractioninthesub-domain.Thedivergenceofthesolidcurvefrom1indicateserroroftheestimationcausedbythefractionestimationFn.ThedashcurveistheevidenceestimatedbytheratiooftheMonteCarlointegrationovertheMCMCsamplefractioninthesub-domain,i.e.,bytheFractionEstimator.ThedashcurveisoverlappedwiththesolidcurveiftheMonteCarlointegrationisaccurateenough.AndthedivergenceofthedashcurvefromthesolidcurveindicatedtheestimationerroroftheplainMonteCarlointegration.NotetheplottingbeginswhenthesizeofVnlargeenoughtondanyMCMCsampleswithinit. 98

PAGE 99

Figure4-3. EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma14-variantNormaldistribution.x-axis:thehalfwidthoftherangeineachdimensioninunitofthestandarddeviation;therangeineachdimensionissymmetricaroundthemode(0).y-axisgivesthevalueoftheevidence.ThetotalnumberofMCMCsamples(NMCMC)andthedrawsoftheMonteCarlointegration(NMCI)aredenotedinthelegend.Thetruevalueoftheevidenceis1.Thesolidcurve:theevidenceestimatedbytheratiooftheanalyticvalueoftheintegrationovertheMCMCsamplefractioninthesub-domain.Thedivergenceofthesolidcurvefrom1indicateserroroftheestimationcausedbythefractionestimationFn.ThedashcurveistheevidenceestimatedbytheratiooftheMonteCarlointegrationovertheMCMCsamplefractioninthesub-domain,i.e.,bytheFractionEstimator.ThedashcurveisoverlappedwiththesolidcurveiftheMonteCarlointegrationisaccurateenough.AndthedivergenceofthedashcurvefromthesolidcurveindicatedtheestimationerroroftheplainMonteCarlointegration.NotetheplottingbeginswhenthesizeofVnlargeenoughtondanyMCMCsampleswithinit. 99

PAGE 100

Figure4-4. EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma21-variantNormaldistribution.x-axis:thehalfwidthoftherangeineachdimensioninunitofthestandarddeviation;therangeineachdimensionissymmetricaroundthemode,whichis0.y-axisgivesthevalueoftheevidence.ThetotalnumberofMCMCsamples(NMCMC)andthedrawsoftheMonteCarlointegration(NMCI)aredenotedinthelegend.Thetruevalueoftheevidenceis1.Thesolidcurve:theevidenceestimatedbytheratiooftheanalyticvalueoftheintegrationovertheMCMCsamplefractioninthesub-domain.Thedivergenceofthesolidcurvefrom1indicateserroroftheestimationcausedbythefractionestimationFn.ThedashcurveistheevidenceestimatedbytheratiooftheMonteCarlointegrationovertheMCMCsamplefractioninthesub-domain,i.e.,bytheFractionEstimator.ThedashcurveisoverlappedwiththesolidcurveiftheMonteCarlointegrationisaccurateenough.AndthedivergenceofthedashcurvefromthesolidcurveindicatedtheestimationerroroftheplainMonteCarlointegration.NotetheplottingbeginswhenthesizeofVnlargeenoughtondanyMCMCsampleswithinit. 100

PAGE 101

Figure4-5. EvidenceestimationusingtheFractionEstimatorbasedontheMCMCsamplesfroma28-variantNormaldistribution.x-axis:thehalfwidthoftherangeineachdimensioninunitofthestandarddeviation;therangeineachdimensionissymmetricaroundthemode,whichis0.y-axisgivesthevalueoftheevidence.ThetotalnumberofMCMCsamples(NMCMC)andthedrawsoftheMonteCarlointegration(NMCI)aredenotedinthelegend.Thetruevalueoftheevidenceis1.Thesolidcurve:theevidenceestimatedbytheratiooftheanalyticvalueoftheintegrationovertheMCMCsamplefractioninthesub-domain.Thedivergenceofthesolidcurvefrom1indicateserroroftheestimationcausedbythefractionestimationFn.ThedashcurveistheevidenceestimatedbytheratiooftheMonteCarlointegrationovertheMCMCsamplefractioninthesub-domain,i.e.,bytheFractionEstimator.ThedashcurveisoverlappedwiththesolidcurveiftheMonteCarlointegrationisaccurateenough.AndthedivergenceofthedashcurvefromthesolidcurveindicatedtheestimationerroroftheplainMonteCarlointegration.NotetheplottingbeginswhenthesizeofVnlargeenoughtondanyMCMCsampleswithinit. 101

PAGE 102

termsofthecomputationcomplexity,wecanconstructtherelationbetweenNMCMCandNMCI. DuringthegenerationofMCMCsamples,thescalesoftheCTPDFarechosentoappropriatevaluessothattheMCMCcanefcientlyexplorethefullparameterspaceandsamplefromthetargetprobabilitydistribution.TheacceptancerateofthechainisdeterminedbythescalesoftheCTPDFanditisusuallysettocloseto0.22[ 13 ].Inordertoreducetheautocorrelationamongthestatesinthechain,thestatesaresavedateverynthstates,wherenisproportionaltothenumberofthemodelparameter.Andanumberofstatesrstgeneratedarediscardedastheburninstagesothattherestofthechainisthestationarystatesthathasconvergedintothetargetdistribution.SoaMCMCsequenceofalengthNMCMCisattheexpenseofatleastnNMCMCevaluationofthemodel,whileforMCIwithNMCIdrawsittakesNMCIevaluationsofthemodel.Therefore,itisreasonabletosetNMCInNMCMCifwespendthesamecomputationexpenseonthecalculationofthecalculationofevidence.Ijustsetn=4forthecomputationsimplicity,andusethisrelationhereafterinthecalculationofmodelevidencewiththeFractionEstimator. 4.4.5.4Principlecomponentsanalysis InthetestoftheFractionEstimatorshowninSection 4.4.5.2 ,multi-variantnormaldistributionswithidentityvariancematrixareusedtosimulatethetargetdistribution,i.e.,theparametersinthemodelareuncorrelated.However,inrealproblems,themodelparametersarecorrelated,e.g.,eand!.WhenusingtheFractionEstimatortocalculatethemodelevidence,sincewedonothavetheknowledgeabouttheoverallshapeoftheposteriorprobabilitydistribution,wehavetodecidetherangeineachdimensionormodelparametertocoverthemode,basedonlyonthemarginalposteriorprobabilitydistributionthatcanbeillustratedwithhistogramplotusingtheMCMCsamples.Foraxedrangeweselectforeachdimension,asanalyzedinSection 4.4.3 thenalpVnis 102

PAGE 103

largerforcorrelatedparameterspacethanforuncorrelated,sincemoresamplesareconcentratedinthemodeiftheparametersarecorrelated. However,inthecasetheparametersarecorrelated,thechoiceofVnasahypercubecanleadtoalargevarianceoftheg(x)inVn.Tosimulateaposteriordistributionasanexample,a2-dgaussianisusedwithazeromeanvectorandavariancematrixof=26410.50.51375. (4) Figure 4-6 showsthesamplesdrawnfromthe2-dgaussian(darkpoints).Themarginaldistributionofeachparameterisstillagaussianwithzeromeanandunitvariance.Asub-domainisdeterminedbysettherangeof2=2.0foreachdimension.Theredpointsarethesamplesincludedinthesub-domain.Obviously,theredsquareincludesaregionwithlargelyvaryingsampledensity,implyinglargevarianceoftheposteriordistribution.ThiscanleadtoalargeestimationvariationofGnintheFractionestimator. Principlecomponentsanalysis(PCA)isamethodtouseorthogonaltransformationtoconvertasetofvaluesofcorrelatedvariablesintoasetofvaluesoflinearlyuncorrelatedvariablescalledprincipalcomponents.Figure 4-7 showsthesamesamplesofFigure 4-6 transformedintothePCAspacealongwiththesub-domainVn.Insteadofchoosingarelativesmoothregionoftheposteriordistribution,Vnincludesaregionwithlargelyvaryingposteriordistribution,whichisindicatedbythevaryingsampledensity. IfwedetermineVninPCAspaceusingthemethodappliedfortheorthogonalmulti-variantgaussianandtransformbacktooriginalspace,thenthevarianceoftheposteriordistributioninVncanbereduced.Figure 4-8 andFigure 4-9 showthisprocess.Figure 4-8 isthetransformedparameterspaceusingPCA.TheredregionisVndeterminedinPCAspacebychoosing2rangearoundthemodeineachdimension. 103

PAGE 104

Figure4-6. Choiceofthesub-domaininoriginalparameterspace.Darkpoints:Thesamplesdrawnfromthe2-dcorrelatedgaussian.Redpoints:thesamplesincludedinthesub-domain.Asub-domainisdeterminedbysettherangeof2=2.0foreachdimensionaroundthemodebasedonthemarginaldistributionofeachparameter,whichisstillagaussianwithzeromeanandunitvariance.Thesub-domaindeterminedthiswayincludesregionswithsmallvalueofthedistributionfunctionandlargevarianceofthefunction. Figure 4-9 showsthesamplesinoriginalspace,andtheredregionisVnchoseninthePCAspace.Obviously,thechoiceofVnisalongtheprinciplecomponents.ComparedwiththechoiceofVnshowninFigure 4-6 ,itincludesaregionwithmoresmoothposteriorPDF.InSection 4.5.2 ,IwillusearealRVdatasettoshowtheeffectofusingPCA. 104

PAGE 105

Figure4-7. Choiceofthesub-domaininoriginalparameterspaceshowninthetransformedparameterspaceusingPCAtechniques.Darkpoints:ThesamesamplesfromFigure 4-6 transformedintoPCAspace.Redpoints:VndeterminedinoriginalparameterspacethentransformedintoPCAspace.Thesub-domaindeterminedthiswayincludesregionswithsmallvalueofthedistributionfunctionandlargevarianceofthefunction. 4.5EvidenceCalculationBasedonMCMCSamplesfromRVData 4.5.1HD88133 Ford&Gregory[ 15 ]usedtheRVobservationofHD88133tocalculatetheoneplanetmodelevidenceusingseveraltypesofestimator.InSection 4.5.2 ,IusethesameRVdatasettoexaminetheFractionEstimator.TheRVdatasetfrom[ 10 ]has17observations,asshowninFigure 4-10 .MCMCsamplesaredrawnfromtheposteriorprobabilitydistributionfollowingtheframeworkdescribeinSection 4.3.1 usingthealgorithmintroducedby[ 14 ]. 105

PAGE 106

Figure4-8. Choiceofthesub-domaininthetransformedparameterspaceusingPCAtechniques.Darkpoints:ThesamesamplesfromFigure 4-6 transformedintoPCAspace.Redpoints:VndeterminedinPCAparameterspace.InthePCAspace,thelinearcorrelationsbetweentheoriginalmodelparametersarereduced,thesub-domaindeterminedthiswaycontainsasmuchregionswithlargedistributionvaluesandsmallvarianceaspossible. IusethesamepriorPDFsofmodelparameterasin[ 15 ],whichislistedTable 4-1 .ThepriorPDFsareregardedasconstantsforep,!p,M0,p,Cj,logPp,logKpandlogjitter.Bytakingstepsintheseparameters,theM-HalgorithmavoidscalculatingthepriorPDFssincetheratioofthepriorPDFsbetweenastateofMCMCandthetrialstateproposedbasedonthecurrentstateisequaltoone.However,whenwecalculatetheevidence,wehavetoincludethevalueofthepriorPDFs.HereafterinthisthesisIwillusethesamepriorPDFforotherapplicationsofmodelevidencecalculation. 106

PAGE 107

ModelParameterVariablePriorDistributionMinimumMaximumNote Amplitudeofjitters(s+s0))]TJ /F14 5.978 Tf 5.75 0 Td[(1 log1.+smax s00m/sKmaxs0=1m=sParametersforeachvelocityreferenceVelocityoffsetCj1 Cmax)]TJ /F4 7.97 Tf 6.59 0 Td[(Cmin-Cmax=KmaxCmin=KmaxKmax=2129m=sParametersforeachplanetOrbitalperiodPpP)]TJ /F14 5.978 Tf 5.75 0 Td[(1p log(Pmax=Pmin)Pmin)=1dPmax=103yrVelocitysemi-amplitudeKi(Kp+K0))]TJ /F14 5.978 Tf 5.76 0 Td[(1 log1.+Kmax K0Pmin Pp1=30m/sKmaxPmin Pp1 3K0=1m=sOrbitaleccentricityep101Argumentofperiastron!i1 202OrbitalphaseM0,p1 202 Table4-1. SAMSIExoplanetWorkingGroupReferencePriors.Thevalueofconstantsarealsothesameasin[ 15 ].ReprintedbypermissionfromFord,E.B.andGregory,P.C.2007,StatisticalChallengesinModernAstronomyIV(Table1). 107

PAGE 108

Figure4-9. Choiceofthesub-domaininthetransformedparameterspaceusingPCAtechniquesshownintheoriginalparameterspace.Darkpoints:ThesamesamplesfromFigure 4-6 inoriginalspace.Redpoints:VndeterminedinPCAparameterspacetransformedbackintooriginalspace.InthePCAspace,thelinearcorrelationsbetweentheoriginalmodelparametersarereduced,thesub-domaindeterminedthiswaycontainsasmuchregionswithlargedistributionvaluesandsmallvarianceaspossible. 4.5.2EvidenceCalculationwithFractionEstimator 4.5.2.1Findingmodeofposteriordistribution AfterthegenerationofMCMCforaRVdatawithcertainmodel,whichisoneplanetmodelinthiscase,thenextsteptoapplytheFractionEstimatoristondthemodeintheparameterspaceusingthemarginalposteriorPDFs.Asdiscussedin[ 15 ],theposteriorPDFforonemodelparametersofRVobservationsofHD88133isclosetoamulti-variantGaussian,soitisstraightforwardtodeterminethelocationofthepossible 108

PAGE 109

Figure4-10. RVobservationsofHD88133.BasedonthedatapublishedinFischeretal.[ 10 ]plottedversusthetimeofobservation.ReprintedbypermissionfromFord,E.B.andGregory,P.C.2007,StatisticalChallengesinModernAstronomyIV(Figure1,leftpanel). 109

PAGE 110

globalmode.However,someoftheparamorsoftheRVmodelrequireparticulartransformationinpractice,i.e.,!pandM0,p. AsshowninFigure 4-11 and 4-13 ,both!pandM0,pareperiodicalandaresampledinrangeof[0,360]deg.Thiscanleadtotheapparentdisplacementofthemode.Byaddingorsubtractingoneperiodtothesamplesononesideofthemode,thePDFof!pandM0,pcanbetransformedintotheformasinFigure 4-12 and 4-14 ,whichismoreconvenientformodedetectionusingprogramming. Anothermethodofparametertransformationformodeselectionistousesin!andsinM0orcos!andcosM0insteadof!andM0,ortouseesin!andecos!toreplaceeand!.ThetransformationoftheparameterspacemightleadtoamoresmoothposteriorPDFwithintheregionwebelievewherethemodeis.ButforthisprojectwekeepusingthestraightforwardformofthemodelparametersoutputtedbyMCMC. 4.5.2.2FractionEstimatorwithHD88133RVData Figure 4-15 showstheplotsoftheoneplanetmodelevidencebasedon17RVobservationsusingtheFractionEstimator,withPCAandwithoutPCA.AsthetotalMCMCsampleusedforcalculationincreases,bothPCAestimationandnon-PCAestimationconvergetoaclosedvalue,whilethePCAestimationisslightlybetterthannon-PCA.AconvergedvalueofpVnis0.22forbothestimations,i.e.,22percentofthetotalMCMCsamplesareincludedinthesub-domainVn. Figure 4-16 showstheinternallyestimatederror,orthestandarddeviationoftheestimationrelativetotheestimation,oftheFractionEstimator,withandwithoutPCA.PCAestimationgivesmorethanoneorderbetterrelativestandarddeviationoftheestimationthannon-PCA.WhichmeansPCAisnecessaryinordertoreducethemodelevidenceestimationerror. 4.5.3ComparisonwithOtherEstimators IfwecomparetheFractionEstimatorwithotherestimators,e.g.,theestimatortestedin[ 15 ],whichisplottedinFigure 4-17 ,theFractionEstimatorperformsaswellas 110

PAGE 111

Figure4-11. Transforming!formodeselection.MarginalposteriorPDFof!inrangeof[0,360]deg.!isperiodicparameter.Itisconvenienttotransform!toanewrangewiththemodeinthecenter. thegoodestimatorsintheplotsintermsofconvergency,e.g.theRestrictedMonteCarloEstimator,ortheImportantSamplingEstimatorfromsingleNormal.NotetheFractionEstimatorissimilartotheRestrictedMonteCarloEstimator,sincebothuseMonteCarloIntegrationoveraselecteddomainandcalculatetheevidencebydividingtheMCIbytheratiooftheselecteddomaintothefulldomain.ForaselectedVn,andthelockedrelationNMCI=NMCMC,theprecisionsoftheestimationofGnandFnbothincreaseintheorderofNMCMC. NotethevalueoftheestimationoftheFractionEstimatorisdifferentfromthatestimatedin[ 15 ]:10)]TJ /F5 7.97 Tf 6.59 0 Td[(34.97'1.110)]TJ /F5 7.97 Tf 6.58 0 Td[(35fortheFractionEstimatorv.s.410)]TJ /F5 7.97 Tf 6.59 0 Td[(35.ThedifferencecanbecausedbythebiasedestimationofFn.Indeed,forthistestof 111

PAGE 112

Figure4-12. Transformed!formodeselection.MarginalposteriorPDFof!innewrangewiththemodeinthecenter.Since!isperiodicparameter,itisconvenienttotransform!toanewrangewiththemodeinthecenter. theFractionEstimator,IuseaseparatelygeneratedMCMCsamples,whichismightdifferentfromtheMCMCsamplesusedin[ 15 ].ThecalculationofFnisbasedontheassumptionthattheMCMCsamplesformtheposteriorPDFofthemodelparametersareindependentandfullycovertheparameterspace.However,inpractice,theMCMCsamples(orstates)areautocorrelated,andthedifferentinputofmodelparameterinitialvaluesandscalesofCTPDFmightcausetheMCMCsamplesconcentratedaroundonemodewhileneglectingothermodes,resultinginanover-estimationofFn.AndfailingtoproperlydiscardtheBurn-instageoftheMCMCsamplescanalsoleadtoanunderestimationofFn 112

PAGE 113

Figure4-13. TransformingM0formodeselection.MarginalposteriorPDFofM0inrangeof[0,360]deg.M0isperiodicparameter.ItisconvenienttotransformM0toanewrangewiththemodeinthecenter. TheestimationerrorbetweentheFractionEstimatorandotherestimatorscanbecomparedusingFigure 4-16 and 4-18 .TheFractionEstimatorwithPCAgivesasimilarestimationerrorastheImportanceSamplingEstimatorortheRatioEstimator.NoteweuseNMCI=4NMCMCfortheMonteCarlointegration.Butinpractice,evenforoneplanetmodel,ifwecanspendonMCIthesamecomputationexpensewiththeMCMC,wecanincreaseNMCIbyoneorderofmagnitude.However,fromFigure 4-16 ,wecantelltheFractionEstimatorwithoutPCAisstillnotabletoyieldthesameprecisionastheestimationwithPCAatcurrentrelationbetweenNMCIandNMCMC. 113

PAGE 114

Figure4-14. TransformedM0formodeselection.MarginalposteriorPDFofM0innewrangewiththemodeinthecenter.SinceM0isperiodicparameter,itisconvenienttotransformM0toanewrangewiththemodeinthecenter. 4.5.4ChoiceoftheSizeofVn RecallthevarianceoftheFractionEstimatorconsistsoftwoterms:variancefromGnandvariancefromFn.BesidestherelationbetweenNMCIandNMCMC,thechoiceofVnalsotradesofffortwoterms.SoitisreasonabletoexaminetheeffectdifferentchoiceofVn.Iuse100000MCMCsamplesfromtheposteriorPDFofoneplanetRVmodelofHD88133andsetNMCI=400000.ByvaryingthesizeofVn,themodelevidenceisestimatedfordifferentpVn. InFigure 4-19 and 4-20 ,therelativeestimationstandarddeviationsareplottedversusthepVn,withoutandwithPCA,respectively.TherelativeestimationstandarddeviationfromGnmonotonicallyincreasesversuspVn,whileitmonotonicallydecreases 114

PAGE 115

Figure4-15. EvidenceestimationofoneplanetmodelforHD88133.EvidencecalculatedusingFractionEstimatorofoneplanetmodelofHD88133with17RVobservations.x-axisgivesthetotalnumberofMCMCsamplesusedforcalculation,NMCMC,andNMCI=4NMCMC;y-axisgivestheevidenceinlogscale.BottomPanel:TherelativeestimationstandarddeviationoftheFractionEstimator.Solidline:theevidenceestimatedusingPCA.Dashline:theevidenceestimatedwithoutusingPCA.BothestimationsconvergeasNMCMC'105,whiletheestimationwithPCAperformbetterwithsmallnumberMCMCsamplesthanthatwithoutPCA. forFn.ThetotalrelativeestimationstandarddeviationiscalculatedaccordingtoEquation 4 .ThePCAhelpstoreducerelativeestimationstandarddeviationfromGnbyaroundtwoorders.Therefore,withPCA,therelativeestimationstandarddeviationsfromGnandfromFnarecomparableandthereisabroadrangeforpVnfrom0.1to0.8tokeepthetotalrelativeestimationstandarddeviationunder1percent.Thatis,aslongaswechoosepVnwithinthisoptimizedregion,theestimationofthemodelevidencecan 115

PAGE 116

Figure4-16. EvidenceestimationuncertaintycalculatedusingFractionEstimatorofoneplanetmodelofHD88133.Basedon17RVobservations.x-axisgivesthetotalnumberofMCMCsamplesusedforcalculation,NMCMC,andNMCI=4NMCMC;y-axisgivestheestimationstandarddeviationrelativetotheestimation.BottomPanel:TherelativeestimationstandarddeviationoftheFractionEstimator.Solidline:estimationusingPCA.Dashline:theevidenceestimatedwithoutusingPCA.WiththevalueofpVn=0.22,theestimationuncertaintywithPCAisoveroneorderofmagnitudesmallerthanthatwithoutPCA. beyieldedwithsufcientlysmallerror.WithoutPCA,therelativeestimationstandarddeviationfromGnrapidlydominatesoverFnresultinginthatthetotalrelativeestimationstandarddeviationincreasesover10percentafterasmalldiptoaround4percentwhentherelativestandarddeviationofGnandFnarecomparable. Inthecasewithhigherdimensionality,forthesameNMCMCandpVn,therelativeestimationstandarddeviationfromFnkeepsthesame.However,toincludethesame 116

PAGE 117

Figure4-17. Comparisonofseveralestimatorsofthemarginalposteriorprobabilityforsingleplanetmodelsandthe17RVobservationsofHD88133.Estimatorsofthemarginalprobabilityforaone-planetmodel.Thelinestylesindicatethealgorithmsusedforeachestimator:RestrictedMonteCarlo(top,solid),PartialLinearization(top,long-dash),HarmonicMean(top,short-dash),Gelfand&DeywithPartialLinearization(top,dotted),ImportanceSamplingfromsingleNormal(bottom,dotted),ImportanceSamplingfromMixtureofNormals(bottom,solid),Gelfand&DeywithMixtureofNormals(bottom,long-dash),RatioEstimator(bottom,short-dash).ReprintedbypermissionfromFord,E.B.andGregory,P.C.2007,StatisticalChallengesinModernAstronomyIV(Figure2,leftpanel). 117

PAGE 118

Figure4-18. Comparisonofseveralestimatorsofthemarginalposteriorprobabilityforsingleplanetmodelsandthe17RVobservationsofHD88133.Internallyestimatedstandarddeviationofeachestimator.Thelinestylesindicatethealgorithmsusedforeachestimator:RestrictedMonteCarlo(top,solid),PartialLinearization(top,long-dash),HarmonicMean(top,short-dash),Gelfand&DeywithPartialLinearization(top,dotted),ImportanceSamplingfromsingleNormal(bottom,dotted),ImportanceSamplingfromMixtureofNormals(bottom,solid),Gelfand&DeywithMixtureofNormals(bottom,long-dash),RatioEstimator(bottom,short-dash).ReprintedbypermissionfromFord,E.B.andGregory,P.C.2007,StatisticalChallengesinModernAstronomyIV(Figure2,rightpanel). 118

PAGE 119

fractionofMCMCsample,Vnincreasesintheorderofdimensiondim,increasingtheestimationvarianceofGn.SettingalargerrelationnumberbetweenNMCIandNMCMC,whichcanbeproportionaltothenumberofthemodelparameter,canreduceestimationerrorbyafactorofp dim.Butasthenumberofthemodelparametersrising,theoptimizedregionofpVnisshortenbyrisingestimationerrorcausedbyGnterm,andtheshapeofestimationerrorcurvewilltransformfromasinFigure 4-20 toFigure 4-19 .ThiscanbefurtherillustratedinhigherdimensionalexamplesinChapter5. Figure4-19. RelativeestimationstandarddeviationcalculatedwiththeFractionEstimatorwithoutPCA.PlottedversuspVn.Solidline:thetotalrelativeestimationstandarddeviation,calculatedusingEquation 4 .Dashline:therelativeestimationstandarddeviationcausedbyGn.Dot-dashline:therelativeestimationstandarddeviationcausedbyFn.Whenthedashanddot-dashcurvecross,thesolidcurvereachesitsminimumvalue,whichmeanwhentheestimationuncertaintiescausedbyGnandFnarecomparable,thechoiceofpVnisoptimized. 119

PAGE 120

Figure4-20. RelativeestimationstandarddeviationcalculatedwiththeFractionEstimatorusingPCA.PlottedversuspVn.Solidline:thetotalrelativeestimationstandarddeviation,calculatedusingEquation 4 .Dashline:therelativeestimationstandarddeviationfromGn.Dot-dashline:therelativeestimationstandarddeviationcausedbyFn.Whenthedashanddot-dashcurvecross,thesolidcurvereachesitsminimumvalue,whichmeanwhentheestimationuncertaintiescausedbyGnandFnarecomparable,thechoiceofpVnisoptimized. 4.6SummaryandDiscussion TheprocessofcalculatingthemodelevidenceusingtheFractionEstimatorforcanbesummarizedasbelow: 1)ObtainMCMCsamplesoftheposteriorPDFofthemodelparametersbasedontheobservations.ReducetheautocorrelationamongthesamletsbytakingsampleseverynthstepanddiscardthesamplesintheBurn-instage,usuallyacertainnumberofthesamplesthatarerstgenerated. 120

PAGE 121

2)FindtheglobalmodeaccordingthemarginalposteriorPDFsofthemodelparameter,whichcanbeplottedusinghistogramsoftheMCMCsamples.Iftherearemorethanoneoutstandingmodes,onecanusetheFractionEstimatorforeachmodeandcombineeachGnandFntoyieldthevalueofmodelevidence. 3)UsingPrincipleComponentsAnalysis(PCA)toreducethevarianceoftheposteriorPDFinchoosesub-domainVnintheparameterspace. 4)DeterminetherelationbetweenNMCIandNMCMC.Intermsofcomputationexpense,NMCIischeaperthanNMCMC,andusuallywecansetNMCI4NMCMC. 5)DeterminethesizeofVnorpVn.WhentheestimationerrorfromGnandFnarecomparable,thetotalestimationerrorisminimized. 6)EstimateGnusingMonteCarlointegrationoverVn;estimateFnbycountingthenumberofsamplesthatareincludedinVn.Byincreasingthetotalnumberofsamplesusedforcalculation,onecanreducetheestimationerror(relativeestimationstandarddeviation). 7)TherelativeestimationstandarddeviationoftheFractionEstimatorcanbecalculatedusingEquation 4 TheFractionEstimatorhasbeenprovensuccessfulinoneplanetmodelwith17RVobservationsofHD88133.However,thereareafewconcernsonthismethodoneshouldnote. First,theMCMCsamplescandeterminetheaccuracyoftheFractionEstimator.AdifferentvaluecanbeyieldedwithdifferentMCMCsample.SinceweusetheMCMCsamplestoestimateuncertaintyofthemeasurementofthemodelparameters,itisreasonabletocalculatethemodelevidencebasedonthesameMCMCsample.Sothepremiseoftheaccuracyofthemodelevidence,isaMCMCsamplesequencethataccuratelysamplefromtheposteriorPDFofthemodelparameter,whichisalsotherequirementforestimatingtheuncertaintyofparameterestimation. 121

PAGE 122

Second,PCAcanonlyeliminatethelinearcorrelationbetweenparameters.Forparametersthathavenolinearcorrelations,theparametersarestillcorrelatedevenafterbeingtransformedusingPCA,resultinglargevarianceoftheposteriorPDFinVn.Forexample,foreand!,itisusefultoconsideralternativeparametersform,i.e.,esin!andecos!toreducethecorrelationbetweeneand!. Third,thenatureoftheFractionEstimatoristoavoidhighdimensionalintegrationinfulldomainbychoosingasmallsub-domainintheparameterspaceandusingtheMCMCsamplestoestimatethefractionofthechosensub-domain.However,theintegrationinthesub-domainstillsuffersfromthecurseofhighdimension,especiallyformodelswithmorethanoneplanet.Sofurtherimprovementoftheintegrationinthesub-domainshouldbeconsidered.1)QuasiMonteCarlointegrationdrawssamplesmoreevenlydistributedintheparameterspacethantheplainMonteCarlointegrationmethod,soastogivebetterestimationofGn.2)TheposteriorPDFmightbebetterapproximatedwithaGaussianfunctionaroundthemodethaninthefullparameterspace,allowingustouseImportantSamplingorotheralgorithmstoestimateGn. 122

PAGE 123

CHAPTER5APPLICATIONS 5.1BackgroundandMotivation InChapter4,IpresentedtheFractionEstimatorforcalculatingthemodelevidencebasedontheMCMCsamplesfromtheposteriorPDFofcertainmodels(howmanyplanets)withtheRVobservationsofexoplanetarysystems.TheprocessofusingtheFractionEstimatorissummarizedbelow. First,theFractionEstimatorrequiresgoodqualityMCMCsamplesoftheposteriorPDFgiventhemodelwhoseevidenceisbeingcalculated.GoodqualitymeanstheMCMChasfullyexploredparameterspaceandhasconvergedtotheposteriorPDF.Fordifferentmodels,afewoftechniquescanbeappliedtoimprovetheefciencyofMCMC.Forthesetechnicalissues,onecanrefer[ 14 ]forexoplanetarymodels. Second,giventheMCMCsamplesdrawnfromtheposteriorPDFofthemodelofinterest,theFractionEstimatorchoosesasub-domaininthemodelparameterspaceandcalculatestheplainMonteCarlointegrationoverthissub-domain.TheratiooftheMonteCarlointegrationtothefractionofthenumbersoftheMCMCsamplesinthissub-domainoverthatinthefulldomainistheestimationofthemodelevidence.Inordertoreducethevarianceoftheestimation,theFractionEstimatorrequiresasub-domainwithsmallvolume,largeaveragevalueofposteriorPDFandsmallvarianceoftheposteriorPDF,leadingtheoptimizedchoiceofthesub-domaintobearoundtheglobalmode. Third,anotherimportantstepistousethetechniqueofspacetransformationusedinprinciplecomponentsanalysis(PCA)todeterminetheorientationofthesub-domaininparameterspacebyreducingthelinearcorrelationbetweenmodelparameters.ThiscanmakethebestuseoftheMCMCsamplesordomainaroundthemode,whichcansignicantlydecreasethevarianceoftheposteriorPDFwithin. 123

PAGE 124

Essentially,thecurrentalgorithmofFractionEstimatorchoosesasub-domainoftheparameterspace,usesMCMCsamplescountingtocalculatethenormalizedintegrationofthesub-domain,andusesplainMonteCarlofortheun-normalizedintegration.Theratioofthetwointegrationsyieldsthenormalizationfactor,i.e.,themarginallikelihood,orthemodelevidence.InChapter4,ithasbeenshownthattheFractionEstimatorisabletoefcientlyestimatetheevidenceofoneplanetmodelbasedonRVobservationswithMCMCsampling.Theoneplanetmodelwith17RVobservationsusedforatesthasonly7parameters.Therefore,itisnecessarytoapplytheFractionEstimatorforsystemswithmoreparameters,i.e.,moreplanets,toexaminetheestimator'sperformance.Thereareseveralissueswecanexpectandshouldpayparticularattentionto. Whenthenumberofmodelparameterincreases,i.e.,thedimensionoftheparameterspaceincreases.First,itrequireslongerMCMCchainstoachievetheacceptableconvergencewithfullexplorationintheparameterspace.Inordertoreducetheauto-correlationbetweenthestateoftheMCMCchains,thefractionofthestatesabandonedduringthegenerationofMCMCshouldalsoincreaseproportionallytothenumberofthemodelparameters.Second,inthehighdimensionalspace,eveninasub-domainaroundtheglobalmode,theplainMonteCarlointegrationwillalsoencounterthecomputationdifcultyduetothecurseofhighdimension.EventhoughthereisnoexplicitdimensionalityterminthevarianceofplainMonteCarlointegration,itwillstillincreaseduetothevarianceofthePDFwithinthesub-domainused.Third,sinceitisdifculttovisualizetheshapeofahigh-dimensionalposteriorPDF,weusethemarginalposteriorPDFforeachmodelparametertodeterminethelocationoftheglobalmode;thesizeofthesub-domainiscontrolledbythechoiceoftherangesofthemarginalposteriorPDFforeachmodelparameter,andasthedimensionincreases,therangesofthesub-domainineachdimensionwillalsoincreaseinordertominimizethe 124

PAGE 125

varianceoftheFractionEstimator(acombinationofthevariancesofthenormalizedintegrationandun-normalizedintegrationofthesub-domainfromFnandGn). InordertoinvestigatetheperformanceoftheFractionEstimatorinhigherdimensionalscenariosandthelimitationcausedbythetechnicalissuesmentionedabove,inChapter5,IexaminetheestimatorwithbothsimulatedRVdataandrealobservedRVdataofexoplanetarysystems.Insection 5.2 ,IchooseaweaklyinteractingtwoplanetsystemgeneratedinChapter3usingN-bodysimulation,anditsmockRVdata,tocalculatetheevidenceofKeplerianmodelsforone,twoandthreeplanets,respectively,usingtheFractionEstimator.Thenbycomparingthemodelevidences,IshowtheprocessofaBayesianmodelselection.Insection 5.3 ,Ichooseseveralrealsystemswithknownplanets,includingtwosystemswiththreeplanetsandonesystemwithfourplanets.TheevidencesofknownnumberofplanetsarecalculatedusingtheFractionEstimatorforexaminingtheperformanceoftheestimator.PleasenoteonlytheKeplerianmodelisconsidereddespitethenumberoftheplanetsassumedintheplanetarysystem. Furthermore,theBayesianmodelselectionmethodispotentiallyableforinstrumentnoiseanalysis.TheFractionEstimatorcanbeappliedfortheBayesianmodelselectionsinceitdoesnotdependonanyspecicmodel.InChapter5,IalsodescribetheattemptofapplyingBayesianmodelselectionontheMARVELSinstrument.AsummaryofthisdissertationisgivenintheendofChapter5alongwiththeprospectivefutureworkforimprovingthisresearch. 5.2ModelSelectionwithSimulatedData IntheprojectintroducedinChapter3,asetoftwoplanetsystemsaregeneratedwithN-bodysimulation.Sincewehavetheaccurateinformationaboutparametersofthecorrectmodel,i.e.,twoplanetmodel,itisconvenienttouseoneofthemtoexaminetheFractionEstimatorandapplyBayesianmodelselection.ForthepurposeofexaminingtheFractionEstimator,wedonothavetoconsidertheplanet-planet 125

PAGE 126

interaction.Therefore,asystemthatcanbeaccuratelyttedwithKeplerianmodelismoresuitable.Asystem'090030br',withinnerplanetof0.9MJandouterplanetof0.3MJbeforeevolvinginto2:1meanmotionresonanceischosen. 150mockRVobservationsover15yearsforthissystemareusedforMCMCanalysis.InChapter3,theKeplerianttingyieldsaccurateestimationsforthissystemastheinteractionbetweenthetwoplanetsisnegligible.Table 5-1 liststheorbitalparameters(P,K,e,!,tp),oftheco-planartwoplanetsystem090030brmeasuredduringtheN-bodysimulation,wheretpisthetimeofperiastronpassage.Thedifferentformoforbitalparametersarethentransformedintotheform(P,K,e,!,M0)astheinitialinputsforMCMC. Table5-1. Orbitalparametersofthetwoplanetsystem090030br.TheparametersweremeasuredduringtheN-bodysimulation.Andthemockobservationtimerangesfrom6538.60295to8345.39832,withoneobservationaboutevery20days.ReferChapter3fordetailsabouttheparametermeasurement. InnerPlanetOuterPlanet P(day)128.26267.46K(m/s)36.259.46e0.00560.0047!(deg)235.68253.82tp(day)719.0010011.00 5.2.1OnePlanetModel UsingtheMCMCalgorithmsdevelopedby[ 14 ],theMCMCsamplesoftheposteriorprobabilitydistributionofoneplanetmodelparametersisgenerated.TheMCMCsamplingtakesstepsin(logP,logK,e,!,M0,Cj,jitter),whereCjisthejthreferencevelocity,andforthisdata,thereisonlyonevelocityreferencesoj=0.Theoutputmodelparametersarealso(P,K,e,!,M0,C,jitter).FromthehistogramsoftheMCMCsamples,themarginalposteriorPDFofthemodelparametersareconcentratedaroundtherealvalueoftheinnerplanet.Itisareasonableresultsincethemassoftheinnerplanetisthreetimesasthatoftheouterplanet,whilethesemi-majoraxisoftheinner 126

PAGE 127

planetislessthanhalfofthethatoftheouterplanet,theinnerplanetcausesaboutoneorderofmagnitudemorechangeofthecentralstar'svelocitythantheouterplanetdoes.Thesignalfromtheouterplanetwillbereectedintheincreasedjitterthatissupposedtobeverysmall. TheevidenceofoneplanetmodelcalculatedusingtheFractionEstimatorisshowninFigure 5-1 ,plottedversusthetotalnumberofMCMCsamplesusedformodelevidencecalculation.ThesamepriorPDFfromTable 4-1 isusedforthemodelparameters.TherelationbetweenthedrawsfortheMCIandthetotalnumberofMCMCsamplesisthesamewiththatinChapter4:NMCI=4NMCMC.AsthenumberofMCMCsamplesincreases,theestimationoftheevidenceofoneplanetmodelforthesystem090030brwithRVobservationsbeginstoconvergetothenalvalueofabout10)]TJ /F5 7.97 Tf 6.58 -.01 Td[(234.38.Thechoiceofthesub-domainmakesthefractionoftheMCMCsamplescountedthesub-domainoverthefullMCMCsamplespVn0.3,asshowninFigure 5-2 .WhenNMCMC'105,theestimationofthefractionconverges.InFigure 5-3 ,thestandarddeviationsoftheevidenceestimation(causedbyFn,Gn,andtotal)areshown.TherelativestandarddeviationsdecreasewiththeincreaseofthetotalnumberofMCMCsampleused,wheretherelativestandarddeviationscausedbyFndominatesthatbyGn.BecausethevarianceofGnincreasesasthesizeofthesub-domainincreasewhilethatofFnincreaseasthesizeofthesub-domain(Moreaccurately,thefractionpVn.Butsincewechoosethesub-domainaroundthemode,pVnispositivecorrelatedwiththesizeofthesub-domain)increases,theoptimizedsizeofthesub-domainiswhatmakestherelativevariancesfromGnandFnarecomparable.Therefore,Figure 5-3 impliesthatthechoiceofsizeofthesub-domainorthefractionpVnissmallerthantheoptimizedvalue.WhenNMCMC'105theevidenceestimationcanachieveaccuracylessthan1percent. 5.2.2TwoPlanetModel Forthetwoplanetmodel,theparametersfortheMCMCsamplesare(Pp,Kp,ep,!p,M0,p,Cj,jitter),wherep=0,1andCj(j=0)isthereferencevelocity,thereforethere 127

PAGE 128

Figure5-1. EvidenceofoneplanetmodelofaN-bodysimulatedtwo-planetsystem.CalculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamplesusedforcalculation.150RVobservationsover15yearsareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC. are12parameters.WithconvergedMCMCsamples,wecanusehistogramtoplotthemarginalposteriorPDFsofeachmodelparametersasinFigure 5-4 .Fromthetoptothebottomandfromthelefttotheright,eachpanelshowstheposteriordistributionoftheorbitalparametersfortheinnerandouterplanet,respectively.ThelasttwopanelsshowtheRVreferenceC,andthelogvalueoftheunknownnoiseincludingjitter.ThisisthesteptoestimatetheuncertaintyoftheorbitalparameterestimationusingMCMC. Figure 5-5 showsthemarginalposteriorPDFineachdimensioninthetransformedspaceusingPCAtechniquesandtheselectedrangeofthesub-domainineachdimensionoftheparameterspace.Theregionbetweenthetwoverticallinesin 128

PAGE 129

Figure5-2. FractionofMCMCsamplesusedforevidencecalculationofoneplanetmodelforaN-bodysimulatedtwoplanetsystem.ThefractionofMCMCsamplesincludedinthesub-domainforplainMonteCarlointegration,plottedversusthetotalnumberofMCMCsamplesusedforcalculation.TheMCMCsamplesisbasedononeplanetmodelofaN-bodysimulatedtwoplanetsystemwith150RVobservationsover15years. eachhistogrampanelistheregionofthehypercubeVnprojectedineachdimension.Currently,IusethemethodtodecidetherangeasusedinChapter4:setarelativeheight,say,0.7,toitspeakvalue(usuallythemode)ofthemarginalposteriordistribution.Searchonbothsidesofthepeak(themode)fortwopositionswherethevaluesofthePDFarecorrespondedtotherelativeheight,andthetwopositionsarethelowerandupperboundoftheregionofthesub-domainprojectedinthisdimension.Thatis,largerrelativeheightwillyieldsmallerrangearoundthemodeineachdimension,thereforesmallerVnandpVn. 129

PAGE 130

Figure5-3. RelativestandarddeviationoftheevidenceestimationofoneplanetmodelforaN-bodysimulatedtwo-planetsystem,calculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamplesused.Solidline:thetotalrelativeestimationstandarddeviation.Dashline:relativeestimationstandarddeviationcausedbyGn.Dotted-dashline:relativeestimationstandarddeviationcausedbyFn.150RVobservationsover15yearsareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .NMCI=4NMCMC. TheevidenceofthetwoplanetmodelcalculatedusingtheFractionEstimatorisshowninFigure 5-6 plottedversusthetotalnumberofMCMCsamples,whereNMCI=4NMCMC.ThevalueoftheevidencestartstoconvergeafterthetotalnumberofMCMCsamplesusedforthecalculationlargerthan105,andthenalestimatedvalueoftheevidenceisabout10)]TJ /F5 7.97 Tf 6.59 0 Td[(173.15whenNMCMC'105.ThefractionoftheMCMCsamplescountedthesub-domainoverthefullMCMCsamplesconvergestoaroundpVn=0.26,asshowninFigure 5-7 .Therelativestandarddeviationoftheevidenceestimationis 130

PAGE 131

Figure5-4. PosteriorPDFsoftwoplanetmodelparameters.BasedonthemockRVobservationsofaN-bodysimulatedtwoplanetsystem.Fromthetoptothebottomandfromthelefttotheright,eachpanelshowstheposteriordistributionoftheorbitalparametersfortheinnerandouterplanet,(Pp,Kp,ep,!p,M0,p),respectively,wherep=0,1.ThelasttwopanelsshowtheRVreferenceC,andthelogvalueoftheunknownnoiseincludingjitterLnJitter. showninFigure 5-8 .WhenNMCMC106,thetotalrelativestandarddeviationisabout10percent.ThedominanceofGnfortherelativeestimationstandarddeviationimpliesthatforgivenrelationconstantbetweenNMCIandNMCMC,thesizeofVn(correspondingpVn=0.26)islargerthantheoptimizedchoiceofpVn.OnewaytoreducetherelativestandarddeviationcausedbyGnistoincreasethedrawsoftheplainMonteCarlointegration.Indeed,asthemodelparameterincrease,theconstantconstrainingNMCIandNMCMCshouldincreasesinceittakesmoreevaluationsofthemodelforoneMCMCsample.However,giventhedifferencebetweentheevidencesoftwomodels(10)]TJ /F5 7.97 Tf 6.59 0 Td[(234.38 131

PAGE 132

Figure5-5. PosteriorPDFsoftwoplanetmodelparameterstransformedintoPCAspace.BasedonthemockRVobservationsofaN-bodysimulatedtwoplanetsystem.Ineachpanel,thetwoverticallinesintersectthePDFwiththesamerelativeheight.TherangebetweenthetwoverticallinesineachhistogramistheregionofthehypercubeVnineachdimension. fortheoneplanetmodelv.s.10)]TJ /F5 7.97 Tf 6.59 0 Td[(173.15forthetwoplanetmodel)for090030br,thecurrentestimationissufcientlyaccuratetotellthetwoplanetmodelismorefavoredbytheobservation. 5.2.3ThreePlanetModel Forthethreeplanetmodel,themodelparametersare(Pp,Kp,ep,!p,M0,p,Cj,jitter),wherep=0,1,2andCj(j=0)isthereferencevelocitythereforethereare17modelparameters.TheevidenceofthethreeplanetmodelcalculatedusingtheFractionEstimatorisshowninFigure 5-9 ,plottedversusthetotalnumberofMCMCsamples.Thevalueofthethreemodelevidenceconvergestoaround10)]TJ /F5 7.97 Tf 6.59 0 Td[(235.81as 132

PAGE 133

Figure5-6. EvidenceofthetwoplanetmodelofaN-bodysimulatedtwoplanetsystem.CalculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamplesusedforthecalculation.150RVobservationsover15yearsareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC. NMCMC'106.AndthefractionoftheMCMCsamplescountedthesub-domainoverthefullMCMCsamplesconvergestopVn=0.14,asshowninFigure 5-10 .TherelativeestimationstandarddeviationoftheevidenceisshowninFigure 5-11 .Similartothecaseofthetwoplanetmodel,thedashline(Gn)isoverlappedwiththesolidline(Total).ThedominanceoftherelativeestimationstandarddeviationfromGnoverFnimpliesthatthechosenVnislargerthantheoptimizedvalueofpVn,whichminimizesthetotalrelativeestimationstandarddeviation.AlthoughwecanincreaseNMCIinsteadofusingtherelationNMCI=4NMCMC,thelargevarianceoftheGnmainlycausedbytheincrease 133

PAGE 134

Figure5-7. FractionofMCMCsamplesforthetwoplanetmodelevidenceofaN-bodysimulatedtwoplanetsystem.ThefractionofMCMCsamplesincludedinthesub-domainforplainMonteCarlointegration,plottedversusthetotalnumberofMCMCsamplesusedforthecalculation.TheMCMCsamplesisbasedontwoplanetmodelofaN-bodysimulatedtwoplanetsystemwith150RVobservationsover15years. ofthespacedimension.However,comparewiththetwoplanetmodel,itisstillobviousthatgiventheobservations,thetwoplanetmodelismorefavored. 5.2.4DiscussionandModelSelection IfwecomparetheevidenceestimationforthreemodelsasshowninFigure 5-1 5-6 and 5-9 ,itisclearlythatthetwoplanetmodel,whichisthetruemodelforthemockRVobservations,ismorelikelythantheoneandthethreeplanetmodel,sincetheevidenceforeachmodelisE1:E2:E3=10)]TJ /F5 7.97 Tf 6.59 0 Td[(234.38:10)]TJ /F5 7.97 Tf 6.59 0 Td[(173.15:10)]TJ /F5 7.97 Tf 6.59 0 Td[(235.81.IttakesmoreMCMCsamplestoconvergeforthemodelwithmoreplanets,i.e.,moremodelparameters, 134

PAGE 135

Figure5-8. RelativestandarddeviationoftheevidenceestimationoftwoplanetmodelforaN-bodysimulatedtwoplanetsystem.CalculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamplesused.Solidline:thetotalrelativeestimationstandarddeviation.Dashline:relativeestimationstandarddeviationcausedbyGn,overlappingwiththesolidline.Dotted-dashline:relativeestimationstandarddeviationcausedbyFn.150RVobservationsover15yearsareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC. sinceasthenumberofthemodelparameterincreases,theFractionEstimatorrequireslargernumberofMCMCsamplestoobtainthesufcientlysmallestimationvariance. ThetrickoftheFractionEstimatorliesinthechoiceofthepVn.AccordingtoEquation 4 ,whentherelationbetweenNMCMCandNMCIisxed,aswehavedoneinthisexample,theonlyfactordeterminingtheestimationvarianceisthechoiceofthesub-domain.InChapter4,wehavediscussedthattheoptimizedpositionofthesub-domainisaroundthemode,sincethetargetdistributiong(x)hasthepeak 135

PAGE 136

Figure5-9. EvidenceofthreeplanetmodelofaN-bodysimulatedtwoplanetsystem.CalculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamplesusedforcalculation.150RVobservationsover15yearsareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC. valueanditresultsinasub-domainwithlargerpVn.AndweusethePCAtechniquestotransformtheparameterspacetooptimizetheorientationofthesub-domain. Alargersub-domainaroundthemodeisrequiredtoreducethevariancecausedbyFn.However,thiscancausealargervarianceofg(x)inthesub-domain,andtheaveragevaluemightbecomesmaller.InthatForeachmodel,thefractionofthetotalMCMCsamplesincludedinthesub-domain,orpVnis0.3,0.25and0.14(Figure 5-2 5-7 and 5-10 ),respectively.ThevaluesofpVnwithsimilarorderofmagnitudearechosenforconvenientcomparison,i.e.,thevariancecausedbyFnfromeachmodelisclose.Sothedifferencebetweentheevidenceestimationvariancesismainlycausedby 136

PAGE 137

Figure5-10. FractionofMCMCsamplesusedforevidenceestimationofthreeplanetmodelofaN-bodysimulatedtwoplanetsystem.ThefractionofMCMCsamplesincludedinthesub-domainforplainMonteCarlointegration,plottedversusthetotalnumberofMCMCsamplesusedforcalculation.TheMCMCsamplesisbasedonthreeplanetmodelofaN-bodysimulatedtwoplanetsystemwith150RVobservationsover15years. Gn.ComparingFigure 5-3 5-8 and 5-11 ,thevariancecausedbyGndominateforallthreemodels.Particularly,forthethreeplanetmodel,thevariancecausedbyGnistwoorderofmagnitudegreaterthanthatcausedbyFn.ThisisbecauseinordertoincludeenoughMCMCsamplesinthesub-domaintoachievethechosenpVn,thesizeofthesub-domainwillincreaseintheorderofthepowerofthenumberofmodelparameters,thereforethevarianceofg(x)willincreaseinthewaywelackaclearknowledgeabout,sinceitishardtodeterminetheshapeofg(x),ortheposteriorPDF. 137

PAGE 138

Figure5-11. RelativestandarddeviationoftheevidenceestimationofthreeplanetmodelforaN-bodysimulatedtwoplanetsystem.CalculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamplesused.Solidline:thetotalrelativeestimationstandarddeviation.Dashline:relativeestimationstandarddeviationcausedbyGn.Dotted-dashline:relativeestimationstandarddeviationcausedbyFn.150RVobservationsover15yearsareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC. TheevidenceestimationcanbemoreaccuratewithlargernumberofMCMCsamplesoramoreaggressivechoiceoftheconstantrelatingNMCIandNMCMC,e.g.largerthan4.Butforthissystem,andformostastronomyproblemsthathavenohighaccuracyrequirement,thecurrentaccuracyissufcientforBayesianmodelselection.However,inthecaseofhighaccuracyrequirementtocomparetwoclosemodels,wecanadjustthepVntoachievetheoptimalsizeofthesub-domain.Forexample,fortheoneplanetmodel,ifwerequire5percentaccuracywithNMCMC105,wecanset 138

PAGE 139

Table5-2. ThreeselectedplanetarysystemsforexaminingtheFractionEstimator.Wechosethreeplanetarysystemswiththree,three,andfourknownplanetsformodelevidenceestimationusingtheFractionEstimator.NumberoftheRVobservations,velocityreferencesandthetimespanoftheobservationsarealsolisted. SystemNpObservationsReferenceObsTimeSpan(days)And32681627461Virgo320621571Ara423432768 smallervalueofpVnsincethecurrentvariancecausedbyFnisfarsmallerthantherequirement.Forthethreeplanet,however,thereisnotmuchspaceforimprovementbyadjustingpVnsincetherelativestandarddeviationcausedbyFnisaboutseveralpercentwhenNMCMC105.Wehaveto1)increasetheconstantconstrainingNMCIandNMCMC,or2)improvethealgorithmtocalculatetheun-normalizedintegration,or3)justincreasetheNMCMC. 5.3EvidenceCalculationwithRealObservations ToexaminetheperformanceoftheFractionEstimator,Ichoosesomepreviouslystudiedplanetarysystemswithawelldeterminednumberofplanets.Table 5-2 liststhenameofthesystems,alongwiththenumberoftheknownplanets,numberoftheRVobservations,andnumbersofobservatoriesforeachsystem.TheMCMCsamplesaregeneratedbyFord[ 14 ]withcorrespondingmodel(numberofplanets)listedinTable 5-2 ,andthesameRVobservationsreferredin[ 14 ]andusedfortheMCMCareusedformodelevidenceestimation.ThiscanguaranteebetterMCMCsamplesthereforereducetheeffectfromthequalityofMCMCsamplesontheevidenceestimation.TherelationNMCI=4NMCMCisstillusedforevidencecalculations.AndthepriorPDFsformodelparametersarethesameaslistedinTable 4-1 5.3.1AndwithThreePlanetsandOneVelocityReference AsofAugust,2012,Andromedahasbeenconrmedwithfourorbitingplanets.ButsincethisresearchusedtheRVdataandtheMCMCsamplesby[ 14 ]forcharacterizing 139

PAGE 140

theuncertaintyoftheorbitalparameters,onlythreeplanetshadbeenconrmedatthemomentandthefourthlongperiodplanetwasdetectedin2010[ 4 ]. ForsystemAndwiththreeplanetmodelandonevelocityreference,thereare17modelparameters.Themodelevidence,thefractionofthenumberofMCMCsamplesincludedinthesub-domain,andtherelativeestimationstandarddeviationcalculatedwiththeFractionEstimatorareplottedversusthetotalnumberofMCMCsamplesinFigure 5-12 5-13 and 5-14 .Theevidenceconvergestoaround10)]TJ /F5 7.97 Tf 6.58 0 Td[(500.295withpVn0.057.AsshowninFigure 5-12 ,thetotalrelativeestimationstandarddeviationislessthan1percentwith106MCMCsamples,andthosecausedbyGnandFnarecomparable,whichmeansthesizeofthesub-domainisveryclosetotheoptimal. 5.3.261VirwithThreePlanetsandTwoVelocityReferences Thesystem61Virhasthreeknownplanets,andtheRVdataweusecomefromtwoobservatoriesthereforetwovelocityreferences.Forthismodel,themodelevidence,thefractionofthenumberofMCMCsamplesincludedinthesub-domain,andtherelativeestimationstandarddeviationarecalculatedwiththeFractionEstimatorandplottedinversusthetotalnumberofMCMCsamplesFigure 5-15 5-16 and 5-17 .Theevidenceconvergestoaround10)]TJ /F5 7.97 Tf 6.59 0 Td[(223.90asNMCMC'105,andpVnstartstoconvergeto0.09asNMCMC'105.Thetotalrelativeestimationerrorislessthan1percentaftertotalnumberofMCMCsamplesislargerthan105.ThecomparablevaluesoftherelativeestimationerrorfromGnandFnindicatesthisisanotherexampleofgoodchoiceofpVn.NoteduetolimitedMCMCsampleforthissystem,theconvergencyofthecurvesinFigure 5-15 and 5-16 isnotveryobvious,butalsonotetherangethecurvesvaryisactuallysmall. 5.3.3ArawithFourPlanetandThreeVelocityReferences ThesystemArahas4knownplanets.TheMCMCsamplesweregeneratedbasedonfourplanetmodelwiththreevelocityreferences,therefore24modelparameters.Forthismodel,theevidence,thefractionofthenumberofMCMCsamplesincludedinthesub-domain,andtherelativeestimationstandarddeviationversusthetotalnumberof 140

PAGE 141

Figure5-12. EvidenceofthreeplanetmodelofAnd.CalculatedusingtheFractionEstimator,isplottedversusthetotalnumberofMCMCsamplesused.268RVobservationsareused.ThepriorPDFsformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC.Theevidenceestimationconvergestoaround10)]TJ /F5 7.97 Tf 6.59 0 Td[(500.295whenNMCMC'106. MCMCsamplesareplottedinFigure 5-18 5-19 and 5-20 ,respectively.Unfortunately,theevidencedoesnotconvergewhenpVn0.013with106MCMCsamples.TheestimationerrorismainlycausedbythetermofGn,asillustratedinFigure 5-20 .Thisimpliesthateventhoughonlyabout1percentofMCMCsamplesareincludedinthesub-domainforevidencecalculation,thecorrespondingpVnisstilllargerthantheoptimizedvalueduetothehighdimensionalityfortheMCI. 5.3.4Discussion ThethreesystemsAnd,61VirandArawithrealRVobservationdatahavebeenpreviouslycarefullystudiedandtheirnumbersoforbitingplanetsgiventhe 141

PAGE 142

Figure5-13. FractionofMCMCsamplesincludedinthesub-domainforAnd,pVn,forplainMonteCarlointegrationofthreeplanetmodel,plottedversusthetotalnumberofMCMCsamplesusedforcalculation.TheestimationofpVnstartstoconvergetothevalueof0.057whenNMCMC'104. detectionprecisionasof2006havebeenwelldetermined.Onlythemodelswiththewelldeterminednumberofplanetforthesesystemsareusedforthemodelevidencecalculation.SincetheMCMCsamplesofthesemodelshavebeencarefullygeneratedbyFord[ 14 ],itreducethepossibilityoftheeffectfromthelowqualityofMCMCsamples.Thenumberofthemodelparametersare17,18and24,respectively.TheFractionEstimatorcanachieveabout1percentrelativestandarddeviationforthemodelevidenceestimationwith105MCMCsamplesintheformertwocases;inthecaseof24modelparameters,theFractionEstimatorshasnotconvergewith106MCMCsamples,anditachievesabout10percentrelativestandarddeviation.With 142

PAGE 143

Figure5-14. RelativestandarddeviationoftheevidenceestimationofthreeplanetmodelforAnd.Asystemwiththreeplanet.CalculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamplesused.Solidline:thetotalrelativeestimationstandarddeviation.Dashline:relativeestimationstandarddeviationcausedbyGn.Dotted-dashline:relativeestimationstandarddeviationcausedbyFn.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC. largernumberofMCMCsamples,oroptimizedchoiceofpVn,theFractionEstimatorisexpectedtoachievetheestimationuncertaintyofafewpercent. TheadvantageoftheFractionEstimatorisitusesnoassumptionabouttheshapeoftheposteriorPDFbuttheMCMCsamplesaresampledfromtheposteriorPDF.WiththePCAtechniquesforreducingthelinearcorrelationsbetweenthemodelparameters,asub-domainaroundthemodewithoptimizedorientationcanbefoundforplainMonteCarlointegration.Thesizeofthissub-domaincanbecontrolledtoobtaintheminimumestimationuncertainty(relativeestimationstandarddeviation),wheretheestimation 143

PAGE 144

Figure5-15. Evidenceofthreeplanetmodelwithtwovelocityreferencesfor61Vir.Withtwovelocityreferences(18modelparameters)for61Vir,calculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamples.206RVobservationsover1571daysareused.ThepriorPDFsformodelparametersarelistedinTable 4-1 .NMCI=4NMCMC.Theestimationofmodelevidencevaryintherangeof10)]TJ /F5 7.97 Tf 6.58 0 Td[(224.0)]TJ /F6 11.955 Tf 11.96 0 Td[(10)]TJ /F5 7.97 Tf 6.59 0 Td[(223.80andstartstoconvergetoaround10)]TJ /F5 7.97 Tf 6.59 0 Td[(223.90asNMCMC'105. uncertaintiesfromtheGnandFnareclose.Fortheintegrationoverthesub-domain,therecouldbemoreefcientalgorithmotherthantheplainMonteCarlo.However,thecurrentalgorithmoftheFractionEstimatorcanworkwellwithfourplanetmodel. 5.4MARVELSReferenceStars TheMARVELSinstrumentishostedinthe2.5mtelescopeApachePointObservatoryasoneofSDSSIIIprojects.Ineverysingleplanetsurveyobservation,itsimultaneouslyobservesabout60targetsincludingseveralreferencestarswithatleastoneknowngiantplanet.ThosereferencestarsusuallyarebrightenoughforMARVELStogain 144

PAGE 145

Figure5-16. FractionofMCMCsamplesinsub-domainfor61Vir.ThefractionofMCMCsamplesincludedinthesub-domain,pVn,forplainMonteCarlointegration,plottedversusthetotalnumberofMCMCsamplesusedforthecalculation.TheMCMCsamplesarebasedonthreeplanetmodelwithtwovelocityreferencesfor61Vir,i.e.,18modelparameters.206RVobservationsover1571daysareused. RVmeasurementswithacceptableerrorsforscienticpurposes.IchoosecertainMARVELSreferencestarsthatalreadyhaveenoughdatapoints,withobservationscoveringatleastoneperiodoftheannouncedplanet,andthesignalscausedbyplanetorbitinglargerthantheRVmeasurementprecision(i.e.detectablebyMARVELS).Twosystemswithoneknowplanetandaunconrmedsecondplanet,HD17156andHD68988. However,duetolimitedobservationprecision,eventhoughtheMARVELSobtaineddatacoveringthepossiblesecondplanet(P=111days)ofHD17156,itssignalistoo 145

PAGE 146

Figure5-17. Relativestandarddeviationoftheevidenceestimationofthreeplanetmodelwithtwovelocityreferencesfor61Vir(18modelparameters).CalculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamples.Solidline:thetotalrelativeestimationstandarddeviation.Dashline:relativeestimationstandarddeviationcausedbyGn.Dotted-dashline:relativeestimationstandarddeviationcausedbyFn.206RVobservationsover1571daysareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .NMCI=4NMCMC.ComparablevaluescausesbyGnandFnindicatethatthechoiceofpVnisclosetotheoptimal. weakcomparedwithmeasurementuncertainty(2.3to10ms)]TJ /F5 7.97 Tf 6.59 0 Td[(1).SothedetectionofthesecondplanetwithMARVELSdataisnotexpected.ForHD68988'sunconrmedsecondplanet,itsestimatedperiodsaremuchlongerthantotalobservationtime.Sooneplanetplustrendandisexpectedtobeafavorablemodel. Theabovestated10ms)]TJ /F5 7.97 Tf 6.58 0 Td[(1detectionuncertaintyisonlyanestimation,sinceMARVELSdatapipelineisunderimprovement.Thevaluedependsonthealgorithms 146

PAGE 147

Figure5-18. EvidenceoffourplanetmodelofafourplanetsystemAra.Withthreevelocityreferences(24modelparameters),calculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamples.234RVobservationsover2768daysareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC.TheestimationfailstoconvergewiththecurrentMCMCsamples,NMCMC106. thedatapipelineuse,whichinvolvesvisibility,phase-velocityscale,etc.Unfortunately,intheattemptsofgeneratingMCMCsamplesfromtheposteriorPDFofHD17156andHD68988withoneplanetmodelbasedonMARVELSRVobservations,itisverydifculttoobtainaconvergedMCMCsequence.Andthesequencefailedtondamodenearthevaluesoftheknownplanetparameters.ItimpliesthattheRVmeasurementuncertaintycouldhavebeenunderestimatedorthereisothernoisesources. 147

PAGE 148

Figure5-19. FractionofMCMCsamplesincludedinthesub-domainforAra.ThefractionofMCMCsamplesincludedinthesub-domain,pVn,fortheplainMonteCarlointegration,plottedversusthetotalnumberofMCMCsamplesusedforcalculation.TheMCMCsamplesarebasedonafourplanetmodelwiththreevelocityreferencesforAra(25modelparameters,with234RVobservationsover2768days.TheestimationofpVnstartstoconvergeto0.013with106MCMCsamples 5.5BayesianModelSelectionforInstrumentAnalysis DuetothecomplicateddesignofMARVELS,includingthelargebeamsplittertoaccommodatemultiplebers,andtheberfeedingsystems,theremightbevarioussourcestocauseobservationuncertainties.Itisdifcultandcouldbebiasedtoestimatetheuncertaintiesonlyusingthepointestimation,especiallyforunknownnoisesources.TheBayesianmodelselectionmethodscanprovideapowerfultooltoanalyzetheinstrument.Forexample,withRVobservationsforreferencestaroreventhereferencelamp,wecangenerateMCMCsamplesfromthemodelswithandwithoutextranoise+ 148

PAGE 149

Figure5-20. RelativestandarddeviationoftheevidenceestimationoffourplanetmodelforAra.Afourplanetsystemwiththreevelocityreferences(24modelparameters),calculatedusingtheFractionEstimator,andplottedversusthetotalnumberofMCMCsamples.Solidline:thetotalrelativeestimationstandarddeviation.Dashline:relativeestimationstandarddeviationcausedbyGn.Dotted-dashline:relativeestimationstandarddeviationcausedbyFn.Thedashlineisoverlappedwiththesolidline,whichindicatedtherelativeestimationstandarddeviationismainlycausedbythatfromGn.234RVobservationsover2768daysareused.ThepriorPDFformodelparametersarelistedinTable 4-1 .AndNMCI=4NMCMC. 149

PAGE 150

orlineartrend,andcalculateevidencesforeachmodel.Thiscandetermineifwebelievethereisanextranoisesource,e.g.theinterferometricfringedrift,andtheuncertaintyoftheextranoise. WecanalsousetheBayesianmethodtoidentifythenoisesourcefromwhichpartoftheinstrument.Furthermore,inChapter2,itismentionedthatthefailureofPZTorthePZThostingcomputerwillcausetheobservationfromdifferentrunsofETnotlinkable.TheBayesianmethodscanalsobeusedtolinktheobservationsbysettinganewvelocityreferenceoncethephaselockingsystemfails. 5.6Summary Intheresearchdescribedinthisdissertation,Ipresentthreeprojects:1)InstrumentationofETandMARVELS,2)Modelselectionandparametercharacterizationforinteractingplanetarysystems,3)Bayesianmodelselectionbasedonamodelevidencecalculationalgorithm. ETandMARVELSaretheinstrumentsdesignedusingonDFDItechniquesforexoplanetdetection.IwasmainlyresponsiblefortheinstrumentationworkofETinitslatestage,andMARVELSforitspilotprogram.InChapter2,besidestheinstrumentcontroldesignandimplementation,Ifocusonthedevelopmentofthephaselockingsystem,whichisusedtostabilizetheinterferometerinbothinstrument. InChapter3,IexaminetheKeplerianmodeltforRVobservationfromasynthesissetofinteractingtwoplanetsystems,whicharegeneratedwithN-bodymodelsimulation.ThestudyimpliesthatitishelpfulforinteractingsystemtoconsiderN-bodymodelinsteadofonlytheKeplerianmodeliftheobservationprecisionisimproved.Withlimitedobservationuncertainty,theN-bodymodeltdoesnotdifferfromtheKeplerianmodelt.ButtheKepleriantstrendtoyieldbiasedestimationoftheorbitalparameters.ItismightbeusefultofurtherinvestigatethesystemsinMMRwiththeorderotherthan2:1. 150

PAGE 151

MCMCanalysisisaBayesianmethodtodeterminetheestimationuncertaintyofmodelparameters.Forexoplanetmodels,themodelevidenceisdifculttocalculateduetothehighdimensionintegrationproblem.Idevelopedanalgorithm,calledtheFractionEstimator,toestimatethemodelevidencebyusingafractionoftheMCMCsamplesinselectedsub-domain.Ithasbeenproventoworkwellformultipleplanetsystem.Andseveraltechnicalissuesareaddressed.However,theintegrationintheselectedsub-domainstillsufferfromthehighdimensionalityformodelswithmorethan4planets.Sothefutureworkcanberelatedwithhowtocalculatedthemulti-dimensionalintegrationinasub-domainnearthemode,possiblesolutioncouldbeusingtheLebesgueintegration,Voronoitessellation,etc.TheBayesianmodelselectionmethodscannotonlybeusedforexoplanetdetectionandcharacterization,butalsocanbeusedforanalyzingtheinstrument. 151

PAGE 152

REFERENCES [1] Bonavita,M.,Desidera,S.2007,A&A,468,721 [2] Bovy,J.,Hogg,D.W.,Roweis,S.T.2011,TheAnnalsofAppliedStatistics,5,2B,1657 [3] Chamber,J.E.1999,MNRAS,304,793 [4] Curiel,S.,Canto,J.,Georgiev,L.,Chavez,C.,Poveda,A.2010,A&A,525,78 [5] Formation,DetectionandDynamics.EditedbyRudolfDvrok,Wiley-VCH,2008 [6] Eggleton,P.P.,Kiseleva-Eggleton,L.2001,ApJ,562,1012 [7] Erskine,D.J.,Ge,J.2000,ASPConf.Ser.195,501 [8] Erskine,D.J.2003,PASP,115,255 [9] Feroz,F.,HobsonM.P.2008,MNRAS,384,449 [10] Fischer,D.A.,etal.2007,ApJ,669 [11] Fleming,S.W.,etal.2010,ApJ,718,1186 [12] Fleming,S.W.,etal.AJ,142,50 [13] Ford,E.B.2005,AJ,129,1706 [14] Ford,E.B.2006,ApJ,642,505 [15] Ford,E.B.,GregoryP.C.2007,StatisticalChallengesinModernAstronomyIV [16] Ge,J.2002,ApJ,571,165 [17] Ge,J.,Erskine,D.J.,Rushford,M.2002,PASP,114,1016 [18] Ge,J.,etal.2006,ApJ,648,683 [19] Ge,J.2006,AAS,20919102G [20] Gregory,P.C.2007,MNRAS,381,1607 [21] Guo,P.C.,Ford,E.B.2012,inprep [22] Haghighipour,N.2009,arXiv:0908.3328v1 [23] Hastings,W.K.1970,Biometrika,57,97 [24] Horne,J.H.,BaliunuasS.L.1986,ApJ,302,757 [25] Jeffreys,H.1961,TheTheoryofProbability(3rded.),Oxford,432 152

PAGE 153

[26] Kiseleva,L.G.,Eggleton,P.P.,Mikkloa,S.1998,MNRAS,300,292 [27] Laughlin,G.,Chambers,J.E.2001,ApJ,551,109 [28] Lee,M.H.,Peale,S.J.2002,ApJ,567,596 [29] Lee,B.L.,etal.2011,ApJ,728,32 [30] Loredo,T.J.,Berger,J.O.,Chernoff,D.F.,Clyde,M.A.,Liu,B.2012,StatisticalMethodology,9,1,101 [31] Lunine,J.I.,Macintosh,B.,Peale,S.2009,PhysicsToday [32] Mahadevan,S.,Ge,J.,Fleming,S.W.,Wan,X.,DeWitt,C.,vanEyken,J.C.,McDavitt,D.2008,PASP,120,1001 [33] Marchy,G.W.,Bulter,R.P.,Fischer,D.,Vogt,S.S.,Lissauer,J.J.,Rivera,E.J.2001,ApJ,556,296 [34] Marcy,G.W.,etal.2005,ApJ,619,570 [35] Marois,C.,etal.2008,Science,322,1348 [36] Metropolis,N.,Rosenbluth,A.W.,Rosenbluth,M.N.,Teller,A.H.,Teller,E.1953,JournalofChemicalPhysics,21,1087 [37] Newman,M.E.J.,Barkema,G.T.,1999,ClarendonPress [38] Payne,M.J.etal.2009,MNRAS,393,1219 [39] vanEyken,J.C.,Ge,J.,Mahadevan,S.,DeWitt,C.2004,ApJ,600,79 [40] vanEyken,J.C.,etal.2010,ApJS,189,156V [41] Wan,X.,Wang,J.,Ge,J.2010,Appl.Opt.49,5645 [42] Weinberg,M.D.2010,arXiv:0911.1777v2 [43] Wolszczan,A.,Frail,D.A.1992,Nature,355,145 [44] Wright,J.T.,Howard,A.W.,2009,ApJ,182,205 153

PAGE 154

BIOGRAPHICALSKETCH HavingsurvivedfromthebirthcontrolpolicyofP.R.China,Peng-Cheng,GuowassuccessfullyborninthecityofJinan(N36.67,E117.00),thecapitalofShandongprovinceinChina.Hegrewupinanartistfamily,whichcultivatedhisfondofartdesign.Hebecameinterestedinscienceandtechnologyduringhisearlyschoollife.AftergraduatingfromPekingUniversitywithaBachelordegreeinAstrophysicsfollowedbythreeyearsworkasanelectronicengineer,hebeganhisAstronomyPhDlifeintheUnitedStates,duringwhicheverythinghadchangedbeyondtheplan.HisworkandresearchfocusturnedfromastronomicalinstrumentationofexoplanetdetectiontothemodelanalysisandBayesianstatisticsonexoplanetarysystemswithachangeoftheresearchtopic.HereceivedthedegreeofMasterofEngineeringinElectronicEngineeringfromUniversityofFloridainAugust,2011,andreceivedhisPhDdegreeinAstronomyfromUniversityofFloridainDecember,2012.Hebecameoneofthefoundersofaninformationtechnologystartupcompanyinhishometown,focusingonthetechniquesofdatamining.Heisafanofmulti-taskandexperienceofdifferentactivities,suchasstockinvestment,parachuting,interiordesign,andjewelrymaking. 154