Detection Limits

Warning

This tutorial is a work in progress.

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 PairPlots

Photometry 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#193"{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:

@planet B Visual{KepOrbit} 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)# ~ UniformCircular()
    mass = system.M_sec
    θ ~ Uniform(0,2pi)
    tp = θ_at_epoch_to_tperi(system,B,57423.0) # epoch of GAIA measurement
end
@system HD91312_pma begin
    M_pri ~ truncated(Normal(0.95, 0.05), lower=0) # Msol
    M_sec ~ LogUniform(0.2, 65) # MJup
    M = system.M_pri + system.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 HGCALikelihood(gaia_id=6166183842771027328) B
model_pma = Octofitter.LogDensityModel(HD91312_pma)

LogDensityModel for System HD91312_pma of dimension 11 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: Determining initial positions and metric using pathfinder
┌ Info: Found a sample of initial positions
└   initial_logpost_range = (-19.066848533790345, -10.883460749011629)
┌ Warning: Invalid log likelihood encountered. (maxlog=1)
│   θ = (M_pri = 0.0017590160075691803, M_sec = 1.1631153538983048, plx = 1.545265778533215e218, pmra = -207.91839355444154, pmdec = -60.30991115391226, M = 0.0028689769897943325, planets = (B = (a = 1.5911088543388132, e = 0.81490014365075, ω = 5.764081435809116, i = 1.9366331422043717, Ω = 0.3461842616738291, θ = 2.031442101059195, mass = 1.1631153538983048, tp = 43908.5603026686),))
│   llike = -Inf
│   θ_transformed =
│    11-element Vector{Float64}:
│     -6.343000712887692
│      ⋮
└ @ 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.25       2.48       5.27    3.4e+08  -1.21e+08          0      0.815      0.926      0.963
        4          0       5.12       4.15      0.481   9.03e+06  -1.55e+05          0      0.701      0.891      0.923
        8          0       5.36        4.6      0.524   9.14e+05  -4.61e+05          0      0.679      0.877      0.915
       16          0       7.29       5.93       1.09   1.73e+06       -228          0      0.573      0.886      0.904
       32          0       7.64       6.67       2.18    3.2e+06       -349   1.7e-298      0.538      0.891      0.907
       64          3       8.08       2.62       5.49   2.28e+08      -17.5   8.81e-27      0.655      0.895      0.906
      128          9       8.42       3.51       9.47   4.09e+08      -17.2   5.43e-06      0.615        0.9       0.91
      256         17       8.93       3.44       18.4   7.86e+08      -17.1   5.21e-05      0.601        0.9       0.91
      512         39       9.26       3.32       37.3   1.59e+09      -17.1     0.0187      0.594      0.902      0.912
 1.02e+03         91       9.43       3.47       74.8   3.18e+09      -17.3     0.0478      0.584      0.903      0.912
 2.05e+03        249       9.82       2.87        150   6.37e+09      -17.2      0.134      0.591      0.901      0.908
  4.1e+03        444       9.93       3.36        302   1.28e+10      -17.3      0.223      0.571      0.902       0.91
────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

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

image_data = ImageLikelihood(
    (band=:L, image=AstroImages.recenter(image), platescale=4.0, epoch=57423.6),
)

OctofitterImages.ImageLikelihood Table with 6 columns and 1 row:
     band  image                 platescale  epoch    contrast              ⋯
   ┌─────────────────────────────────────────────────────────────────────────
 1 │ L     [-6.80461e-8 -6.247…  4.0         57423.6  65-element extrapol…  ⋯

@planet B Visual{KepOrbit} 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, B.mass)
    # Calculate absolute magnitude
    abs_mag_L = sonora_temp_mass_L(B.tempK, B.mass)
    # Deal with out-of-grid values by clamping to grid max and min
    abs_mal_L′ = if isfinite(B.abs_mag_L)
        B.abs_mag_L
    elseif B.mass > 10
        8.2 # jump to absurdly bright
    else
        16.7 # jump to absurdly dim
    end
    # Calculate relative magnitude
    rel_mag_L = B.abs_mal_L′ - system.rel_mag + 5log10(1000/system.plx)
    # Convert to contrast (same units as image)
    L = 10.0^(B.rel_mag_L/-2.5)

    θ ~ Uniform(0,2pi)
    tp = θ_at_epoch_to_tperi(system,B,57423.6)
end image_data

@system HD91312_img begin
    # age ~ truncated(Normal(40, 15),lower=0, upper=200)
    age = 10
    M_pri ~ truncated(Normal(0.95, 0.05), lower=0) # 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 = system.M_pri + system.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  B
model_img = Octofitter.LogDensityModel(HD91312_img)

LogDensityModel for System HD91312_img of dimension 9 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×21×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

@planet B Visual{KepOrbit} 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, B.mass)
    # Calculate absolute magnitude
    abs_mag_L = sonora_temp_mass_L(B.tempK, B.mass)
    # Deal with out-of-grid values by clamping to grid max and min
    abs_mal_L′ = if isfinite(B.abs_mag_L)
        B.abs_mag_L
    elseif B.mass > 10
        8.2 # jump to absurdly bright
    else
        16.7 # jump to absurdly dim
    end
    # Calculate relative magnitude
    rel_mag_L = B.abs_mal_L′ - system.rel_mag + 5log10(1000/system.plx)
    # Convert to contrast (same units as image)
    L = 10.0^(B.rel_mag_L/-2.5)

    # L ~ Uniform(0,1)

    θ ~ Uniform(0,2pi)
    tp = θ_at_epoch_to_tperi(system,B,57423.6)
end image_data

@system HD91312_both begin
    # age ~ truncated(Normal(40, 15),lower=0, upper=200)
    age = 10
    M_pri ~ truncated(Normal(0.95, 0.05), lower=0) # 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 = system.M_pri + system.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 HGCALikelihood(gaia_id=6166183842771027328) B
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×23×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ⱼᵤₚ]"
    )
)