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 Hipparcos-Gaia astrometry
- Any combination of the above
The recipe is the same in every case: build a model in which the companion's brightness is derived from its mass through an evolutionary model, fit the data, and read off which masses survive.
The star is Gaia DR3 6166183842771027328 (HIP 65808), a ~32 pc solar-type dwarf. It is chosen because its proper motion is a non-detection: the HGCA χ² is 17.8 and every Hipparcos/Hipparcos–Gaia/DR3 channel sits within ~3σ of the long-term trend, so the PMA posterior below is a genuine upper limit — the broad envelope the other data sets are compared against. Do not retarget this page onto a star with a real astrometric companion (HD 91312, used by Fit Proper Motion Anomaly, has a 16σ anomaly): the PMA posterior then pins itself to the top of the mass prior and the page stops demonstrating anything.
So that the page builds offline it reads the Hipparcos–Gaia catalog subset and the cached scan-law forecast that ship with Octofitter's tests, together with the example L′ contrast image from Fitting Images — which is a stand-in, not a real observation of this star. For a real target you simply drop the catalog= and forecast_table= keywords and Octofitter fetches both itself.
using Octofitter
using OctofitterImages
using Distributions
using AstroImages
using CairoMakie
using PairPlots
using Statistics
using PigeonsPhotometry 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#628"{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, Float64}) (generic function with 1 method)There is a single mass unit throughout, M⊙, and these interpolators take M⊙ by default. mjup is a plain multiplicative constant, so a Jupiter-mass threshold is written 10mjup, not 10 — a bare 10 means ten solar masses, which lands off the end of the grid and quietly returns NaN. Pass mass_unit=:Mjup to the constructors if you would rather keep a script in Jupiter masses.
Proper Motion Anomaly Data
We start by defining and sampling from a model that only includes the Hipparcos-Gaia proper motion anomaly.
HGCAObs is a helper over the joint Gaia-Hipparcos likelihood G23HObs, restricted to the six HGCA channels. It models the actual Gaia and Hipparcos scan epochs; freeze_epochs=true below fixes the epoch selection, which is much faster and is what makes a detection-limit sweep practical.
A = Body(
name="A",
variables=@variables begin
mass = system.M_pri
flux_G = 1.0 # the host defines the contrast scale in each band
flux_Hp = 1.0
end
)
B = Body(
name="B",
about=A,
variables=@variables begin
mass = system.M_sec # M⊙
a ~ LogUniform(1, 65)
e ~ Uniform(0,0.9)
ω ~ Uniform(0,2pi)
i ~ Sine() # The Sine() distribution is defined by Octofitter
Ω ~ Uniform(0,pi)
θ ~ Uniform(0,2pi)
epoch = 57423.0 # epoch of the Gaia measurement
flux_G = 0.0 # dark to Gaia...
flux_Hp = 0.0 # ...and to Hipparcos
end
)
gaia_id = 6166183842771027328 # HIP 65808 -- see the note above
pma = HGCAObs(;
gaia_id = gaia_id,
target = A,
blends = (B,),
ref = Barycentre,
freeze_epochs = true, # faster but approximate; set false for real use
)┌ Info: Removed forecast transits in data gaps of every applicable release.
└ n_removed = 1
┌ Info: Count of missed or rejected transits:
└ dr3 = 3
┌ Info: DR2/DR3 epoch selection
│ n2_win = 13
│ n_tail = 14
│ n_dr2_total = 13
└ n_dr2_distinct_range = (13, 13)The system block supplies the frame. G23HObs (and therefore HGCAObs) needs the reference point's proper motion, which means a full absolute frame: plx, ra, dec, pmra, pmdec, rv, ref_epoch. Declaring only some of them is a build-time error: Octofitter will not guess at half a frame.
cat = pma.catalog
HD_pma = System(
name="limits_pma",
bodies=[A, B],
observations=[pma],
variables=@variables begin
M_pri ~ truncated(Normal(0.95, 0.05), lower=0.1) # M⊙
M_sec ~ LogUniform(0.2mjup, 65mjup) # M⊙ -- note the units!
plx ~ truncated(Normal(cat.parallax, cat.parallax_error), lower=0.1)
ra = $(cat.ra)
dec = $(cat.dec)
# Wide priors on the centre-of-mass proper motion: the data constrain these.
pmra ~ Uniform(cat.pmra_dr3 - 100, cat.pmra_dr3 + 100)
pmdec ~ Uniform(cat.pmdec_dr3 - 100, cat.pmdec_dr3 + 100)
rv = $(isnan(cat.radial_velocity) ? 0.0 : cat.radial_velocity * 1e3) # m/s
ref_epoch = $(Octofitter.meta_gaia_DR3.ref_epoch_mjd)
end
)
model_pma = Octofitter.LogDensityModel(HD_pma)LogDensityModel for System limits_pma of dimension 11 and 6 epochs with fields .ℓπcallback and .∇ℓπcallback
Sample:
init_pma = initialize!(model_pma)
chain_pma, pt = octofit_pigeons(model_pma, n_chains=16, n_chains_variational=16, n_rounds=12)┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Info: Found sample of initial positions
│ logpost_range = (-7.799110919466683, 1.965947425909647)
└ mean_logpost = -0.5553473122308529
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [limits_pma] observing_geometry = true (user)
[ Info: [limits_pma] barycentric_lighttime = true (user)
[ Info: [limits_pma] observing_geometry = true (user)
[ Info: [limits_pma] barycentric_lighttime = true (user)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 2.7 2.66 1.57 1.55e+09 -3.38e+05 0 0.827 1 1
4 0 3.15 4.81 0.784 2.95e+09 -1.97e+04 0 0.743 1 1
8 0 4.47 5.67 1.55 5.73e+09 -1.74e+03 0 0.673 1 1
16 0 5.95 6.09 3.04 1.14e+10 -151 0 0.612 1 1
32 0 6.01 6.66 6.73 2.24e+10 -28.8 5.87e-21 0.591 1 1
64 3 6.21 2.5 13.5 4.62e+10 -12 1.35e-14 0.719 0.999 1
128 16 6.94 2.65 27.6 9.27e+10 -12.4 0.00357 0.691 1 1
256 25 7.28 2.47 54.1 1.84e+11 -12.2 0.0152 0.685 1 1
512 71 8.01 2.45 109 3.66e+11 -12 0.0834 0.662 1 1
1.02e+03 160 8.28 2.32 205 7.3e+11 -12.2 0.179 0.658 1 1
2.05e+03 301 8.25 2.41 397 1.46e+12 -12 0.324 0.656 1 1
4.1e+03 630 8.18 2.4 758 2.91e+12 -12.1 0.434 0.659 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────Plot the marginal mass vs. semi-major axis posterior with contours using PairPlots.jl. Note that B_mass is in solar masses, so we convert for the plot:
pairplot(
PairPlots.Series(
(;
sma=log.(chain_pma[:B_a][:],),
mass=log.(chain_pma[:B_mass][:] ./ mjup),
),
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
Now the same star, with an L′-band contrast image instead. This is where the flux/band unification pays off: the companion's brightness is a body variable, flux_L, derived from its mass — so the image constrains the mass directly.
download(
"https://github.com/sefffal/Octofitter.jl/raw/main/docs/image-examples-1.fits",
"image-examples-1.fits"
)
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),
])B_img = Body(
name="B",
about=A,
variables=@variables begin
mass = system.M_sec
# Calculate companion temperature from the cooling track and its mass
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.
# NB: the threshold is `10mjup`, not `10` -- masses are in M⊙.
abs_mag_L′ = if isfinite(abs_mag_L)
abs_mag_L
elseif mass > 10mjup
8.2 # jump to absurdly bright
else
16.7 # jump to absurdly dim
end
# Calculate relative magnitude
rel_mag_L = abs_mag_L′ - system.rel_mag + 5log10(1000/system.plx)
# Convert to contrast -- the same units as the image. This
# 10^(-0.4 Δmag) step is what makes `flux_L` a *linear* flux.
flux_L = 10.0^(rel_mag_L/-2.5)
a ~ LogUniform(1, 65)
e ~ Uniform(0,0.9)
ω ~ Uniform(0,2pi)
i ~ Sine()
Ω ~ Uniform(0,pi)
θ ~ Uniform(0,2pi)
epoch = 57423.6
end
)
image_data = ImageObs(
img_dat_table,
targets = (B_img,),
ref = A,
band = :L,
name = "imgdat-sim",
variables=@variables begin
# 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
)[ Info: Measuring contrast from imageA_img = Body(
name="A",
variables=@variables begin
mass = system.M_pri
end
)
HD_img = System(
name="limits_img",
bodies=[A_img, B_img],
observations=[image_data],
variables=@variables begin
# age ~ truncated(Normal(40, 15),lower=0, upper=200)
age = 10 # Myr
M_pri ~ truncated(Normal(0.95, 0.05), lower=0.1) # M⊙
# Mass of the secondary.
# Make sure to pick only a mass range that is covered by your models.
M_sec ~ LogUniform(0.55mjup, 65mjup) # M⊙
plx ~ truncated(Normal(cat.parallax, cat.parallax_error), lower=0.1)
rel_mag = 5.65
end
)
model_img = Octofitter.LogDensityModel(HD_img)LogDensityModel for System limits_img of dimension 9 and 1 epochs with fields .ℓπcallback and .∇ℓπcallback
init_img = initialize!(model_img)
chain_img, pt = octofit_pigeons(model_img, n_chains=5, n_chains_variational=5, n_rounds=7)┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 9
└ num_fixed = 0
┌ Info: Found sample of initial positions
│ logpost_range = (-10.424966939147886, 18.619440213736848)
└ mean_logpost = -5.2424150442189115
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [limits_img] observing_geometry = true (user)
[ Info: [limits_img] barycentric_lighttime = true (user)
[ Info: [limits_img] observing_geometry = true (user)
[ Info: [limits_img] barycentric_lighttime = true (user)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 1.09 1.44 0.874 1.02e+07 -282 0 0.719 1 1
4 0 1.53 2 0.022 9.87e+04 -394 0 0.607 1 1
8 0 2.13 2.27 0.0392 1.72e+05 10.9 0.029 0.511 1 1
16 0 1.93 2.09 0.0838 3.12e+05 5.29 0.000227 0.553 1 1
32 0 2.17 2.26 0.175 5.69e+05 7.81 0.000296 0.507 1 1
64 2 2.04 1.74 0.791 2.01e+07 93 0.0337 0.58 1 1
128 3 2.2 2.36 0.539 1.81e+07 8.49 0.0159 0.494 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────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][:] ./ mjup),
),
label="PMA",
color=Makie.wong_colors()[1],
)=>vis_layers,
PairPlots.Series(
(;
sma=log.(chain_img[:B_a][:],),
mass=log.(chain_img[:B_mass][:] ./ mjup),
),
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
Combining the two is just a matter of listing both observations on the system. Observations are a flat list and each names its own references, so there is nothing else to wire up.
pma_joint = HGCAObs(;
gaia_id = gaia_id,
target = A,
blends = (B_img,),
ref = Barycentre,
freeze_epochs = true, # faster but approximate; set false for real use
)
HD_both = System(
name="limits_both",
bodies=[A, B_img],
observations=[image_data, pma_joint],
variables=@variables begin
age = 10
M_pri ~ truncated(Normal(0.95, 0.05), lower=0.1)
M_sec ~ LogUniform(0.55mjup, 65mjup)
rel_mag = 5.65
plx ~ truncated(Normal(cat.parallax, cat.parallax_error), lower=0.1)
ra = $(cat.ra)
dec = $(cat.dec)
pmra ~ Uniform(cat.pmra_dr3 - 100, cat.pmra_dr3 + 100)
pmdec ~ Uniform(cat.pmdec_dr3 - 100, cat.pmdec_dr3 + 100)
rv = $(isnan(cat.radial_velocity) ? 0.0 : cat.radial_velocity * 1e3)
ref_epoch = $(Octofitter.meta_gaia_DR3.ref_epoch_mjd)
end
)
model_both = Octofitter.LogDensityModel(HD_both)
init_both = initialize!(model_both)
chain_both, pt = octofit_pigeons(model_both, n_chains=5, n_chains_variational=5, n_rounds=10)┌ Info: Removed forecast transits in data gaps of every applicable release.
└ n_removed = 1
┌ Info: Count of missed or rejected transits:
└ dr3 = 3
┌ Info: DR2/DR3 epoch selection
│ n2_win = 13
│ n_tail = 14
│ n_dr2_total = 13
└ n_dr2_distinct_range = (13, 13)
[ Info: [limits_both] observing_geometry = true (auto): no observation can report comparable predictions (ImageObs, G23HObs); nothing to measure, so keeping the correction on — no draws needed (seed 0xc70f177e5000001)
[ Info: [limits_both] barycentric_lighttime = true (auto): no observation can report comparable predictions (ImageObs, G23HObs); nothing to measure, so keeping the correction on — no draws needed (seed 0xc70f177e5000001)
[ Info: Preparing model
┌ Info: Determined number of free variables
└ D = 11
┌ Info: Determined number type
└ T = Float64
ℓπcallback(θ): 0.000068 seconds (42 allocations: 72.547 KiB)
∇ℓπcallback(θ): 7.091796 seconds (25.10 M allocations: 1.063 GiB, 1.04% gc time, 99.99% compilation time)
┌ Info: Starting values not provided for all parameters! Guessing starting point using global optimization:
│ num_params = 11
└ num_fixed = 0
┌ Info: Found sample of initial positions
│ logpost_range = (-6.171135183887338, 5.591735188162785)
└ mean_logpost = 1.9068800102976986
[ Info: Sampler running with multiple threads : true
[ Info: Likelihood evaluated with multiple threads: false
[ Info: [limits_both] observing_geometry = true (user)
[ Info: [limits_both] barycentric_lighttime = true (user)
[ Info: [limits_both] observing_geometry = true (user)
[ Info: [limits_both] barycentric_lighttime = true (user)
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
scans restarts Λ Λ_var time(s) allc(B) log(Z₁/Z₀) min(α) mean(α) min(αₑ) mean(αₑ)
────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ────────── ──────────
2 0 2.27 2.22 1.6 4.33e+08 -1.27e+06 0 0.501 1 1
4 0 2.37 2.23 0.399 8.08e+08 -2.75e+05 0 0.489 1 1
8 0 2.9 3.51 0.673 1.58e+09 -5.34e+04 0 0.287 1 1
16 0 3.11 3.5 1.39 3.21e+09 -1.06e+04 0 0.265 1 1
32 0 3.51 3.65 2.71 6.43e+09 -3.5e+03 0 0.204 1 1
64 0 3.61 2.47 6.22 1.34e+10 -8.22 0 0.325 1 1
128 1 3.7 3.19 12 2.83e+10 -5.57 0 0.235 1 1
256 8 3.67 2.78 24.3 5.7e+10 -4.34 0 0.283 1 1
512 23 3.73 2.58 48.8 1.13e+11 -4.16 0 0.299 1 1
1.02e+03 54 3.78 2.63 92.4 2.24e+11 -4.06 7.12e-72 0.287 1 1
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────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][:] ./ mjup),
),
label="PMA",
color=Makie.wong_colors()[1],
)=>vis_layers,
PairPlots.Series(
(;
sma=log.(chain_img[:B_a][:],),
mass=log.(chain_img[:B_mass][:] ./ mjup),
),
label="IMG",
color=Makie.wong_colors()[2],
)=>vis_layers,
PairPlots.Series(
(;
sma=log.(chain_both[:B_a][:],),
mass=log.(chain_both[:B_mass][:] ./ mjup),
),
label="IMG + PMA",
color=Makie.wong_colors()[3],
)=>vis_layers,
labels=Dict(
:sma=>"log Semi-major axis [au]",
:mass=>"log Mass [Mⱼᵤₚ]"
)
)
For a systematic version of this analysis — injecting companions on a grid of mass and separation and measuring what fraction are recovered — see Detection Completeness Mapping.