Fit Proper Motion Anomaly
Octofitter.jl supports fitting orbit models to astrometric motion in the form of GAIA-Hipparcos proper motion anomaly (HGCA; https://arxiv.org/abs/2105.11662). These data points are calculated by finding the difference between a long term proper motion of a star between the Hipparcos and GAIA catalogs, and their proper motion calculated within the windows of each catalog. This gives four data points that can constrain the dynamical mass & orbits of planetary companions (assuming we subtract out the net trend).
If your star of interest is in the HGCA, all you need is it's GAIA DR3 ID number. You can find this number by searching for your target on SIMBAD.
For this tutorial, we will examine the star and companion HD 91312 A & B discovered by SCExAO. We will use their published astrometry and proper motion anomaly extracted from the HGCA.
We will also perform a model comparison: we will fit the same model to four different subsets of data to see how each dataset are impacting the final constraints. This is an important consistency check, especially with proper motion / absolute astrometry data which be susceptible to systematic errors.
The first step is to find the GAIA source ID for your object. For HD 91312, SIMBAD tells us the GAIA DR3 ID is 756291174721509376.
Fitting Astrometric Motion Only
Initial setup:
using Octofitter, Distributions, RandomWe begin by finding orbits that are consistent with the astrometric motion. Later, we will add in relative astrometry to the fit from direct imaging to further constrain the planet's orbit and mass.
Compared to previous tutorials, we will now have to add a few additional variables to our model. The first is a prior on the mass of the companion, called mass. The units used on this variable are Jupiter masses, in contrast to M, the primary's mass, in solar masses. A reasonable uninformative prior for mass is Uniform(0,1000) or LogUniform(1,1000) depending on the situation.
For this model, we also want to place a prior on the host star mass rather than system total mass. For exoplanets there is litte difference between these two values, but in this example we have a reasonably informative prior on the host mass, and know from the paper that the companion is has a non-neglible effect on the total system mass.
To make this parameterization change, we specify priors on both masses in the @system block, and connect it to the planet.
Planet Model
@planet b Visual{KepOrbit} begin
a ~ LogUniform(0.1,20)
e ~ Uniform(0,0.999)
ω ~ UniformCircular()
i ~ Sine() # The Sine() distribution is defined by Octofitter
Ω ~ UniformCircular()
mass = system.M_sec
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,57423.0) # epoch of GAIA measurement
endPlanet b
Derived:
ω, Ω, mass, θ, tp,
Priors:
a Distributions.LogUniform{Float64}(a=0.1, b=20.0)
e Distributions.Uniform{Float64}(a=0.0, b=0.999)
ωy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
ωx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
i Sine()
Ωy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Ωx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
θy Distributions.Normal{Float64}(μ=0.0, σ=1.0)
θx Distributions.Normal{Float64}(μ=0.0, σ=1.0)
Octofitter.UnitLengthPrior{:ωy, :ωx}: √(ωy^2+ωx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:Ωy, :Ωx}: √(Ωy^2+Ωx^2) ~ LogNormal(log(1), 0.02)
Octofitter.UnitLengthPrior{:θy, :θx}: √(θy^2+θx^2) ~ LogNormal(log(1), 0.02)
System Model & Specifying Proper Motion Anomaly
Now that we have our planet model, we create a system model to contain it.
We specify priors on plx as usual, but here we use the gaia_plx helper function to read the parallax and uncertainty directly from the HGCA catalog using its source ID.
We also add parameters for the star's long term proper motion. This is usually close to the long term trend between the Hipparcos and GAIA measurements. If you're not sure what to use here, try Normal(0, 1000); that is, assume a long-term proper motion of 0 +- 1000 milliarcseconds / year.
@system HD91312_pma begin
M_pri ~ truncated(Normal(1.61, 0.1), lower=0) # Msol
M_sec ~ LogUniform(0.5, 1000) # MJup
M = system.M_pri + system.M_sec*Octofitter.mjup2msol # Msol
plx ~ gaia_plx(gaia_id=756291174721509376)
# Priors on the center of mass proper motion
pmra ~ Normal(-137, 10)
pmdec ~ Normal(2, 10)
end HGCALikelihood(gaia_id=756291174721509376) b
model_pma = Octofitter.LogDensityModel(HD91312_pma)LogDensityModel for System HD91312_pma of dimension 14 with fields .ℓπcallback and .∇ℓπcallback
After the priors, we add the proper motion anomaly measurements from the HGCA. If this is your first time running this code, you will be prompted to automatically download and cache the catalog which may take around 30 seconds.
Sampling from the posterior (PMA only)
Because proper motion anomaly data is quite sparse, it can often produce multi-modal posteriors. If your orbit already has several relative astrometry or RV data points, this is less of an issue. But in many cases it is recommended to use the Pigeons.jl sampler instead of Octofitter's default. This sampler is less efficient for unimodal distributions, but is more robust at exploring posteriors with distinct, widely separated peaks.
To install and use Pigeons.jl with Octofitter, type using Pigeons at in the terminal and accept the prompt to install the package. You may have to restart Julia.
octofit_pigeons scales very well across multiple cores. Start julia with julia --threads=auto to make sure you have multiple threads available for sampling.
We now sample from our model using Pigeons:
using Pigeons
chain_pma, pt = octofit_pigeons(model_pma, n_rounds=13, explorer=SliceSampler())
display(chain_pma)[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-52.2943597532602, -44.62558646620697)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.98 4.22 0.331 7.01e+06 -1.53e+04 0 0.736 1 1
4 0 5.52 5.2 0.118 1.6e+05 -431 0 0.654 1 1
8 0 6.11 5.25 0.222 2.49e+05 -1.05e+03 0 0.634 1 1
16 0 6.48 5.46 0.438 4.38e+05 -53.7 7.7e-45 0.615 1 1
32 0 7.97 7.32 0.846 6.37e+05 -27.5 1.93e-11 0.507 1 1
64 1 7.56 5.66 1.87 6.04e+07 -20.7 0.0166 0.574 1 1
128 2 8.43 5.86 3.32 1.09e+08 -22.7 0.0191 0.539 1 1
256 6 8.84 5.93 6.53 2.15e+08 -23 0.134 0.524 1 1
512 22 8.74 5.87 13.2 4.32e+08 -22.3 0.305 0.529 1 1
1.02e+03 43 8.8 6.3 26.4 8.7e+08 -22.6 0.343 0.513 1 1
2.05e+03 88 8.72 6.39 52.3 1.71e+09 -22.8 0.389 0.513 1 1
4.1e+03 177 8.71 6.6 104 3.4e+09 -22.6 0.396 0.506 1 1
8.19e+03 361 8.72 6.8 208 6.76e+09 -22.5 0.406 0.499 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Note that octofit_pigeons took somewhat longer to run than octofit typically does; however, as we will see, it sampled successfully from severally completely disconnected modes in the posterior. That makes it a good fit for sampling from proper motion anomaly and relative astrometry with limited orbital coverage.
Analysis
The first step is to look at the table output above generated by MCMCChains.jl. The rhat column gives a convergence measure. Each parameter should have an rhat very close to 1.000. If not, you may need to run the model for more iterations or tweak the parameterization of the model to improve sampling. The ess column gives an estimate of the effective sample size. The mean and std columns give the mean and standard deviation of each parameter.
The second table summarizes the 2.5, 25, 50, 75, and 97.5 percentiles of each parameter in the model.
Pair Plot
If we wish to examine the covariance between parameters in more detail, we can construct a pair-plot (aka. corner plot).
# Create a corner plot / pair plot.
# We can access any property from the chain specified in Variables
using CairoMakie: Makie
using PairPlots
octocorner(model_pma, chain_pma, small=true)
Notice how there are completely separated peaks? The default Octofitter sample (Hamiltonian Monte Carlo) is capabale of jumping 2-3σ gaps between modes, but such widely separated peaks can cause issues (hence why we used Pigeons in this example).
Posterior Mass vs. Semi-Major Axis
Given that this posterior is quite unconstrained, it is useful to make a simplified plot marginalizing over all orbital parameters besides semi-major axis. We can do this using PairPlots.jl:
using CairoMakie, PairPlots
pairplot(
(; a=chain_pma["b_a"][:], mass=chain_pma["b_mass"][:]) =>
(
PairPlots.Scatter(color=:red),
PairPlots.MarginHist(),
PairPlots.MarginConfidenceLimits()
),
labels=Dict(:mass=>"mass [Mⱼᵤₚ]", :a=>"sma. [au]"),
axis = (;
a = (;
scale=Makie.pseudolog10,
ticks=2 .^ (0:1:6)
)
)
)
Model: PMA & Relative Astrometry
The first orbit fit to only Hipparcos/GAIA data was very unconstrained. We will now add six epochs of relative astrometry (measured from direct images) gathered from the discovery paper.
astrom_like = PlanetRelAstromLikelihood(
(epoch=mjd("2016-12-15"), ra=133., dec=-174., σ_ra=07.0, σ_dec=07., cor=0.2),
(epoch=mjd("2017-03-12"), ra=126., dec=-176., σ_ra=04.0, σ_dec=04., cor=0.3),
(epoch=mjd("2017-03-13"), ra=127., dec=-172., σ_ra=04.0, σ_dec=04., cor=0.1),
(epoch=mjd("2018-02-08"), ra=083., dec=-133., σ_ra=10.0, σ_dec=10., cor=0.4),
(epoch=mjd("2018-11-28"), ra=058., dec=-122., σ_ra=10.0, σ_dec=20., cor=0.3),
(epoch=mjd("2018-12-15"), ra=056., dec=-104., σ_ra=08.0, σ_dec=08., cor=0.2),
)
scatter(astrom_like.table.ra, astrom_like.table.dec)
We use the same model as before, but now condition the planet model B on the astrometry data by adding astrom_like to the end of the @planet defintion.
using OctofitterRadialVelocity
rvlike = PlanetRelativeRVLikelihood(
(epoch=mjd("2008-05-01"), rv=1300, σ_rv=150, inst_idx=1),
(epoch=mjd("2010-02-15"), rv=700, σ_rv=150, inst_idx=1),
(epoch=mjd("2016-03-01"), rv=-2700, σ_rv=150, inst_idx=1),
)
@planet b Visual{KepOrbit} begin
a ~ LogUniform(0.1,400)
e ~ Uniform(0,0.999)
ω ~ UniformCircular()
i ~ Sine()
Ω ~ UniformCircular()
mass = system.M_sec
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,57737.0) # epoch of astrometry
end astrom_like # Note the relative astrometry added here!
@system HD91312_pma_astrom begin
M_pri ~ truncated(Normal(1.61, 0.1), lower=0)
M_sec ~ LogUniform(0.5, 1000) # MJup
M = system.M_pri + system.M_sec*Octofitter.mjup2msol
plx ~ gaia_plx(gaia_id=756291174721509376)
# Priors on the centre of mass proper motion
pmra ~ Normal(-137, 10)
pmdec ~ Normal(2, 10)
end HGCALikelihood(gaia_id=756291174721509376) b
model_pma_astrom = Octofitter.LogDensityModel(HD91312_pma_astrom,autodiff=:ForwardDiff,verbosity=4)
chain_pma_astrom, pt = octofit_pigeons(model_pma_astrom, n_rounds=12, explorer=SliceSampler())[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 14
ℓπcallback(θ): 0.000015 seconds
┌ Info: Tuning autodiff
│ chunk_size = 14
└ t = 1.2433e-5
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 14
∇ℓπcallback(θ): 0.000038 seconds
[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-204.76484306360447, -197.80112130815783)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.39 3.41 0.343 4.03e+06 -2.82e+03 0 0.781 1 1
4 0 5.75 4.22 0.255 1.55e+05 -9.87e+03 0 0.678 1 1
8 0 7.21 6.53 0.5 2.24e+05 -345 5.82e-176 0.557 1 1
16 0 7.21 7.78 0.982 3.88e+05 -151 3.3e-113 0.516 1 1
32 0 9.13 8.62 1.83 5.96e+05 -130 7.73e-23 0.427 1 1
64 0 9.86 7.88 3.64 6.59e+07 -94.3 0.018 0.428 1 1
128 0 9.9 7.39 6.74 1.15e+08 -76.4 0.00186 0.442 1 1
256 0 9.24 7.59 13.2 2.29e+08 -73.8 0.0272 0.457 1 1
512 9 9.6 7.17 25.6 4.54e+08 -60.7 0.0295 0.459 1 1
1.02e+03 18 9.91 7.48 50.8 9.09e+08 -74.1 0.0601 0.439 1 1
2.05e+03 55 10.1 6.66 101 1.82e+09 -74.3 0.145 0.459 1 1
4.1e+03 101 10.1 6.87 202 3.63e+09 -73.7 0.255 0.451 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Model: PMA & Relative Astrometry & RVs
We now add in three additional epochs of stellar RVs.
using OctofitterRadialVelocity
rvlike = StarAbsoluteRVLikelihood(
(epoch=mjd("2008-05-01"), rv=1300, σ_rv=150, inst_idx=1),
(epoch=mjd("2010-02-15"), rv=700, σ_rv=150, inst_idx=1),
(epoch=mjd("2016-03-01"), rv=-2700, σ_rv=150, inst_idx=1),
)
@planet b Visual{KepOrbit} begin
a ~ LogUniform(0.1,400)
e ~ Uniform(0,0.999)
ω ~ UniformCircular()
i ~ Sine()
Ω ~ UniformCircular()
mass = system.M_sec
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,57737.0) # epoch of astrometry
end astrom_like # Note the relative astrometry added here!
@system HD91312_pma_rv_astrom begin
M_pri ~ truncated(Normal(1.61, 0.1), lower=0)
M_sec ~ LogUniform(0.5, 1000) # MJup
M = system.M_pri + system.M_sec*Octofitter.mjup2msol
plx ~ gaia_plx(gaia_id=756291174721509376)
# Priors on the centre of mass proper motion
pmra ~ Normal(-137, 10)
pmdec ~ Normal(2, 10)
jitter ~ truncated(Normal(0,400),lower=0)
rv0 ~ Normal(0,10e3)
end HGCALikelihood(gaia_id=756291174721509376) rvlike b
model_pma_rv_astrom = Octofitter.LogDensityModel(HD91312_pma_rv_astrom,autodiff=:ForwardDiff,verbosity=4)
chain_pma_rv_astrom, pt = octofit_pigeons(model_pma_rv_astrom, n_rounds=12, explorer=SliceSampler())
display(chain_pma_rv_astrom)[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 16
ℓπcallback(θ): 0.000024 seconds (8 allocations: 400 bytes)
┌ Info: Tuning autodiff
│ chunk_size = 16
└ t = 1.7473e-5
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 16
∇ℓπcallback(θ): 0.000051 seconds (8 allocations: 1.203 KiB)
[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-2477.7949595937976, -2465.511272781165)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 7.43 7.17 0.408 2.02e+07 -3.89e+04 0 0.529 1 1
4 0 2.59 3.45 0.356 2.89e+07 -417 0 0.805 1 1
8 0 5.6 5.45 0.642 5.12e+07 -607 0 0.644 1 1
16 0 7.95 7.28 1.33 1.02e+08 -391 2.87e-218 0.509 1 1
32 0 10.4 8.59 2.53 2.02e+08 -165 1.05e-41 0.387 1 1
64 0 10.8 8.04 5.3 4.87e+08 -113 0.0346 0.392 1 1
128 0 10.3 7.76 9.99 9.48e+08 -100 0.0362 0.418 1 1
256 0 10.9 8.29 20.1 1.91e+09 -99.3 0.136 0.382 1 1
512 3 10.6 7.41 40 3.81e+09 -70.3 0.205 0.417 1 1
1.02e+03 8 10.7 8.14 79.6 7.58e+09 -100 0.217 0.394 1 1
2.05e+03 61 10.6 6.33 157 1.5e+10 -102 0.227 0.455 1 1
4.1e+03 140 10.6 6.18 315 3.02e+10 -100 0.261 0.458 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────The mass vs. semi-major axis posterior is now much more constrained:
using CairoMakie, PairPlots
pairplot(
(; a=chain_pma_rv_astrom["b_a"][:], mass=chain_pma_rv_astrom["b_mass"][:]) =>
(
PairPlots.Scatter(color=:red),
PairPlots.MarginHist(),
PairPlots.MarginConfidenceLimits()
),
labels=Dict(:mass=>"mass [Mⱼᵤₚ]", :a=>"sma. [au]"),
)
It is now useful to display the orbits projected onto the plane of the sky using octoplot. This function produces a nine-panel figure showing posterior predictive distributions for velocity (in three dimensions), projected positions vs. time in the plane of the sky, and various other two and three-dimensional views.
octoplot(model_pma_rv_astrom, chain_pma_rv_astrom, show_mass=true)
Model: Relative Astrometry & RVs (no PMA)
There is a final model we should consider: one using the RV and astrometry data, but not the proper motion anomaly:
using OctofitterRadialVelocity
@planet b Visual{KepOrbit} begin
a ~ LogUniform(0.1,400)
e ~ Uniform(0,0.999)
ω ~ UniformCircular()
i ~ Sine()
Ω ~ UniformCircular()
mass = system.M_sec
θ ~ UniformCircular()
tp = θ_at_epoch_to_tperi(system,b,57737.0) # epoch of astrometry
end astrom_like # Note the relative astrometry added here!
@system HD91312_rv_astrom begin
M_pri ~ truncated(Normal(1.61, 0.1), lower=0)
M_sec ~ LogUniform(0.5, 1000) # MJup
M = system.M_pri + system.M_sec*Octofitter.mjup2msol
plx ~ gaia_plx(gaia_id=756291174721509376)
jitter ~ truncated(Normal(0,400),lower=0)
rv0 ~ Normal(0,10e3)
end rvlike b
model_rv_astrom = Octofitter.LogDensityModel(HD91312_rv_astrom,autodiff=:ForwardDiff,verbosity=4)
chain_rv_astrom, pt = octofit_pigeons(model_rv_astrom, n_rounds=12)[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 14
ℓπcallback(θ): 0.000016 seconds (8 allocations: 400 bytes)
┌ Info: Tuning autodiff
│ chunk_size = 14
└ t = 1.1351e-5
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 14
∇ℓπcallback(θ): 0.000027 seconds (8 allocations: 1.109 KiB)
[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-99.42140707996221, -92.19905845133069)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 4.07 4.3 5.91 4.1e+08 -186 0 0.73 0.942 0.967
4 0 4.77 4.6 0.954 6.9e+07 -89.2 8.69e-61 0.698 0.902 0.943
8 0 5.49 6.05 1.32 1.06e+08 -99.7 6.78e-20 0.628 0.914 0.938
16 0 7.56 7.27 2.65 2.09e+08 -85.7 8.82e-05 0.521 0.917 0.932
32 0 8.08 7.56 5.43 4.22e+08 -86.2 0.132 0.496 0.923 0.935
64 0 7.89 4.36 11.9 1.06e+09 -87 0.314 0.605 0.909 0.931
128 4 7.8 4.4 22.2 2.11e+09 -87 0.318 0.606 0.91 0.933
256 11 7.95 4.65 43.9 4.2e+09 -87.2 0.392 0.594 0.915 0.931
512 27 7.89 4.83 88 8.41e+09 -87.1 0.404 0.59 0.913 0.932
1.02e+03 41 7.86 5.4 177 1.69e+10 -87 0.433 0.572 0.917 0.934
2.05e+03 77 8.05 5.65 347 3.32e+10 -87 0.41 0.558 0.916 0.933
4.1e+03 177 7.99 5.55 691 6.62e+10 -87 0.447 0.563 0.916 0.934
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Model Comparison
Let's now display the constraints provided by each data set in a single corner plot
# Create a corner plot / pair plot.
using CairoMakie: Makie
using PairPlots
octocorner(
model_pma,
chain_pma,
chain_pma_astrom,
chain_rv_astrom,
chain_pma_rv_astrom,
small=true, # pass small=false to show all variables
axis=(;
b_a = (;lims=(low=0, high=25))
),
viz=(
PairPlots.MarginDensity(),
PairPlots.Scatter()
)
)
We see that the constraints provided by the PMA, the astrometry, and the radial velocity data all individually overlap, and agree with the joint model constraint. This is means that none of the datasets are in tension with each other, which might suggest an issue with the data or with the modelling assumptions (e.g. single planet).