PEXO v2.0
Fabo Feng & Maksym Lisogorskyi
16 September, 2020
1 Introduction
PEXO is a package for making precise exoplanetology. As compared with previous models and packages, PEXO is significantly advanced and accounts for the orbital dynamics of binary motion and stellar reflex motions induced by planetary companions. PEXO treats both classic and relativistic effects such as the Roemer, Shapiro, and Einstein time delays both in the Solar System and in the target system
PEXO is able to model timing to a precision of 1 ns, astrometry to a precision of 1 microarcsecond, and radial velocity to a theoretical precision of 1 \(\mu\)m/s and a realistic precision of 1 cm/s. PEXO was bechmarked with the pulsar timing package TEMPO2. Theoretical and computational details of the code are described in the paper by Feng et al (2019).
2 Installation
The code is written in R and depends on several libraries. To install R on Linux, download it here or, in Ubuntu:
sudo apt-get install r-basein MacOS:
brew install rClone this repository:
git clone https://github.com/phillippro/pexo.gitInstall missing R libraries “cd pexo” then
Rscript install_dependencies.R(might require su privileges to be installed): or just install them manually using install.packages(‘package_name’) in R console.
3 Usage
3.1 Command Line
To use PEXO, one needs to go to the directory pexo/code/ and run command lines such as
Rscript pexo.R -m emulate -p 'GJ551' -t '2440000 2450000 10' -i 'APF' -m 'TAR'which will run a simulation of the astrometry and radial velocity variation of GJ551 from JD[UTC]=2440000 to JD[UTC]=2450000 with a time step of 10 days.
By runing
Rscript pexo.R -m fit -p HD10700 -N 100 -d '../input/HD10700'The code will fit a single star model to the data saved in the ‘../input/HD10700’ directory.
One can also simulate binary systems such as alpha Cen by using
Rscript pexo.R -m emulate -p 'HD128620' -t '2440000 2450000 10' -i 'ESO' -o '../input/HD128620.par' -C 1One can fit a binary model to various timing, RV and astrometry data using
Rscript pexo.R -m fit -p 'HD128620' -N 100 -C 1 -d '../input/HD128620'3.2 Input parameters
The revelant command line arguments are listed as follows.
| Short name | Full name | Meaning |
|---|---|---|
| -m | --mode | PEXO mode: emulate or fit [optional; default=emulate] |
| -c | --component | PEXO model component: timing (T), astrometry (A), radial velocity (R) and their combinations [optional; default=TAR] |
| -i | --ins | Instrument or observatory [mandatory for emulation mode; default=NA] |
| -P | --par | Parameter file: parameters for astrometry and observatory [default: automatically obtained from simbad if not find parameter file] |
| -N | --Niter | Sample size of MCMC [optional; default=1000] |
| -C | --Companion | Companion number [optional; default=0] |
| -g | --geometry | geometric orbit or relativistic orbit [optional; default=TRUE] |
| -n | --ncore | Number of cores [optional; default=4] |
| -c | --component | PEXO model component: timing (T), astrometry (A), radial velocity (R) and their combinations [optional; default=TR] |
| -t | --time | Two options are possible. 1. Timing file: epochs or times could be in 1-part or 2-part JD[UTC] format; 2. Format of "Start End By" [mandatory if mode=emulate] |
| -p | --primary | primary star name [mandatory] |
| -s | --secondary | secondary star name [default=NA] |
| -M | --mass | Mass of primary in unit of solar mass [optional; default=1] |
| -d | --data | Data directory: directory with timing, RV or astrometry data files [default= ‘../input/primary’ where primary is the star name] |
| -v | --var | Output variables [optional; default=NULL] |
| -o | --out | Output file name: relative or absolute path [optional; default=out.txt] |
| -f | --figure | Output figure and verbose: FALSE or TRUE [optional; default= TRUE] |
| -V | --verbose | Verbose: FALSE or TRUE [optional; default= FALSE] |
Since the astrometry and radial velocity modeling depends on the the output of timing model, T should always be included in the -c or --component argument.
3.3 Input timing file (only for emulate mode)
For emulate mode, the ‘-t’ or ‘-time’ argument is mendatory. It could either be a timing file or a string with "Start End By" format.
The timing file could be two-part or one-part JD or MJD (MJD=JD-2400000.5) in UTC time standard. The former can store epochs with precision of \(10^{-14}\) second while the the latter can store epoch with precision of \(10^{-6}\) second or microsecond in a 64-bit computer.
The "Start End By" format timing argument is composed of the start epoch (Start), the end epoch (End) and the time step (By). For example, a run of PEXO with -t "2456640.5 2458462.5 0.5" will simulate the system from JD2456640.5 to JD2458462.5 by a time step of 0.5 days. The -t argument could also be in MJD format such as -t "56640 58462 0.5". The times generated from the sequence would be transformed into 2-part JD format for high precision emulation.
3.4 orbital parameter file
The orbital parameter file is specified by assigning a file path to -o or –orbit. By default, pexo will look for orbital parameters in ../input/primary.par (primary is the value of -p or –primary). A typically parameter file is as follows:
secondary HD128621
mass 1.1
logmC -0.09761283
logP 10.281
e 0.5179
I 79.32
omegaT 232
Omega 205
Tp 2435328.96It can redefine the secondary name if --secondary is not given through commandline. It will give mass of the primary and the initial/optical values of logarithmic mass of companion (logmC), logarithmic orbital period in days (logP), eccentricity (e), inclination (I in degree), argument of periastron for the primarty (omegaT in degree), longitude of ascending node (Omega in degree), and the priastron epoch (Tp in Julian days). In addition to these orbital parameters which is only for binaries, extra parameters of astrometry and observatory should be given for both singles and binaries. These additional parameters could either be automatically found by PEXO getting astrometry from Simbad (using get_cor_from_simbad.py) or manually given by specifying the parameter file (-P or --parfile). If specifying these parameters manually, the user should create a file with the following format. I use instrument APF/LICK and alpha Cen A for an example:
phi.APF 37.3425
elong.APF -121.63825
height.APF 1.274
epoch 2457206.375
ra 219.902058333333
dec -59.1660075
pmra -3679.25
pmdec 473.67
plx 743
rv -21.4
logmC -0.09761283 -10 10 U
logP 10.281 -10 20 U
e 0.5179 0 1 U
I 79.32 0 180 U
omegaT 232 0 360 U
Omega 205 0 360 U
Tp 2435328.96 2335328.96 2535328.96 U The first three parameters are the latitude, longitude and height of the instrument or observatory. For fitting mode with multiple data sets, coordinates for multiple instruments are needed. The format to specify an instrument is using xxx.yyy where xxx is the coordinate name and yyy is the instrument name. Then ra, dec, pmra, pmdec, plx, and rv are the astrometric parameters from astrometric surveys such as Gaia and Hipparcos. The epoch parameter is the reference astrometry epoch for a given survey. The other parameteters define the priors for orbital parameters and are only necessary for fitting mode for stars with companions (i.e. –Companion>0) For example, the prior for e is uniform from 0 to 1. The initial values are taken from the orbital parameters given in ‘–orbfile’ and are for initializing MCMC.
The other model related parameters are given in ../input/basic.par. A summary definition of parameters are as follows. The bold-faced values are default ones.
| parameter | unit | options or examples | meaning |
|---|---|---|---|
| RefType | - | none, refro, refco, refcoq | computation method for atmospheric refraction |
| EopType | - | 2006, 2000B | type of Earth rotation model and corresponding Earth orientation parameters |
| TaiType | - | instant, scale | UTC to TAI method |
| TtType | - | BIPM, TAI | TAI to TT method |
| unit | - | TCB, TDB | output quantities compatible with TCB or TDB time standard |
| DE | - | 430, 430t, 438, … | JPL ephemerides |
| TtTdbMethod | - | eph, FB01, FBgeo | TT to TDB method |
| SBscaling | - | FALSE, TRUE | linear scaling between tB and tS due to relativistic effects |
| PlanetShapiro | - | TRUE, FALSE | planetary shapiro delay |
| CompareT2 | - | FALSE, TRUE | calculate uSB using TEMPO2 method for comparison |
| RVmethod | - | analytical, numerical | the method used for RV modeling, numerical is used only for comparison |
| LenRVmethod | - | T2, PEXO | the method used to derive RV lensing, T2 is used by default to be consistent with shapiro delay model in PEXO |
| BinaryModel | - | none, DDGR, kepler | binary model |
| ellipsoid | - | WGS84, GRS80, WGS72 | ellipsoidal (normal) Earth Gravitational Model |
| epoch | JD or MJD | 2448349.06250 | epoch when the astrometry and position of the target is measured |
| observatory | - | CTIO | observatory name |
| xtel | metre | 1814985.3 | geocentric position of the telescope in the International Terrestrial Reference Frame (ITRF) |
| ytel | metre | -5213916.8 | geocentric position of the telescope in ITRF |
| ztel | metre | -3187738.1 | geocentric position of the telescope in ITRF |
| tdk | K | 278 | ambient temperature at the observer |
| pmb | millibar | 1013.25 | pressure at the telescope |
| rh | - | 0.1 | relative humidity at the observer (range 0-1) |
| wl | \(\mu\)m | 0.5 | effective wavelength of the source |
| tlr | K/metre | 0.0065, any value>0 | Temperature lapse rate in the troposphere |
| g | 1, 0, any other values >0 | one of the PPN parameters | |
| mT | \(M_\odot\) | 1.1055 | target mass |
| mC | \(M_\odot\) | 0.9373 | companion mass |
| ra | degree | 219.9175253 | right ascension (RA) of the barycenter (TSB) |
| dec | degree | -60.8371344 | declination (DEC) of the barycenter (TSB) |
| plx | mas | 747.1700008 | parallax of the barycenter (TSB) |
| pmra | mas/yr | -3649.4980522 | proper motion in RA of the barycenter (TSB) |
| pmdec | mas/yr | 624.7691720 | proper motion in DEC of the barycenter (TSB) |
| rv | km/s | -22.3929553 | radial velocity of the barycenter (TSB) |
| aT | au | 10.80332 | semi-major axis of the barycentric motion of the target |
| P | year | 79.929 | orbital period of the target |
| e | - | 0.5208 | eccentricity |
| I | degree | 79.32 | inclination |
| omegaT | degree | 52.006 | argument of periastron |
| Omega | degree | 205.064 | longitude of ascending node |
| Tp | JD or MJD | 2435328.96 | periastron epoch |
In the above table, epoch is the time when the position and astrometry of the target sysstem is measured. It is tpos defined in the PEXO paper.
3.5 Observatory data
PEXO will first look for observatory data by finding xtel, ytel, and ztel from the parameter file. If it does not find these parameters, it will look for elong (longitude in degree), phi (latitude in degree) and height (altitude in km). If these parameters are not given, PEXO will look for the observatory name (observatory) and code (ObsCode). It will look for the observatory data in observatories/observatory_MPC.csv or in observatories/satellite_list.csv. For space-based observatory, the atmospheric refraction and delay are zero and would not be implemented by PEXO. For ground-based telescope, the RefType parameter should be refro (recommended), refco, or refcoq for the calculation of refraction in astrometry modeling. If RefType is none, the atmospheric refraction is zero. The tropospheric delay and its time derivative are automatically implemented for ground-based observatories and thus do not depend on the choice of RefType.
For space-based observatory, “code/GetSpaceObsEph.py” will be used to find the telescope’s ephemerides. So you need to install python as well as python packages sys, astropy, numpy, and astroquery.
3.6 PEXO output
3.6.1 Emulate mode
A diagram for propagation of the light ray from the target star to the observer is shown below to aid the understanding of the ouptut quantities.
Fig. 1. Illustration of a light ray emitted from the target (T) and observed by the observer (O). The target system is composed of the target (T) and its companion (C). The binary barycenter is denoted by B. The observer is in the solar system with a barycenter at S.
There are outputs from four functions in PEXO:
OutObs <- time_Utc2tb(utc,Par)utc is the input 2-part JD epochs, Par is the input and derived parametes. The output OutObs is a list of variables related to the transformation from JD[UTC] to JD[TCB] or JD[TDB].
OutTime <- time_Ta2te(OutObs,Par)This function uses OutBary and Par to transform JD[TCB] to BJD[TCB] to light emission time. Thus OutBary and OutTime are the outputs from the timing models.
OutAstroT <- astro_FullModel(OutObs,OutTime,Par,Mlens=Par$mC,component='T')
OutAstroC <- astro_FullModel(OutObs,OutTime,Par,Mlens=Par$mT,component='C')OutAstroT and OutAstroC are outputs of the astrometry modeling of the T and C component in the target system. Since the astrometry function astro_FullModel calls OutBary and OutTime, these astrometry outputs depend on outputs of timing model.
OutRv <- rv_FullModel(OutObs,OutTime,Par)OutRv is the output of radial velocity modeling and also depends on the outputs of timing model.
These output lists will be combined as OutAll to be saved as ascii file if the output variables -v are specified in the command line.
We list the output variables in OutAll, their unit and meaning in the following table.
| variable | unit | meaning |
|---|---|---|
| AbeTarget | second | target aberration delay |
| BJDtcb | day | BJD[TCB] |
| BJDtdb | day | BJD[TDB] |
| BT | au; au/yr | Position and velocity vectors from TSB to T |
| DefEarth | rad | Deflection vector due to Earth lensing |
| DefJupiter | rad | Deflection vector due to Jupiter lensing |
| DefMars | rad | Deflection vector due to Mars lensing |
| DefMercury | rad | Deflection vector due to Mercury lensing |
| DefMoon | rad | Deflection vector due to Moon lensing |
| DefNeptune | rad | Deflection vector due to Neptune lensing |
| DefSaturn | rad | Deflection vector due to Saturn lensing |
| DefSun | rad | Deflection vector due to Sun lensing |
| DefUranus | rad | Deflection vector due to Uranus lensing |
| DefVenus | rad | Deflection vector due to Venus lensing |
| delevation | rad/day | time derivative of elevation angle |
| delevationT2 | rad/day | time derivative of elevation angle computed by TEMPO2 method |
| DirObs | rad | observed right ascension and declination of the target |
| dl.all | rad | Light deflection vector due to all effects |
| dl.woRef | rad | Light deflection vector due to all effects except for atmospheric refraction |
| dTCB.dTT | - | dTCB/dTT |
| dTDB.dTT | - | dTDB/dTT |
| dzenith | rad/day | Time derivative of zenith: dzenith/dt |
| EinsteinIS | second | Einstein delay due to relative motion between TSB and SSB |
| EinsteinTarget | second | Einstein delay in the target system |
| elevation | rad | elevation angle |
| elevationT2 | rad | elevation angle calculated using TEMPO2 method |
| emrat | - | Earth-Moon mass ratio |
| Eph | - | a list of ephermerides of solar system objects |
| EphEarth | km; km/s | Earth ephemeris in the Barycentric celestial reference system (BCRS) frame; units are denoted by columns names |
| EphJupiter | km; km/s | Jupiter ephemeris in BCRS |
| EphMars | km; km/s | Mars ephemeris in BCRS |
| EphMercury | km; km/s | Mercury ephemeris in BCRS |
| EphMoon | km; km/s | Moon ephemeris in BCRS |
| EphNeptune | km; km/s | Neptune ephemeris in BCRS |
| EphSaturn | km; km/s | Saturn ephemeris in BCRS |
| EphSun | km; km/s | Sun ephemeris in BCRS |
| EphUranus | km; km/s | Uranus ephemeris in BCRS |
| EphVenus | km; km/s | Venus ephemeris in BCRS |
| GM | km; km/s | Position and velocity vector from the geocenter to the Moon |
| GO | km; km/s | Position and velocity vector from the geocenter to the observer/telescope |
| JDtai | JD | JD[TAI] or TAI |
| JDtcb | JD | JD[TCB] or TCB |
| JDtcg | JD | JD[TCG] or TCG |
| JDtdb | JD | JD[TDB] or TDB |
| JDtt | JD | JD[TT] or TT |
| JDut1 | JD | JD[UT1] or UT1 |
| leap | second | leap second |
| li | - | unit vector or direction of the incident or pre-refraction light ray |
| limll | - | li - ll |
| ll | - | direction of the light ray after leaving the target system |
| llmle | - | ll - le |
| lo | - | direction of the light ray at the telescope before being observed |
| lomli | - | lo-li |
| MO | km; km/s | Position and velocity vector from the Moon to the observer |
| OffAbe | arcsecond | offset due to aberration in (dRA*, dDEC) |
| OffAbe1 | arcsecond | offset due to first-order aberration |
| OffAbe2 | arcsecond | offset due to second-order aberration |
| OffAll | arcsecond | offset due to all effects |
| OffLenS | arcsecond | offset due to all lensing in the solar system |
| OffLenT | arcsecond | offset due to all lensing in target system |
| OffRef | arcsecond | offset due to atmospheric refraction |
| OL | - | a list of observer to solar system body (lens) vectors |
| OutBT | - | a list of outputs from binary models |
| rBT | au | position vector from TSB to T |
| RBT | au | length of rBT |
| ref | rad | refraction vector |
| Ref | rad | refraction angle |
| rOB | pc | position vector from the observer to the TSB |
| rOC | pc | position vector from the observer to the companion (C) |
| Roemer1 | second | first order Roemer delay in the solar system |
| Roemer2 | second | second order Roemer delay in the solar system |
| Roemer3 | second | third order Roemer delay in the solar system |
| RoemerOrder | second | a combined list of Roemer1, Roemer2 and Roemer3 |
| RoemerSB | second | Roemer delay using SB rather than ST as the reference direction (only for comparison) |
| RoemerSolar | second | total Roemer delay in the solar system |
| RoemerT2 | second | Roemer delay calculated using the TEMPO2 method (including the total effects and effects for different terms) |
| RoemerTarget | second | Roemer delay in the target system |
| rOT | pc | position vector from the observer to the target |
| rSB | pc | position vector from SSB to TSB |
| rSC | pc | position vector from SSB to the companion |
| rST | pc | position vector from SSB to the target |
| rTC | au | position vector from the target to the companion |
| RvBT | m/s | radial velocity for TSB to T |
| RvGO | m/s | radial velocity for geocenter to observer |
| RvgsO | m/s | general and special relativistic effect on RV at the observatory or in the solar system |
| RvgT | m/s | general and special relativistic effect on RV in the target system |
| RvlO | m/s | lensing RV in the solar system |
| RvLocal | m/s | all RV effects in the solar system |
| RvlT | m/s | lensing RV in the target system |
| RvRemote | m/s | all RV effects in the target system |
| RvSB | m/s | RV due to motion of TSB w.r.t. SSB |
| RvSG | m/s | RV due to motion of the geocenter w.r.t. SSB |
| RvSO | m/s | RV due to motion of the observer w.r.t. SSB |
| RvsT | m/s | special relativitistic effect on RV in the target system |
| RvST | m/s | RV due to motion of the target w.r.t. SSB |
| RvTot | m/s | total RV |
| RvTropo | m/s | tropospheric RV |
| SB | pc; au/yr | position and velocity vectors from the SSB to TSB |
| SG | km; km/s | position and velocity vectors from the SSB to the geocenter |
| ShapiroEarth | second | Shapiro delay due to Earth |
| ShapiroJupiter | second | Shapiro delay due to Jupiter |
| ShapiroMars | second | Shapiro delay due to Mars |
| ShapiroMercury | second | Shapiro delay due to Mercury |
| ShapiroMoon | second | Shapiro delay due to Moon |
| ShapiroNeptune | second | Shapiro delay due to Neptune |
| ShapiroPlanet | second | a combined list of Shapiro delays due to solar system objects |
| ShapiroSaturn | second | Shapiro delay due to Saturn |
| ShapiroSolar | second | Shapiro delay in the solar system |
| ShapiroSun | second | Shapiro delay due to Sun |
| ShapiroTarget | second | Shapiro delay in the target system |
| ShapiroUranus | second | Shapiro delay due to Uranus |
| ShapiroVenus | second | Shapiro delay due to Venus |
| SO | km; km/s | position and velocity vectors from the SSB to the observer |
| SolarDef | rad | deflection angle (vector) due to lensing in the solar system |
| SolarDefList | rad | a list of deflection angles due to lensing by solar system objects |
| TargetDelay | second | total delay in the target system |
| tauE | JD | proper emission time |
| tB | JD | coordinate light arrival time at TSB |
| TDBmTTgeo | second | TDB-TT at the geocenter |
| TropoDelay | second | tropospheric delay |
| TropoDelayT2 | second | tropospheric delay calculated using uSB(t=tpos) or ub (see paper) as the reference direction as done in TEMPO2t |
| tS | JD | same as BJD[TCB]; coordinate ligth arrival time at SSB |
| U | rad | eccentric anomaly |
| uBT | - | unit vector for rBT |
| uo | - | observed direction of the target |
| uOB | - | unit vector for rOB |
| uOC | - | unit vector for rOC |
| uommlo | - | uo+lo |
| uommlo1 | - | uo+lo1 |
| uommlo2 | - | uo+lo2 |
| uommlo3 | - | uo+lo3 |
| uOT | - | unit vector for rOT |
| uSB | - | unit vector for rSB |
| uSB.T2 | - | uSB calculated using the TEMPO2 method (ignoring third order effects) |
| uST | - | unit vector for rST |
| VacuumIS | - | vacuum delay in interstellar medium |
| vBT | au/yr | velocity of T w.r.t. TSB |
| vGO | km/s | velocity of the observer w.r.t. the geocenter |
| vOB | au/yr | velocity of TSB to the observer |
| vOT | au/yr | velocity of target to the observer |
| vSB | au/yr | velocity of TSB to SSB |
| vST | au/yr | velocity of the target to SSB |
| xp | rad | parameter for polar motion of the Earth |
| yp | rad | parameter for polar motion of the Earth |
| ZB | - | barycentric correction of Doppler shift |
| ZBwe | - | barycentric correction of Doppler shift using Wright & Eastman 2014 method |
| Zcomb | - | combined list of all doppler shifts |
| ZenIn | rad | zenith angle |
| ZenInT2 | rad | zenith angle using uSB(t=tpos) or ub (see paper) |
| zenith | rad | zenith vector |
| ZgO | - | doppler shift due to general relativistic effect in the solar system |
| ZgsO | - | doppler shift due to relativistic effects in the solar system |
| ZgsO.de | - | doppler shift due to relativistic effects in the solar system calculated using JPL ephemerides |
| ZgSS | - | combined list of gravitational doppler shifts due to solar system objects |
| ZgsT | - | doppler shift due to relativistic effects in the target system |
| ZgTk | - | doppler shift due to general relativistic effect in the target system |
| ZkpO | - | doppler shift due to parallax delay in the solar system |
| ZkpT | - | doppler shift due to parallax delay in the target system |
| Zlensing | - | combined list of doppler shifts due to lensing by solar system objects |
| ZlO | - | doppler shift due to solar system lensing |
| Zlocal | - | local doppler shift |
| ZlT | - | doppler shift due to target system lensing |
| Zremote | - | local doppler shift |
| ZsO | - | special relativistic doppler shift in the solar system |
| ZSO | - | doppler shift due to the motion of SSB w.r.t. the observer |
| ZsT | - | special relativistic doppler shift in the target system |
| ZST | - | doppler shift due to the motion of target w.r.t. SSB |
| ZST0 | - | doppler shift due to the motion of target w.r.t. SSB using uSB.T |
| ZST0 | - | doppler shift due to the motion of target w.r.t. SSB using uSB.T |
| zTDBmTTgeo | - | doppler shift corresponding to the time derivative of TDB-TT at the geocenter |
| zTDBmTTobs | - | doppler shift corresponding to the time derivative of the observer term in TDB-TT |
| zTDBmTTobsR | - | zTDBmTTobs due to rGO |
| zTDBmTTobsV | - | zTDBmTTobs due to vGO |
| Ztot | - | total doppler shift |
| Ztropo | - | tropospheric doppler shift |
4 Examples
4.1 Simulation of single star systems
PEXO can be used to simulate the timing, RV and astrometry variation of a single star. PEXO will automatically find astrometry data from Simbad if it is not given manualy. To explain the output of PEXO, tau Ceti is simulated using the following commandline:
Rscript pexo.R -p HD10700 -t '2450000 2460000 10' -i 'ESO' -c 'TAR'In the output of PEXO, the various timing, RV and astrometric effects are printed. Summary plots for these effects are saved as pdfs with names of ../results/AllTimes_HD10700_timing_kepler_Ntime1001.pdf, ../results/absolute_HD10700_astrometry_kepler_Ntime1001_none.pdf, and ../results/paper_HD10700_RV_kepler_Ntime1001_none.pdf. The user can also choose which variables to be saved in an output file with default name of out_pexo.txt. The user can specify -o and -v to define the output file name and variable names, respectively. For example, if the user is interested in barycentric correction of a star, he/she could add -v 'BJDtdb ZB' to the commandline and then get BJDtdb and ZB from the output file. The barycentric RV is ZB*c where c is the speed of light.
4.2 Simulation of stars with companions
The main advantage of PEXO is to model the timing, RV and astrometry for binaries to a high precision. alpha Cen A and B are a good exmample to demonstrate this. For exmaple, the following commandline
Rscript pexo.R -p HD128620 -t '2450000 2460000 10' -i 'PFS' -c 'TAR' -C 1will simulate alpha Cen A and B from JD2450000 to JD2460000 by a time step of 10 days. The instrument is PFS and the components in the simulation include timing (T), astrometry (A), and RV (R). As for the case of single stars, one will get summary plots and magnitudes of various effects as well as output ascii file (named out_pexo.txt by default).
4.3 Fitting for single star systems
For emulation/simulation mode, one cannot optimize the astrometric and orbital parameters. By fitting these parameters to various data, one is able to find the optimal astrometric and orbital solution. To do that, the data files for fitting should be named according to the following format if HD10700 is the primary star and xxx is an instrument name:
Timing data (typically is data for time delay): HD10700_xxx.delay
RV data: HD10700_xxx.rv
Absolute astrometry data: HD10700_xxx.abs
Relative astrometry data: HD10700_xxx.relFor timing data, there are at least two columns: JD[UTC], delay [s], (delay error [s]). For RV data, there are three columns: JD[UTC], RV [m/s], RV error [m/s]. Note that the time and RV in RV data should not be barycentrically corrected. For absolute and relative astrometry data, there are at least five columns: JD[UTC], ra [deg], ra error [mas], dec [deg], dec error [mas], (lambda [\(\mu\)m]). Lambda is the wavelength of an ground-based astrometric measurement. If not given, PEXO will assumes 0.55 \(\mu\)m wavelength which will be used to determine the amount of atmospheric refraction.
Proxima Cen or GJ551 is taken as an example to explain how to use the PEXO fitting mode. For example, the following commandline
Rscript pexo.R -m fit -p GJ551 -N 100or you can specify the output Robj name, e.g.
Rscript pexo.R -m fit -p GJ551 -N 100 -o 'test.Robj'will launch 4 MCMC and each has a length of 100. The pexo outputs (exclude the variation of barycentrically corrected RV files) are explained as follows:
| file name | meaning |
|---|---|
| ../results/GJ551_barycorrected.rv | Barycentrically corrected combined RV |
| ../results/GJ551_0companion_llmax-2810322_N100_einsteinTRUE.pdf | Summary plots for the model prediction (red) and data (black) as well as residuals and mcmc trace plots etc. |
| ../results/GJ551_0companion_llmax-2810322_N100_einsteinTRUE_fancy.pdf | Fancy plots for PEXO-simulated optimal fit to raw data as well as small panels showing residuals |
| ../results/GJ551_0companion_llmax-2810322_N100_einsteinTRUE.Robj | All outputs of PEXO are saved as an R object (xxx.Robj). The path of this file can be defined by the -o input parameter. |
In the Robj file, the MCMC related variables are explained as follows.
| variable name | meaning |
|---|---|
| Data | A data frame saving all input data, first column: UTC time; second column: RV[m/s], or delay [second], or RA [deg], or dRA [mas]; third column: error of the second column; fourth column: NA for RV and delay data and DEC[deg] or dDEC[mas] for absolute and relative astrometry respectively; fifth column: error of fourth column; sixth column: star name; seventh column: data type; eighth column: instrument; ninth column: wavelength[\(\mu\)m] for astrometry data |
| model | Model prediction at the raw data epochs, the columns have the same meanings. |
| pred | Model prediction at simulated epochs (typically with higher cadence and with constant time steps sampling over the whole data timespan), the columns have the same meanings. |
| ParStat | Statistics of parameters infered from MCMC posterior samples. Columns are model parameters and rows are the statistics for them, including xopt: optimal value at the maximum a posteriori, x1per: 1% quantile; x99per: 99% quantile; x10per: 10% quantile; x90per: 90% quantile; xminus.1sig: xopt-1sigma; xplus.1sig: xopt+1sigma; and other statistics are easy to understand. The parameters of raOff [mas], decOff [mas], pmraOff [mas/yr], pmdecOff [mas/yr], plxOff [mas] are offsets added to the initial astrometry of the barycenter of the target system. The parameters of logjitterRv.star.instrument and bRv.star.instrument are log jitter and intercept for a given instrument and star. For binary or star-companion system, the orbital parameters are also given with the number of companions added to the end of each parameter name. In specific, logmC1 is the log mass of the first companion in units of solar mass; logP1 is the log orbital period in units of day; e1 is eccentricity; I1 is inclination in radiant; omegaT1 is the primary star argument of periastron in radiant; Omega1 is the longitude of ascending node in units of radiant; Mo1 is the mean anomaly in units of radiant. |
4.4 Fitting for stars with companions
To fit a binary model to the data, you can run the following commandline (e.g., Kruger 60 or HD239960)
Rscript pexo.R -m fit -p HD239960 -N 100 -C 1Then you will get similar output as for the case of single stars.
5 Run PEXO using shell script
A simple shell wrapper is included for convenience. The arguments are identical to pexo.R, type ./pexo.sh --help so see help.
To run PEXO via this script, you need to set an environment variable $PEXODIR to a path to the PEXO repository. It is also recommended to create an alias for this script to run it from anywhere in the terminal. To do that, add
export PEXODIR=/example/path/to/pexo
alias pexo="/example/path/to/pexo/pexo.sh"to your ~/.bashrc or ~/.bash_profile if you’re using bash, or
setenv PEXODIR /example/path/to/pexo
alias pexo /example/path/to/pexo/pexo.shto ~/.tcshrc if you’re using tcsh. You’ll need to open a new terminal or type source ~/.bashrc / source ~/.tcshrc to apply this. You may also source ~/.bash_profile to permanently save the settings. Note that the path to pexo should be absolute. Then you can open a new terminal window to use PEXO everywhere in your computer.
For example, you have the following two options to run a simulation of Tau Ceti.
You may either go to the pexo/code directory and run
Rscript pexo.R -p HD128620 -t '2450000 2460000 10' -i 'PFS' -c 'TAR' -C 1or run the following shell script in your current directory by providing paths for timing and parameter files
pexo -p HD128620 -t '2450000 2460000 10' -i 'PFS' -c 'TAR' -C 1By running the above command line, PEXO will look for the parameter and data files in the ../input/ directory and save output file if you specify the -v and -o arguments.
6 Future development
PEXO v2 is able to emulate a system or fit a Keplerian model to timing, astrometry and RV data sets. A python wrapper are also developed for python users. PEXO v3 will provide functionality for pulsar timing and gravitational wave searches. Feedback from PEXO users and contribution from the astronomical community are appreciated and are important to improve the software.