Detection Limits
This guide shows how to calculate detection limits, in mass, or in photometry, as a function of orbital parameters for different combinations of data.
There are a few use cases for this:
- Mass limit vs semi-major axis given one or more images and/or contrast curves
- Mass limit vs semi-major axis given an RV non-detection
- Mass limit vs semi-major axis given proper motion anomaly from the GAIA-Hipparcos Catalog of Accelerations
- Any combination of the above
We will once more use some sample data from the system HD 91312 A & B discovered by SCExAO.
using Octofitter
using OctofitterImages
using OctofitterRadialVelocity
using Distributions
using Pigeons
using CairoMakie
using PairPlotsPhotometry Model
We will need to decide on an atmosphere model to map image intensities into mass. Here we use the Sonora Bobcat cooling and atmosphere models which will be auto-downloaded by Octofitter:
const cooling_tracks = Octofitter.sonora_cooling_interpolator()
const sonora_temp_mass_L = Octofitter.sonora_photometry_interpolator(:Keck_L′)(::Octofitter.var"#model_interpolator#396"{Interpolations.FilledExtrapolation{Float64, 2, Interpolations.ScaledInterpolation{Float64, 2, Interpolations.BSplineInterpolation{Float64, 2, Matrix{Float64}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Tuple{StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}}, Interpolations.BSpline{Interpolations.Linear{Interpolations.Throw{Interpolations.OnGrid}}}, Float64}, Float64, Float64, Float64, Float64}) (generic function with 1 method)Proper Motion Anomaly Data
We start by defining and sampling from a model that only includes proper motion anomaly data from the HGCA:
B = Planet(
name="B",
basis=Visual{KepOrbit},
likelihoods=[],
variables=@variables begin
a ~ LogUniform(1, 65)
e ~ Uniform(0,0.9)
ω ~ Uniform(0,2pi)
i ~ Sine() # The Sine() distribution is defined by Octofitter
Ω ~ Uniform(0,pi)
mass = system.M_sec
θ ~ Uniform(0,2pi)
tp = θ_at_epoch_to_tperi(θ, 57423.0; M=system.M, e, a, i, ω, Ω) # epoch of GAIA measurement
end
)
HD91312_pma = System(
name="HD91312_pma",
companions=[B],
likelihoods=[HGCAInstantaneousLikelihood(gaia_id=6166183842771027328)],
variables=@variables begin
M_pri ~ truncated(Normal(0.95, 0.05), lower=0.1) # Msol
M_sec ~ LogUniform(0.2, 65) # MJup
M = M_pri + M_sec*Octofitter.mjup2msol # Msol
plx ~ gaia_plx(gaia_id=6166183842771027328)
# Priors on the center of mass proper motion
pmra ~ Normal(0, 1000)
pmdec ~ Normal(0, 1000)
end
)
model_pma = Octofitter.LogDensityModel(HD91312_pma)LogDensityModel for System HD91312_pma of dimension 11 and 4 epochs with fields .ℓπcallback and .∇ℓπcallback
Sample:
using Pigeons
chain_pma, pt = octofit_pigeons(model_pma, n_chains=16, n_chains_variational=16, n_rounds=12);[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: Performing global optimization with 11 parameters (0 initial parameter value provided).
┌ Info: Found sample of initial positions
│ logpost_range = (-24.340362583978493, -15.749320753160012)
└ mean_logpost = -18.223385123478558
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 2.89 2.67 0.235 4.15e+06 -1.85e+08 0 0.821 1 1
4 0 4.22 5.26 0.0705 1.43e+05 -5.36e+04 0 0.694 1 1
8 0 4.74 7.83 0.136 2.35e+05 -704 0 0.595 1 1
16 0 4.85 6.93 0.27 4.1e+05 -1.16e+04 0 0.62 1 1
32 0 5.52 7.4 0.54 7.06e+05 -123 0 0.583 1 1
64 3 5.93 2.1 1.23 4.41e+07 -17.2 0 0.741 1 1
128 6 6.89 3.03 2.1 7.41e+07 -17.3 0 0.68 1 1
256 15 7.39 3.83 4.21 1.47e+08 -17.5 5e-290 0.638 1 1
512 50 7.83 3.06 8.31 2.9e+08 -17.3 5.74e-25 0.649 1 1
1.02e+03 93 8.05 3.09 16.8 5.81e+08 -17.2 0.00196 0.641 1 1
2.05e+03 199 8.54 2.98 33.7 1.16e+09 -17.1 0.00209 0.629 1 1
4.1e+03 367 8.88 3.38 67.8 2.33e+09 -17.3 0.00422 0.605 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Plot the marginal mass vs. semi-major axis posterior with contours using PairPlots.jl:
pairplot(
PairPlots.Series(
(;
sma=log.(chain_pma[:B_a][:],),
mass=log.(chain_pma[:B_mass][:]),
),
label="PMA",
color=Makie.wong_colors()[1],
)=>(
PairPlots.Scatter(markersize=3,alpha=0.35),
PairPlots.Contour(sigmas=[1,3]),
PairPlots.MarginStepHist(),
),
labels=Dict(
:sma=>"log Semi-major axis [au]",
:mass=>"log Mass [Mⱼᵤₚ]"
)
)
Image Data
using AstroImages
# download(
# "https://github.com/sefffal/Octofitter.jl/raw/main/docs/image-examples-1.fits",
# "image-examples-1.fits"
# )
# Or multi-extension FITS (this example)
image = AstroImages.load("image-examples-1.fits").*2e-7 # units of contrast
img_dat_table = Table([
(image=AstroImages.recenter(image), platescale=4.0, epoch=57423.6),
])
image_data = ImageLikelihood(
img_dat_table,
name="imgdat-sim",
variables=@variables begin
# Planet flux in image units -- could be contrast, mags, Jy, or arb. as long as it's consistent with the units of the data you provide
flux = planet.L
# The following are optional parameters for marginalizing over instrument systematics:
# Platescale uncertainty multiplier [could use: platescale ~ truncated(Normal(1, 0.01), lower=0)]
platescale = 1.0
# North angle offset in radians [could use: northangle ~ Normal(0, deg2rad(1))]
northangle = 0.0
end
)OctofitterImages.ImageLikelihood Table with 5 columns and 1 row:
image platescale epoch contrast ⋯
┌───────────────────────────────────────────────────────────────────
1 │ [-6.80461e-8 -6.247… 4.0 57423.6 65-element extrapol… ⋯B = Planet(
name="B",
basis=Visual{KepOrbit},
likelihoods=[image_data],
variables=@variables begin
a ~ LogUniform(1, 65)
e ~ Uniform(0,0.9)
ω ~ Uniform(0,2pi)
i ~ Sine() # The Sine() distribution is defined by Octofitter
Ω ~ Uniform(0,pi)
mass = system.M_sec
# Calculate planet temperature from cooling track and planet mass variable
tempK = $cooling_tracks(system.age, mass)
# Calculate absolute magnitude
abs_mag_L = $sonora_temp_mass_L(tempK, mass)
# Deal with out-of-grid values by clamping to grid max and min
abs_mal_L′ = if isfinite(abs_mag_L)
abs_mag_L
elseif mass > 10
8.2 # jump to absurdly bright
else
16.7 # jump to absurdly dim
end
# Calculate relative magnitude
rel_mag_L = abs_mal_L′ - system.rel_mag + 5log10(1000/system.plx)
# Convert to contrast (same units as image)
L = 10.0^(rel_mag_L/-2.5)
θ ~ Uniform(0,2pi)
tp = θ_at_epoch_to_tperi(θ, 57423.6; M=system.M, e, a, i, ω, Ω)
end
)
HD91312_img = System(
name="HD91312_img",
companions=[B],
likelihoods=[],
variables=@variables begin
# age ~ truncated(Normal(40, 15),lower=0, upper=200)
age = 10
M_pri ~ truncated(Normal(0.95, 0.05), lower=0.1) # Msol
# Mass of secondary
# Make sure to pick only a mass range that is covered by your models
M_sec ~ LogUniform(0.55, 65) # MJup
M = M_pri + M_sec*Octofitter.mjup2msol # Msol
plx ~ gaia_plx(gaia_id=6166183842771027328)
# Priors on the center of mass proper motion
# pmra ~ Normal(0, 1000)
# pmdec ~ Normal(0, 1000)
rel_mag = 5.65
end
)
model_img = Octofitter.LogDensityModel(HD91312_img)LogDensityModel for System HD91312_img of dimension 9 and 1 epochs with fields .ℓπcallback and .∇ℓπcallback
using Pigeons
chain_img, pt = octofit_pigeons(model_img, n_chains=5, n_chains_variational=5, n_rounds=7)(chain = MCMC chain (128×26×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))Plot mass vs. semi-major axis posterior:
vis_layers = (
PairPlots.Contour(sigmas=[1,3]),
PairPlots.MarginStepHist(),
)
pairplot(
PairPlots.Series(
(;
sma=log.(chain_pma[:B_a][:],),
mass=log.(chain_pma[:B_mass][:]),
),
label="PMA",
color=Makie.wong_colors()[1],
)=>vis_layers,
PairPlots.Series(
(;
sma=log.(chain_img[:B_a][:],),
mass=log.(chain_img[:B_mass][:]),
),
label="IMG",
color=Makie.wong_colors()[2],
)=>vis_layers,
labels=Dict(
:sma=>"log Semi-major axis [au]",
:mass=>"log Mass [Mⱼᵤₚ]"
)
)
Image and PMA data
B = Planet(
name="B",
basis=Visual{KepOrbit},
likelihoods=[image_data],
variables=@variables begin
a ~ LogUniform(1, 65)
e ~ Uniform(0,0.9)
ω ~ Uniform(0,2pi)
i ~ Sine() # The Sine() distribution is defined by Octofitter
Ω ~ Uniform(0,pi)
mass = system.M_sec
# Calculate planet temperature from cooling track and planet mass variable
tempK = $cooling_tracks(system.age, mass)
# Calculate absolute magnitude
abs_mag_L = $sonora_temp_mass_L(tempK, mass)
# Deal with out-of-grid values by clamping to grid max and min
abs_mal_L′ = if isfinite(abs_mag_L)
abs_mag_L
elseif mass > 10
8.2 # jump to absurdly bright
else
16.7 # jump to absurdly dim
end
# Calculate relative magnitude
rel_mag_L = abs_mal_L′ - system.rel_mag + 5log10(1000/system.plx)
# Convert to contrast (same units as image)
L = 10.0^(rel_mag_L/-2.5)
# L ~ Uniform(0,1)
θ ~ Uniform(0,2pi)
tp = θ_at_epoch_to_tperi(θ, 57423.6; M=system.M, e, a, i, ω, Ω)
end
)
HD91312_both = System(
name="HD91312_both",
companions=[B],
likelihoods=[HGCAInstantaneousLikelihood(gaia_id=6166183842771027328)],
variables=@variables begin
# age ~ truncated(Normal(40, 15),lower=0, upper=200)
age = 10
M_pri ~ truncated(Normal(0.95, 0.05), lower=0.1) # Msol
# Mass of secondary
# Make sure to pick only a mass range that is covered by your models
M_sec ~ LogUniform(0.55, 65) # MJup
M = M_pri + M_sec*Octofitter.mjup2msol # Msol
plx ~ gaia_plx(gaia_id=6166183842771027328)
# Priors on the center of mass proper motion
pmra ~ Normal(0, 1000)
pmdec ~ Normal(0, 1000)
rel_mag = 5.65
end
)
model_both = Octofitter.LogDensityModel(HD91312_both)
using Pigeons
chain_both, pt = octofit_pigeons(model_both,n_chains=5,n_chains_variational=5,n_rounds=10)(chain = MCMC chain (1024×28×1 Array{Float64, 3}), pt = PT(checkpoint = false, ...))Compare all three posteriors to see limits:
vis_layers = (
PairPlots.Contour(sigmas=[1,3]),
PairPlots.MarginStepHist(),
)
pairplot(
PairPlots.Series(
(;
sma=log.(chain_pma[:B_a][:],),
mass=log.(chain_pma[:B_mass][:]),
),
label="PMA",
color=Makie.wong_colors()[1],
)=>vis_layers,
PairPlots.Series(
(;
sma=log.(chain_img[:B_a][:],),
mass=log.(chain_img[:B_mass][:]),
),
label="IMG",
color=Makie.wong_colors()[2],
)=>vis_layers,
PairPlots.Series(
(;
sma=log.(chain_both[:B_a][:],),
mass=log.(chain_both[:B_mass][:]),
),
label="IMG + PMA",
color=Makie.wong_colors()[3],
)=>vis_layers,
labels=Dict(
:sma=>"log₂ Semi-major axis [au]",
:mass=>"log₂ Mass [Mⱼᵤₚ]"
)
)