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.
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 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 = Octofitter.LogDensityModel(HD91312)LogDensityModel for System HD91312 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
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
model = Octofitter.LogDensityModel(HD91312)
chain, pt = octofit_pigeons(model, n_rounds=10)
display(chain)[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-18.285422824717745, -3.6412792008190777)
┌ Warning: Invalid log likelihood encountered. (maxlog=1)
│ θ = (M_pri = 6.046564855054674e11, M_sec = 997.1088383096163, plx = 3.885293999452655e168, pmra = -146.04005537197912, pmdec = 6.228727291579075, M = 6.04656485506419e11, planets = (B = (a = 0.1, e = 0.06987081365056605, ωy = 0.7446947674171549, ωx = -2.059939627306561, i = 1.5901369116596182, Ωy = -58.73457457102934, Ωx = -77.7618354629104, θy = 59.05842566992233, θx = 86.03271073726886, ω = -1.2239020606857682, Ω = -2.2176894970389105, mass = 997.1088383096163, θ = 0.9692122829155981, tp = 57422.99999242064),))
│ llike = -Inf
│ θ_transformed =
│ 14-element Vector{Float64}:
│ 27.127926341169072
│ ⋮
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/logdensitymodel.jl:105
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 3.3 4.69 5.22 3.42e+08 -1.4e+04 0 0.742 0.923 0.963
4 0 3.73 5.75 0.605 8.91e+06 -89.5 0 0.694 0.902 0.937
8 0 5.67 6.94 0.82 6.13e+05 -15.6 1.28e-06 0.593 0.887 0.935
16 0 7.24 7.91 1.69 1.14e+06 -18.9 0.000202 0.511 0.911 0.934
32 0 8.46 8.14 3.76 2.11e+06 -19.3 0.0927 0.464 0.908 0.925
64 0 8.54 5.15 8.35 2.19e+08 -22.2 0.152 0.559 0.912 0.929
128 1 8.51 6 15.2 4.04e+08 -26.9 0.209 0.532 0.906 0.93
256 7 8.73 6.13 30.1 7.47e+08 -26.1 0.264 0.521 0.905 0.925
512 11 8.63 6.69 68.6 1.76e+09 -23.1 0.34 0.506 0.905 0.921
1.02e+03 33 8.76 6.58 138 3.52e+09 -23 0.383 0.505 0.907 0.922
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────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, chain, 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["B_a"][:], mass=chain["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)
)
)
)
Adding 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.
@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 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 = Octofitter.LogDensityModel(HD91312)LogDensityModel for System HD91312 of dimension 14 with fields .ℓπcallback and .∇ℓπcallback
Sample as before:
using Pigeons
model = Octofitter.LogDensityModel(HD91312)
Random.seed!(1)
chain, pt = octofit_pigeons(model, n_rounds=10)
display(chain)[ Info: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└ initial_logpost_range = (-49.32472211299783, -41.63786747733925)
┌ Warning: Invalid log likelihood encountered. (maxlog=1)
│ θ = (M_pri = 2.105840941449194e46, M_sec = 999.9999986279773, plx = 1.2170691597235583e156, pmra = -155.10274442742778, pmdec = 27.014658675895895, M = 2.105840941449194e46, planets = (B = (a = 0.1, e = 0.45823566788527886, ωy = -2.3222539364418036, ωx = 10.008325050738009, i = 2.220446049927349e-16, Ωy = 188.22872389266922, Ωx = 224.5568096427597, θy = 182.82852632140217, θx = 242.3636193547763, ω = 1.7987939781158442, Ω = 0.8731791685230385, mass = 999.9999986279773, θ = 0.9245130025204944, tp = 57737.0),))
│ llike = -Inf
│ θ_transformed =
│ 14-element Vector{Float64}:
│ 106.66362916220822
│ ⋮
└ @ Octofitter ~/work/Octofitter.jl/Octofitter.jl/src/logdensitymodel.jl:105
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 4.3 4.17 5.38 3.33e+08 -9.92e+03 0 0.727 0.947 0.972
4 0 4.71 5.84 1 8.9e+06 -1.08e+03 0 0.66 0.921 0.944
8 0 6.34 6.65 1.42 6.21e+05 -274 6.29e-195 0.581 0.917 0.938
16 0 7.73 6.63 2.97 1.16e+06 -224 1.62e-146 0.537 0.919 0.934
32 0 8.18 8.28 5.85 2.1e+06 -157 7.69e-81 0.469 0.911 0.93
64 1 8.91 4.45 12.3 2.12e+08 18.8 0.00614 0.569 0.918 0.933
128 1 9.25 5.91 24.1 3.92e+08 -75.1 0.00156 0.511 0.919 0.934
256 11 9.6 6.19 45.8 7.48e+08 -74.8 0.00498 0.491 0.913 0.932
512 15 9.81 6.23 88.9 1.44e+09 -74.1 0.0409 0.483 0.914 0.932
1.02e+03 42 9.93 6.13 182 2.94e+09 -74.9 0.0977 0.482 0.913 0.931
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────The mass vs. semi-major axis posterior is now much more constrained:
using CairoMakie, PairPlots
pairplot(
(; a=chain["B_a"][:], mass=chain["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, chain)
Finally we display the posterior as a 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, chain, small=true)