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.
Model: PMA 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.
Retrieving the HGCA
To start, we retrieve the HGCA data for this object.
hgca_like = HGCALikelihood(gaia_id=756291174721509376)For this example, we will speed things up by treating the HGCA measurements as instantaenous. The default is to average measurements 25 times over the duration of each mission. For real inference, you'll want to increase this or leave it at the default value, unless you have reason to believe the companion has an orbital period much greater than 20 years or so.
hgca_like = HGCALikelihood(gaia_id=756291174721509376, N_ave=1)HGCALikelihood Table with 3 columns and 4 rows:
epoch meas inst
┌────────────────────
1 │ 48345.7 ra hip
2 │ 48500.8 dec hip
3 │ 57423.6 ra gaia
4 │ 57503.4 dec gaiaPlanet Model
@planet b Visual{KepOrbit} begin
a ~ LogUniform(0.1,20)
e ~ Uniform(0,0.999)
ω ~ Uniform(0, 2pi)
i ~ Sine() # The Sine() distribution is defined by Octofitter
Ω ~ Uniform(0, 2pi)
mass = system.M_sec
θ ~ Uniform(0, 2pi)
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)
ω Distributions.Uniform{Float64}(a=0.0, b=6.283185307179586)
i Sine()
Ω Distributions.Uniform{Float64}(a=0.0, b=6.283185307179586)
θ Distributions.Uniform{Float64}(a=0.0, b=6.283185307179586)
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.1) # 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 hgca_like b
model_pma = Octofitter.LogDensityModel(HD91312_pma)LogDensityModel for System HD91312_pma of dimension 11 and 4 epochs 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)┌ Warning: Module Octofitter with build ID fafbfcfd-a27b-b7ca-0000-006fbc8fe2e9 is missing from the cache.
│ This may mean Octofitter [daf3887e-d01a-44a1-9d7e-98f15c5d69c9] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2011
┌ Warning: Module ForwardDiff with build ID fafbfcfd-dc20-ae10-0000-0075d95f115c is missing from the cache.
│ This may mean ForwardDiff [f6369f11-7733-5829-9624-2563aa707210] does not support precompilation but is imported by a module that does.
└ @ Base loading.jl:2011
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: Determining initial positions using pathfinder, around that location.
┌ Warning: ErrorException("Could not find starting point within 0.001 - 0.999 quantiles of priors.")
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/initialization.jl:302
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-14.466201297627096, -14.466201297627096)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 4.02 1.29 0.892 2.92e+07 -2.76e+04 0 0.829 1 1
4 0 4.95 7.61 0.0908 1.3e+05 -25.5 0 0.595 1 1
8 0 6.8 8.51 0.17 2.16e+05 -20.4 0.0662 0.506 0.989 0.999
16 0 7.06 7.67 0.335 3.71e+05 -24.8 6.01e-08 0.525 0.998 1
32 0 7.65 8.39 0.68 6.22e+05 -21.2 0.037 0.483 1 1
64 2 8.18 5.29 1.59 5.38e+07 -23 0.0545 0.566 0.999 1
128 5 7.95 5.41 2.7 9.07e+07 -20.3 0.326 0.569 1 1
256 10 8.62 5.41 5.41 1.79e+08 -21.2 0.359 0.547 1 1
512 21 8.24 5.62 10.6 3.5e+08 -21 0.395 0.553 1 1
1.02e+03 46 8.33 5.87 21.1 7.02e+08 -21.1 0.371 0.542 1 1
2.05e+03 106 8.32 5.7 41.7 1.38e+09 -21.2 0.393 0.548 1 1
4.1e+03 223 8.33 5.55 84 2.78e+09 -21.2 0.417 0.552 1 1
8.19e+03 462 8.29 5.75 168 5.54e+09 -21.1 0.427 0.547 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.
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
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).
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),
jitter=:jitter_1,
instrument_name=""
)
@planet b Visual{KepOrbit} begin
a ~ LogUniform(0.1,400)
e ~ Uniform(0,0.999)
ω ~ Uniform(0, 2pi)
i ~ Sine()
Ω ~ Uniform(0, 2pi)
mass = system.M_sec
θ ~ Uniform(0, 2pi)
tp = θ_at_epoch_to_tperi(system,b,57737.0) # epoch of astrometry
jitter_1 = 0.0
end astrom_like # Note the relative astrometry added here!
@system HD91312_pma_astrom begin
M_pri ~ truncated(Normal(1.61, 0.1), lower=0.1)
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 hgca_like b
model_pma_astrom = Octofitter.LogDensityModel(HD91312_pma_astrom,autodiff=:ForwardDiff,verbosity=4)
using Pigeons
chain_pma_astrom, pt = octofit_pigeons(model_pma_astrom, n_rounds=12, explorer=SliceSampler())[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 11
ℓπcallback(θ, args...): 0.000025 seconds (9 allocations: 1.688 KiB)
┌ Info: Tuning autodiff
│ chunk_size = 11
└ t = 1.035e-5
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 11
∇ℓπcallback(θ): 0.000055 seconds (1 allocation: 144 bytes)
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: Determining initial positions using pathfinder, around that location.
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-65.47425181523383, -55.24407487145618)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 5.07 2.53 0.322 6.35e+06 -1.89e+04 0 0.755 1 1
4 0 5.02 4.23 0.186 1.45e+05 -256 0 0.702 0.992 1
8 0 5.73 6.47 0.342 2.34e+05 -131 0 0.606 0.996 1
16 0 7.61 7.42 0.66 3.57e+05 -181 2.95e-98 0.515 0.998 1
32 0 8.56 8.73 1.32 5.47e+05 -76.6 0.000521 0.442 0.999 1
64 1 8.7 5.05 2.77 5.32e+07 -64.4 2.25e-10 0.557 0.997 0.999
128 1 8.9 6.5 5.17 9.3e+07 -77.8 0.00167 0.503 0.999 1
256 6 9.58 6.72 10.2 1.85e+08 -73.7 0.0158 0.474 1 1
512 16 9.52 6.71 20.4 3.68e+08 -71.9 0.0848 0.476 1 1
1.02e+03 24 9.8 6.93 39.9 7.2e+08 -71.9 0.122 0.46 1 1
2.05e+03 71 9.78 6.65 79.3 1.44e+09 -72.2 0.21 0.47 1 1
4.1e+03 167 9.84 6.51 160 2.9e+09 -72.1 0.309 0.472 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────octoplot(model_pma_astrom, chain_pma_astrom)
Model: PMA & Relative Astrometry & RVs
We now add in three additional epochs of stellar RVs.
using OctofitterRadialVelocity
rvlike = MarginalizedStarAbsoluteRVLikelihood(
(epoch=mjd("2008-05-01"), rv=1300, σ_rv=150),
(epoch=mjd("2010-02-15"), rv=700, σ_rv=150),
(epoch=mjd("2016-03-01"), rv=-2700, σ_rv=150),
instrument_name="SOPHIE",
jitter=:jitter_1
)
@planet b Visual{KepOrbit} begin
a ~ LogUniform(0.1,400)
e ~ Uniform(0,0.999)
ω ~ Uniform(0, 2pi)
i ~ Sine()
Ω ~ Uniform(0, 2pi)
mass = system.M_sec
θ ~ Uniform(0, 2pi)
tp = θ_at_epoch_to_tperi(system,b,57737.0) # epoch of astrometry
end astrom_like ObsPriorAstromONeil2019(astrom_like) # Note the relative astrometry added here!
@system HD91312_pma_rv_astrom begin
M_pri ~ truncated(Normal(1.61, 0.1), lower=0.1)
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_1 = 0.0
end hgca_like 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 = 11
ℓπcallback(θ, args...): 0.000036 seconds (9 allocations: 1.969 KiB)
┌ Info: Tuning autodiff
│ chunk_size = 11
└ t = 1.4216e-5
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 11
∇ℓπcallback(θ): 0.000060 seconds (1 allocation: 144 bytes)
┌ Warning: This model has priors that cannot be sampled IID.
└ @ OctofitterPigeonsExt ~/work/Octofitter.jl/Octofitter.jl/ext/OctofitterPigeonsExt/OctofitterPigeonsExt.jl:104
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: Determining initial positions using pathfinder, around that location.
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-68.77372786169532, -63.30752129965552)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.49 2.11 0.346 3.83e+06 -2.48e+04 0 0.819 1 1
4 0 5.48 6.61 0.306 1.43e+05 -279 0 0.61 1 1
8 0 6.2 9.19 0.556 2.1e+05 -99.3 0 0.504 1 1
16 0 7.33 9.83 1.09 3.38e+05 -102 2.22e-118 0.446 1 1
32 0 8.02 10.3 2.12 5.18e+05 -104 4.39e-102 0.409 1 1
64 1 9.05 3.54 3.99 5.42e+07 -126 2.98e-13 0.594 1 1
128 8 9.4 3.81 7.92 9.91e+07 -98.4 3.1e-13 0.574 1 1
256 16 9.94 3.89 15.8 1.99e+08 -98.5 0.000982 0.554 0.998 1
512 35 10.1 3.75 31.3 3.94e+08 -98.6 0.0104 0.555 1 1
1.02e+03 52 10.3 5.24 61.7 7.76e+08 -98.2 0.0227 0.5 1 1
2.05e+03 98 10.4 5.39 122 1.53e+09 -98.5 0.0774 0.491 1 1
4.1e+03 197 10.7 5.33 241 3.04e+09 -98.5 0.163 0.484 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)
ω ~ Uniform(0, 2pi)
i ~ Sine()
Ω ~ Uniform(0, 2pi)
mass = system.M_sec
θ ~ Uniform(0, 2pi)
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.1)
M_sec ~ LogUniform(0.5, 1000) # MJup
M = system.M_pri + system.M_sec*Octofitter.mjup2msol
plx ~ gaia_plx(gaia_id=756291174721509376)
jitter_1 = 0.0
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 = 9
ℓπcallback(θ, args...): 0.000029 seconds (9 allocations: 608 bytes)
┌ Info: Tuning autodiff
│ chunk_size = 9
└ t = 8.235e-6
┌ Info: Selected auto-diff chunk size
└ ideal_chunk_size = 9
∇ℓπcallback(θ): 0.000043 seconds (1 allocation: 128 bytes)
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: Determining initial positions using pathfinder, around that location.
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-75.95577756347852, -68.34018577641103)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 2.61 2.2 0.354 8.64e+06 -1.36e+04 0 0.845 1 1
4 0 5.36 4.64 0.143 1.43e+05 -103 7.04e-41 0.677 1 1
8 0 4.9 5.53 0.263 2.45e+05 -85 4.24e-10 0.664 1 1
16 0 6.17 7.08 0.495 3.79e+05 -87.8 0.0354 0.573 1 1
32 0 7.24 6.84 0.954 6.42e+05 -88 0.119 0.546 1 1
64 2 7.36 4.16 2.07 4.35e+07 -89.4 0.242 0.628 0.999 1
128 4 7.38 4.8 3.76 7.4e+07 -90.2 0.378 0.607 0.999 1
256 14 7.3 4.61 7.51 1.48e+08 -89.9 0.407 0.616 0.999 1
512 40 7.25 4.4 14.8 2.93e+08 -90.3 0.474 0.624 0.999 1
1.02e+03 80 7.11 4.12 29.7 5.86e+08 -90.2 0.493 0.638 1 1
2.05e+03 160 7.22 4.48 59.5 1.17e+09 -90 0.482 0.623 1 1
4.1e+03 329 7.22 4.5 119 2.34e+09 -90 0.503 0.622 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────octoplot(model_rv_astrom, chain_rv_astrom)
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=false,
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).