Octofitter: Fast, Flexible, and Accurate Orbit Modeling to Detect Exoplanets

As next-generation imaging instruments and interferometers search for planets closer to their stars, they must contend with increasing orbital motion and longer integration times. These compounding effects make it difficult to detect faint planets but also present an opportunity. Increased orbital motion makes it possible to move the search for planets into the orbital domain, where direct images can be freely combined with the radial velocity and proper motion anomaly, even without a confirmed detection in any single epoch. In this paper, we present a fast and differentiable multimethod orbit-modeling and planet detection code called Octofitter. This code is designed to be highly modular and allows users to easily adjust priors, change parameterizations, and specify arbitrary function relations between the parameters of one or more planets. Octofitter further supplies tools for examining model outputs including prior and posterior predictive checks and simulation-based calibration. We demonstrate the capabilities of Octofitter on real and simulated data from different instruments and methods, including HD 91312, simulated JWST/NIRISS aperture masking interferometry observations, radial velocity curves, and grids of images from the Gemini Planet Imager. We show that Octofitter can reliably recover faint planets in long sequences of images with arbitrary orbital motion. This publicly available tool will enable the broad application of multiepoch and multimethod exoplanet detection, which could improve how future targeted ground- and space-based surveys are performed. Finally, its rapid convergence makes it a useful addition to the existing ecosystem of tools for modeling the orbits of directly imaged planets.


INTRODUCTION * 51 Pegasi b Fellow
Instruments for directly studying exoplanets are steadily improving in sensitivity.Current facilities are now accessing planets less than 10 AU from their stars.Below these separations, orbital motion can become sig-nificant over mere months.This will be especially true for facilities with high angular resolving power thanks to their larger apertures and/or shorter operating wavelengths.This is an advantage for those who wish to determine the orbits of already-known planets but presents a significant challenge when searching for new companions.
When planets move from observation to observation, naive image stacking causes their signals to blur out and fade away.Reaching a necessary integration in a single epoch hits practical scheduling constraints and eventually physical limitations-for a sufficiently faint planet, it would not be possible to detect a significant number of photons before it moves by a full resolution element.These constraints apply equally to images as they do to integral field units and interferometers, including aperture masking interferometry (AMI) on JWST (Sivaramakrishnan et al. 2023), VLTI-GRAVITY (Collaboration et al. 2017), and even ALMA (Wootten & Thompson 2009) in the case of accreting protoplanets.
A number of projects have sought to solve this challenge for image data only by compensating for orbital motion between epochs.These include a search for planets around Sirius B (Skemer & Close 2011), K-Stacker (Nowak et al. 2018;Le Coroller et al. 2020), PACOME (Dallant et al., submitted), the search for planets around ε Eri (Mawet et al. 2019;Llop-Sayson et al. 2021), and the search for additional HR 8799 planets by Thompson et al. (2023).These have now led to promising evidence for α Cen AB b (Le Coroller et al. 2022) and HR 8799 f (Thompson et al. 2023).Moving the analysis of direct imaging data into the orbital domain enables a further extension: joint models of both images or interferometric observables with indirect exoplanet detection techniques, including radial velocity, astrometric motion, and transit (not directly considered in this paper).These have previously been explored in Mawet et al. (2019) and Llop-Sayson et al. (2021).This opens many possible scenarios.In addition to combining images to search for planets despite orbital motion, this allows one to freely combine Doppler or astrometric velocimetry with image data.This can then be used to constrain orbits using images with tentative detections or non-detections, improve photometry accuracy or limits, better constrain a planet's mass, and/or detect planets where any individual kind of data fails to reach significance (e.g.Llop-Sayson et al. 2021).
These scenarios are possible since all exoplanet detection methods provide orbital constraints that at least partially overlap.Imaging, RV, and transit all provide the orbital period (P ); proper motion anomaly (PMA) and RV constrain eccentricity (e), the argument of pe-riapsis (ω), and either mass (m) or m sin(i) where i is the orbital inclination; and multi-epoch imaging/interferometry constrains all orbital parameters up to a ± ambiguity on the longitude of the ascending node (Ω).These connections could, in principle, allow information from all methods to flow into a single orbit model and, ultimately, the detection of a new planet.
To apply these ideas broadly, the community will need a tool capable of modelling all different types of exoplanet data.The orbitize!(Blunt et al. 2020) and orvara (Brandt, T. D. et al. 2021) packages come close: they support Bayesian modelling of relative astrometry, RV, and PMA.We needed a publicly available and generally applicable package that goes further to directly model image and interferometer data with or without independent detections at each epoch.
It is, generally, challenging to accurately compute orbital posteriors since the traditional Campbell orbital elements (a, e, i, ω, Ω, t peri ) possess complex codependencies and degeneracies (e.g. when e = 0 or i = 0).This task becomes even more challenging when working with short orbital arcs because they lack the constraining power to independently determine each Campbell element, meaning orbit posteriors are typically complex and very sensitive to their priors (e.g.O'Neil et al. 2019).Introducing image and interferometer data to the model exacerbates issues further, as they produce multi-modal posteriors that are challenging to traverse.Any inaccuracies in the calculation of an orbit posterior can lead to errors in mass and/or photometry and a spurious detection.
These challenges motivate the development of our new orbit-and data-modelling framework, "Octofitter."1Named for the eight types of data through which it aims to grasp new planets, Octofitter is designed from the ground up to be a flexible platform for modelling and experimentation while providing very high computational performance.These advances are thanks to our implementation of a pure-Julia (Bezanson et al. 2012), fully differentiable, and non-allocating modelling language, and our use of the higher-order No-U-Turn Hamiltonian Monte Carlo sampler (Xu et al. 2020).We further present the use of Simulation-Based Calibration (SBC; Talts et al. 2020) to confirm the accuracy of orbital posteriors, which is made practical by Octofitter's speed.Figure 1 shows a schematic of a few of the ways Octofitter can be applied to exoplanet data.
As a publicly available code with these advances, Octofitter will enable the wide application of multiepoch, multi-instrument, and multi-method exoplanet detection and modelling.These approaches have the ability to improve how direct surveys are completed and improve the yield of the upcoming Nancy Grace Roman Space Telescope Coronagraphic Instrument (Kasdin et al. 2020), the Habitable Worlds Observatory (HWO), and future facilities for extremely large telescopes.

DATA MODELS
This section describes how the eight kinds of exoplanet data can be modelled by Octofitter.These types of data are relative astrometry, images, interferometric observables, star and planet radial velocity (absolute and relative), and proper motion anomaly.
The Octofitter framework is structured around three concepts: likelihood functions for observations, system models to tie observations to parameters, and generative functions to create synthetic observations.These break the problem down into orthogonal components that can be freely combined to solve a wide range of orbit modelling problems.
Each kind of observation in Octofitter is supported by its own Julia data type.Every observation type is a wrapper for a data table of observations with a preset list of required columns.The observational data, in turn, can be provided to its associated type directly in the code, loaded from a local CSV or Arrow file, or obtained from a remote SQL database.For each observation type, a method of the Octofitter.ln_likefunction is provided that computes the likelihood of the data it contains given a specific set of parameters.We now describe the observation types and likelihood functions included in Octofitter.

Relative Astrometry
The position measured between a directly imaged planet and its host star is one of the most fundamental measurements gathered from direct observations.It can be extracted from any image where the planet is robustly detected in a single epoch.Relative astrometry can be expressed either as separation in milliarcseconds and position angle in degrees, or as offsets in milliarcseconds in the R.A.-Dec.tangent plane.
Relative astrometry measurements can be provided to Octofitter using AstrometryLikelihood, which accepts a table with columns for epoch, the date of the measurement in units of modified Julian days; sep, the separation between primary and secondary in mas; pa, the position angle of the secondary measured East from North;  1.
Based on these quantities, the likelihood function is simply taken as the Gaussian likelihood that the residual between the model and measurement for each parameter would be seen, given the provided measurement uncertainty.For relative positions measured along the R.A. and Dec. axes, this is  where ∆RA and ∆DEC are the separation from the star along the R.A. and Dec. axes, σ ∆RA and σ DEC are the corresponding uncertainties, and where the tilde distinguishes measured quantities from values calculated from model parameters.A more complex expression is used when the correlations between uncertainties are nonzero.

Images
Directly modelling point sources in images is one of the key features of Octofitter.This is accomplished with the same approach described in Ruffio et al. (2018) and applied in Mawet et al. (2019), Llop-Sayson et al. (2021), and Thompson et al. (2023).In the case where a planet is robustly detected in each image and the uncertainties are well-approximated by Gaussian uncertainties, it is mathematically equivalent to extracting relative astrometry and photometry and then modelling those measurements.Where this approach goes beyond the two-step process, however, is when a planet is not detectable in a single epoch or when one wants to constrain the orbit of a planet based on a non-detection.The first case may arise when searching for planets that are too faint to detect before they exhibit orbital motion.The second case may arise any time a monitored planet passes too close to its star to detect.In this case, it's not possible to extract astrometry.By directly modelling the images, large swathes of the orbital parameter space may still be ruled out.This can improve constraints on orbital parameters and improve predictions of the planet's future location.
In Octofitter, image data can be modelled across any number of photometric bands, instruments, and epochs.Data can be provided using the ImageLikelihood observation type which accepts a data table with the columns epoch, band, platescale, image, and contrast.A sample input table is presented in Table 2.The platescale and parallax system parameters (which may be fixed or fit to the data) are used to map an orbit to a pixel location in the image.The contrast entry for each image allows one to pass an arbitrary function of position that gives 1σ contrast.If not provided, Octofitter calculates a contrast curve automatically from the image itself.This is sufficient for images with no clear detection, but should be avoided if the image contains a bright planet as the contrast curve will overestimate the noise at the planet's separation.Extended emission from disks is not currently considered.
Given this data, the likelihood function used by Octofitter for each image by Octofitter is that of Ruffio et al. (2018): where f B is the model flux parameter for the photometric band B, x is the position in the image determined from orbital parameters, i is the epoch, σ is the uncertainty, and fB is the measured flux at that same location.This likelihood function assumes that flux is constant from epoch to epoch, but could easily be adapted to use an orbital phase function for planets imaged in reflected light (e.g.Pogorelyuk et al. 2022).Taken to the extreme, the flux can be fit independently at each epoch (e.g.Nowak et al. 2018) at the expense of reduced constraining power.For both images and relative astrometry, we do not currently consider uncertainty in the instrument's North-angle or platescale.These might need to be considered when combining data from multiple instruments, in which case they could be added straightforwardly.We also note that this likelihood function assumes Gaussian distributed noise in each image.Where this is not a good assumption, standard techniques for inflating uncertainties could be used before applying Octofitter.Furthermore, spatial correlations between pixels are not currently handled.Finally, this likelihood function is only specified up to a constant and is not suitable for techniques like nested sampling.

Interferometric Observables
Just as with imaging, combining interferometric observations across multiple epochs using orbital modelling removes the need to detect a companion in a single observation.We implement a model assuming an unresolved point source primary and N unresolved point source companions.The complex visibilities of this model are given by where f i is the (companion/primary) contrast of the ith companion, ∆RA i and ∆DEC i are the right ascension and declination of the ith companion, and u and v are the Fourier domain coordinates, with magnitudes given by the interferometer baseline lengths divided by the observing wavelength (e.g., Berger 2003;Kammerer et al. 2023).A sample input table is presented in Table 3.
The squared visibilities are calculated from the squared moduli of the complex visibilities, and the closure phases are calculated by summing the phases of the three complex visibilities calculated from triangles of stations in the interferometer.We construct likelihood functions for the squared visibilities and the closure phases separately, assuming Gaussian noise statistics and diagonal covariance matrices.We also assume that there is at least a moderate contrast between the primary and any companions such that we can neglect phase wrapping in the closure phase (Éric Thiébaut & Young 2017).These likelihood functions are given by and where CP i is the ith closure phase, V 2 i is the ith squared visibility and σ CP,i , and σ V 2 ,i are the uncertainties in the ith squared visibility and closure phase, respectively.

Radial Velocity
Similarly to other packages, Octofitter allows one to model radial velocity data in combination with other observation types.If combined with relative astrometry, proper motion anomaly, or image data, this allows one to directly model the dynamical mass of a planet.
RV data can be specified for the star or the planet using RadialVelocityLikelihood, which accepts a table with columns for epoch in MJD, rv in m/s, σ_rv, the uncertainty on rv in the same units, and inst_idx.inst_idx is an integer between 1 and 4 used to specify which instrument the measurement corresponds to.The zero point rv0_i and jitter jitter_i must be specified as variables in the model where i corresponds to the instrument index.
A sample input table is presented in Table 4.The input format is the same for both cases.Depending on how the zero point is modelled, it is possible to use either relative or absolute RVs.The zero point can also modelled as an arbitrary function of other variables, allowing one to fit, for example, linear trends.
When combining radial velocity data with direct imaging modelling, it is possible (though not required) to connect the planet's photometry with its dynamical mass by using a user-supplied model.This model can either map a mass and/or system age to photometry in each band, or vice-versa.Connecting these two variables may be useful in cases such as when the orbit is determined by radial velocity up to the inclination degeneracy but has not yet been detected with direct imaging.
In a similar manner to relative astrometry, we define a radial velocity likelihood function that allows us to fit orbital parameters to radial velocity data.This function is where v r,i is the measured radial velocity, ṽr,i is the maximum likelihood estimate of the radial velocity, and σ vr,i is the uncertainty in the radial velocity, all at the epoch i.
Octofitter supports multiple instruments, each with their own RV zero-point and jitter term and all with independently selectable priors.This is accomplished via the inst_idx column which associates each RV measurement to a particular instrument, zero point, and jitter.As a convenience, Octofitter includes a helper function to load RV curves from the publicly available HARPS RV Bank (Trifonov et al. 2020).Accessing this data will prompt the user to accept its license and will automatically fetch the RV curves.

Proper Motion Anomaly
Many systems that are candidates for direct imaging have their positions and proper motions measured ac-  In Octofitter, proper motion anomaly must be loaded directly from the HGCA by specifying the host star's Gaia identifier (the DR3 version at the time of writing).As with the radial velocity loaders, this will prompt the user to accept a license for the catalog and automatically download the HGCA.
In order to connect proper motion anomaly with the orbital image modelling described above, we define a likelihood function based on the HGCA data and the same orbital parameters: where v RA and ṽDEC are the R.  2021) we adopt a simpler model of the Gaia and Hipparcos data.Rather than use the exact epochs each mission scanned the star (or estimate using the Gaia Observation Forecast Tool 2 ), we presume that the position and proper motion of the star were measured independently 25 times over the duration of each mission.For orbits with periods longer than the observation windows of Gaia (≈ 668 days) and Hipparcos (≈ 1227 days), using these 25 epochs approximates the smearing effect caused by the moving photocenter during observations.Future work could refine this model further following their more sophisticated approach by re-implementing the calculations of HTOF (Brandt, G. M. et al. 2021) in Julia.Once it is available, it should also be possible to use the full intermediate GAIA astrometry catalogue in Octofitter.

Assessing Detections
This section describes how we choose to interpret the results of a model fit with Octofitter, though it should be noted that the modelling code makes no such prescriptions.For evaluating upper limits and detections from image or visibility data alone, we follow the conventions of Ruffio et al. (2018).
As set out in Ruffio et al. (2018), we calculate upper photometry or mass limits by finding where the cumulative distribution function equals some threshold e.g.97.7%.That is, CDF (f lim ) = 97.7%.This can be stated as an overall value for the system to say, for example, that we believe with 97.7% confidence that there are no planets present at all with photometry above some f lim .It can also be calculated over small ranges of orbital parameters, in order to, for example, calculate a flux upper limit as a function of semi-major axis.This value will be driven by the brightest speckles close to the star, so a 2 https://gaia.esac.esa.int/gost/much more useful metric is to provide f lim as a function of orbital parameters.
This model does presuppose that a planet exists and follows a Keplerian orbit around the star; however, as long as the photometry prior is broad, the model is free to drive the photometry towards zero.We can therefore use the photometry posterior to calculate a signal to noise ratio (SNR) analogous to the typical metric used in direct imaging.This approach has the benefit of being familiar to those with experience evaluating direct images for detections.
The procedure we adopt for assessing detections is to 1. sample from the posterior, 2. marginalize over all orbital parameters, 3. calculate the SNR as the mean of the marginal photometry posterior divided by its standard deviation, 4. and then compare this SNR with a pre-selected threshold based on a tolerable false positive rate, e.g.5σ.
The results of this procedure are visualized in Figure 2.This gives a single SNR value for existence of any planet in that dataset consistent with the user's choice of priors.This is somewhat different than the typical SNR detection thresholds used because this SNR is calculated by marginalizing over all credible orbits, or equivalently, if applied to a single image, marginalizing over all credible positions.By contrast, the standard 5σ threshold is applied to the position in an image with the highest SNR, rather than a weighted combination of all positions.
In general, SNRs obtained using this approach assume that 1. the noise is approximately Gaussian distributed, 2. the noise is uncorrelated along orbital tracks between images, 3. the images contain only a single planet, 4. and that planet closely follows a Keplerian orbit.
Deviations from assumptions (1-2) may increase the false positive rate, while deviations from (3-4) may lower the recovered photometry and SNR.
Multi-planet systems could be accommodated easily by introducing additional planets to the model and/or restricting the priors to include only a subset of the parameter space.
For models that combine direct and indirect data, the relationship between mass and photometry variables Posterior density of photometry versus semimajor axis for a simulated detection and non-detection.The solid lines mark the mean of the marginal photometry posterior, and the dashed lines mark ±1σ.We consider a planet detected when the SNR based on the photometry marginalized over all other parameters is greater than some chosen threshold, e.g.5σ.
may be complex, and this simple scheme may not be appropriate.In this case, it may be better to calculate a Bayes factor between planet and no-planet models.Such a Bayes factor can be treated as a direct measurement of our belief that a planet exists and follows a Keplerian orbit.This can be carried out in a limited fashion in Octofitter by calculating the Savage-Dickey density ratio of the mass or photometry marginal posterior (Dickey (1971); Koop (2003); for a pedagogical text, see Wagenmakers et al. (2010)).This is more flexible since it does not require a uniform prior-any prior that includes zero mass or photometry would suffice-and because it does not presume that the marginal posterior is Gaussian distributed.

Orbital Bases and Priors
Octofitter supports a wide range of different orbital bases for use in different situations.These include traditional Campbell elements (a, e, i, ω, Ω), Thiele-Innes elements (e, A, B, F, G, Hartkopf et al. 1989;Wright & Howard 2009;O'Neil et al. 2019), Cartesian elements (x, y, z, v x , v y , v z , Ferrer-Chávez et al. 2021a), and a reduced basis set for modelling radial velocity only (a, e, m sin(i), ω).Users can specify priors using arbitrary distributions from Distributions.jl 3 and functions of those distributions.
For the analyses presented in this paper, we adopt either Campbell elements or Thiele-Innes elements with the following priors and modifications.When using 3 https://juliastats.org/Distributions.jlCampbell elements, we adopt a log-uniform prior on semi-major axis, a uniform prior on eccentricity, a sine prior on inclination, a Gaussian prior on host mass, and uniform priors on the remaining angular parameters.Note Octofitter reports the argument of periastron of the planet (and not the star) as ω and adopts +z increasing away from the observer.These conventions match those used by orvara and orbitize!(Householder & Weiss 2023).The full conventions used in Octofitter are described in Appendix A. When using Thiele-Innes elements, we adopt a log uniform prior on a "scale" parameter and multiply this with standard Gaussian priors on the constants A, B, F, and G.This maintains a log-uniform prior on semi-major axis and sine prior on inclination.For both cases, we replace the parameter τ , which gives the position of a planet along its orbit (Blunt et al. 2020) , with θ, the position angle at the average epoch of the observations.This improves convergence relative to sampling from the τ parameter directly since θ is directly constrained by relative astrometry and is not sensitive to other orbital parameters.A derivation of τ from position angle θ is presented in Appendix B.

Modelling Language
To specify the structure of a system model, Octofitter provides a domain-specific modelling language.This simple language allows for the parameterizations and observations associated with each planet and the host star to be independently specified.For example, we could attach separate radial velocity measurements to each object in a system.Orbital parameters can be drawn from arbitrary prior distributions, fixed to particular values, or computed from combinations of other system parameters.
This modelling language makes it convenient to work with simple systems, like fitting the orbit of one planet to relative astrometry measurements, as well as more complex multi-planet systems.
The following is an example of a planet model using traditional Campbell orbital parameters, user-specified priors, and relative astrometry data: § ¤

¦ ¥
The planet and system definition blocks contain pairs of variable names and values which can be constants, prior distributions, or arbitrary functions of other variables.Variables drawn from priors are specified by whereas variables defined as constants or calculated as a function of other variables are defined with =.In this example, the convenience function UniformCircular creates two independent variables with standard normal priors and computes the angle between them using the arctangent.This, in effect, creates a circular domain where samples can smoothly wrap past 0 and 2π.Users are free to create their own parameterizations and likelihoods and combine them arbitrarily.The only requirements are that they are differentiable, smooth, and return a finite value at all points in the domain given by the priors.
This modelling approach, by being declarative rather than imperative, as in exoplanet.py(Foreman-Mackey et al. 2021), allows us to transform and evaluate the model in several ways.One key restriction is that each prior is proper, meaning it is a true probability distribution that integrates to unity.This is in contrast to some of the defaults used by orvara, which adopts fully scaleindependent, but improper (in the statistical sense of the word) log-uniform priors.Octofitter requires proper priors to support tasks that require sampling directly from the prior distributions, such as simulation-based calibration, as will be discussed in Section 3.7.

Numerical Methods
From a model definition, Octofitter can generate efficient machine code using runtime-generated functions and Julia's just-in-time compiler.This code generation step allows Octofitter to support a rich variety of observation types without paying any runtime overhead for features that aren't used.
For the purposes of sampling from the posterior, Octofitter begins by remapping all variables from their possibly limited support (for example, eccentricity constrained between 0 and 1) into new variables defined across the entire real line.This makes it so that by construction, invalid parameter values like negative masses or semi-major axes are not possible to construct and will not be hit by a sampler.This transformation is performed by the Julia package Bijectors.jl4 .
Next, a log-prior function is created for the model that simply evaluates the log-probability density of each prior distribution given a set of parameters.To preserve the user-specified prior distributions in place of the automatic bijections, a correction is applied based on the Jacobian of the transformation.
Similarly, a log-likelihood function is created based on the provided model and observations.Various constants are pre-calculated and reused between orbit solutions.
To enable the use of higher-order samplers, including Hamiltonian Monte Carlo, forward mode automatic differentiation is used to differentiate through the generated log-prior and log-likelihood functions (Revels et al. 2016).This provides the gradient of the log-posterior density, in addition to the value itself, without the overhead of calculating finite differences.The overhead of calculating both the log-posterior density and its gradient using forward mode automatic differentiation can be as low as just 2× compared to 10.2× for finite differences.
Special care was taken to remove all dynamic memory allocations from the generated log-density and gradient functions to prevent overhead from the Julia garbage collector.
Octofitter implements the Julia LogDensityProblems interface so that user models can be sampled from a wide variety of Julia-based MCMC sampler packages, including AdvancedMH, AdvancedHMC, and MCMCTempering.This allows users to select the best sampler for their particular problem and to compare results against samplers used by other popular packages.
The No-U-Turn Sampler (NUTS) variant of Hamiltonian Monte Carlo (Hoffman & Gelman 2014) provided by AdvancedHMC is the default used by Octofitter.It allows the code to explore complex posterior distributions with many fewer log-posterior density evaluations by simulating physical trajectories across the posterior landscape.The performance of NUTS in Octofitter will be discussed in Section 4.1.
When using the default NUTS sampler, we use a dense mass matrix, a jittered leapfrog integrator, a multinomial trajectory sampler, and allow the user to specify a maximum tree depth.We initialize the sampler by drawing 250,000 samples from the priors and selecting the value with the highest posterior density as the starting point.We also initialize the diagonal elements of the mass matrix using the interquartile range of the priors.
We found this procedure was more robust than trying to determine the maximum a posteriori value with an optimizer for multi-modal posteriors like those found when modelling images.In particular, this procedure is less likely than an optimization pass to get stuck in a local optimum or other pathological position.A downside of this approach is that it may be inefficient if the priors are narrow and far from the posterior.For this case, the user can supply a starting point manually as in other tools.Next, we adapt the mass matrix and step size according to AdvancedHMC's implementation of the Stan windowed adaptation strategy.Finally, sampling proceeds until a preset number of accepted proposals are found-typically on the order of 1,000 to 10,000.

Kepler Solver
To sample from planet models, one must map from the orbital parameters specified in the system model to simulated observations, such as a planet's position over time in the plane of the sky or the radial velocity of the host star.In all cases, this requires solving Kepler's equation at every observation epoch and parameter draw.Kepler's equation (Eq.A3) connects eccentricity and mean anomaly, a pseudo-angle somewhat analogous to the amount a planet has moved around its orbit, into eccentric anomaly, which can be used to find the actual location of the planet in its orbit.
Kepler's equation is transcendental and can't be solved analytically outside of special cases.The traditional approach to solving the equation is to use an iterative procedure like Newton's method with an initial guess chosen carefully from the mean anomaly and eccentricity.Octofitter's strategy is to wrap many different pluggable Kepler solvers that can be useful in different scenarios and use a fast non-iterative method as a default.
A non-exhaustive list of solvers supported by Octofitter are Markley (Markley 1995), Goat (Philcox et al. 2021), and Newton.Newton's method and many other available root finding algorithms are provided by the Roots.jl5julia package.
Since these solvers are all implemented in pure Julia, there is no overhead from calling between Python and a C/Cython6 library and full performance is achieved with or without vectorization.Additionally, the solvers support changing the numerical precision between 16bit, 32-bit, 64-bit, and arbitrary precision.Currently, no particular effort has been made to exploit hardware SIMD vectorization across epochs.
The Markley algorithm, the default choice, converges to nearly machine precision for all bound orbits when used with 64-bit floating point values.The Newton method can be combined with arbitrary precision floating point numbers to achieve arbitrarily tight tolerances if needed for niche applications.The default Markley 1995 based method executes in just 90 ns on 64-bit floating point values for any valid eccentricity and mean anomaly (benchmarked on a Skylake Intel Xeon processor).This is slightly slower than the state-of-the-art results reported by Brandt et al. (2021); however, their results of ≈ 60 ns are only achieved when the solver is vectorized over many epochs.
Applying automatic differentiation through an iterative Kepler equation solver would lead to poor performance.Even though Kepler's equation has no closed solution, the derivatives of eccentric anomaly with respect to mean anomaly and eccentricity can be found analytically using implicit differentiation (Eqs.A5 and A6); that is, if we have already solved Kepler's equation for eccentric anomaly, we can calculate its gradient inexpensively.We supply these equations as a manual rule to the automatic differentiation library.

Analysis and Visualization
Once sampling is complete, Octofitter supports the user in testing the convergence of their chain or chains.The sampling results are returned as Chains objects from MCMCChains.jl.This table-like structure includes entries for each accepted MCMC proposal as well as metadata about the sampling process, such as the compute time used.All variables are returned, including those that were fixed to a constant or calculated deterministically from combinations of other variables.An automatic summary is output that includes plausible intervals, expectation values, effective sample sizes (ESS), and R statistics for each parameter.The user can test the convergence of their chain or chains by assessing those ESS and R values, creating a trace plot, and/or by calculating figures like the Gelman, Rubin, and Brooks diagnostic (Gelman & Rubin 1992;Brooks & Gelman 1998) using tools provided by MCMCChains.jl and MCMCDiagnosticTools.jl.
Octofitter can be used to visualize orbit fits in several ways.The orbits represented by chains can be plotted in the plane of the sky, in physical system coordinates (i.e.AU), or as a time series (e.g. for Doppler or astrometric velocimetry).
When visualizing an orbital posterior, a common challenge is ensuring enough data points are used to create a smooth arc.This becomes especially challenging with eccentric orbits or ensembles of orbits with widely vary- Sample results of running the simulation-based calibration procedure on a model consisting of a single planet parameterized with Thiele-Innes elements (A, B, F , G, etc.) and position angle θ.Each count represents a model fit to a different simulated system.The horizontal band gives a 99% range around the expected value for a perfect sampling procedure.The observation epochs and uncertainties are taken from Table 1.
ing periods.Tracing out orbits in equal time steps will lead to most points clustering closer to apoastron, where the planet is moving the most slowly.The strategy used by Octofitter is to trace orbits in equal steps of eccentric anomaly, as suggested in Berry & Healy (2002).This places points in regions of greater curvature, creating smooth arcs with fewer points.
Octofitter includes functions to generate plots that visualize the posterior and compare it to the input data.This includes astrometry, separation, position angle, radial velocity, and proper motion anomaly.It also includes a function for visualizing orbits in spatial coordinates (units of AU) in one, two, or three dimensions.Examples of these plots are shown throughout the text and documented with the online tutorials.

Simulation-Based Calibration
Simulation-based calibration (SBC; Cook et al. 2006;Talts et al. 2020) is a technique that allows one to verify that the output of a Bayesian modelling procedure is unbiased.This includes any part of the computation, from model specification to the sampling procedure after the choice of priors.Verifying that this choice of priors is reasonable is the domain of other procedures like prior and posterior predictive tests.
A conceptually related procedure was carried out in Ferrer-Chávez et al. (2021b) in which the authors systematically tested the orbitize!package, parameterization, and default priors for biases.Our contribution here is to present a procedure that is well explored in the statistics literature, is specific to a given set of observation epochs and measurement uncertainties, and can be applied in an automated fashion.
The calibration procedure requires a generative model-that is, a way to take a given set of parameters and create a simulated observation.A familiar example from direct imaging is the injection of fake point sources to test image processing algorithms.To follow this procedure, one repeatedly draws a set of parameters from the priors, creates simulated observations based on those parameters, samples from the resulting posterior, and then compares the true parameter values (originally sampled from the priors) with the resulting posterior.By doing this many times, one creates a histogram of rank statistics that can reveal many sources of biases present in the model and sampling process.To be clear, this does not evaluate how the choice of priors impacts the posterior.
We implemented support for performing simulationbased calibration automatically in Octofitter.In this procedure, Octofitter takes a given model's observational data, discards the actual measurements, and keeps the epochs and uncertainties associated with each measurement.Octofitter then repeatedly creates simulated data at each epoch by drawing from the priors and performs SBC on this simulated data.SBC should, in general, be applied to each model to confirm that it is working as expected.Figure 3 shows the results of the simulation-based calibration procedure applied to a model of a planet parameterized by Thiele-Innes elements and position angle at the average epoch.For the most part, these histograms do not reveal any systematic biases in the Octofitter sampling procedure.One exception is the "bump" in the position angle histogram.This shape indicates that Octofitter is under-confident and the sampled posterior distribution is wider than the true posterior.By contrast, seemingly small tweaks to the choice of priors can result in histograms with strong biases.For example, drawing the Thiele-Innes constants (A,B) and (F,G) from log-uniform prior distributions rather than Gaussian distributions centred around a single scale, itself drawn from a log-uniform prior, results in noticeable issues with the SBC histograms.A guide to interpreting the results of the SBC procedure is available in Talts et al. (2020).
These tests illustrate the value of performing the simulation-based calibration procedure for each new model and data combination a user wishes to use.Given the complexity of orbit models and the difficulty of sampling from them, we do not expect our sampling procedure to be entirely unbiased.Rather, we hope that by creating an easy way to diagnose these biases, users of Octofitter can interpret their results with an appropri-ate level of skepticism in accordance with the level of bias detected.

Other Packages
Fitting observations of planets and binary stars to orbits has been widely addressed in the literature, dating back to Kepler's seminal work.More recently, a variety of software packages have been released following both frequentist and Bayesian approaches.Some of these packages include EXOFAST (Eastman et al. 2019), rvfit (Iglesias-Marzoa et al. 2015), radvel (Fulton et al. 2018) 4 compared between three packages: orvara, orbitize!, and Octofitter.Semi-major axis is plotted on a log scale to reveal how the sampler behaviour differs in the long tail.Though we expect all packages to eventually converge to near-identical results, we find that there are small differences in the ω, a, and e marginal posteriors that persist even after running for multiple days.
exception of exoplanet, the majority of these codes have used first-order samplers.That is to say, they rely only on evaluating the posterior density and do not calculate or make use of gradient information.

Demonstration with Relative Astrometry
We begin our demonstration of Octofitter by fitting orbits to simulated relative astrometry measurements, which is one of the simplest use cases for the package.We considered a simulated orbit of a single planet and generated astrometry at 27 epochs (Table 1).For the sake of comparison, we performed the same fit using two other popular orbit fitting packages: orvara and orbitize!.
We followed the best practices laid out in the tutorials provided with each package.orvara and orbitize!use slightly different priors by default and cannot be made to match each other exactly without code modifications.Orvara uses an √ e sin(ω) and √ e cos(ω) parameterization while orbitize!uses separate uniform priors on both e and ω.Therefore, we ran our comparisons twice, first adopting priors similar to orvara and then priors similar to orbitize!In all cases, we used a log-uniform prior on semi-major axis between 0.1 and 1000 AU.We ran the orvara and orbitize!packages with settings recommended by their authors.These were 4 temperatures with 100 walkers for orvara and 20 temperatures with 1000 walkers for orbitize!.We ran Octofitter with 16 independent chains.We initialized the walkers used by orvara and orbitize! in a Gaussian ball around the true orbit.For orbitize!, we improved the convergence by randomly initializing half the walkers on the second mode of the ω marginal posterior.For Octofitter, on the other hand, we initialized it automatically using our procedure of drawing from the priors and selecting the parameters with the highest posterior density.
We drew 20,000,000 samples from each walker using orvara and 100,000 from each walker using orbitize!.For Octofitter, we adapted the step size and mass matrix for 5,000 iterations and then drew a further 15,000 samples.The resulting posteriors are compared in Figure 4. To remove the burn-in phase, We discarded the first halves of the orvara and orbitize!chains and the first 5,000 samples of the Octofitter chains (the adaptation phase).
To measure how long each package takes to converge to a steady distribution, we followed a similar procedure to Ferrer-Chávez et al. (2021a).We divided each chain into 50 segments and assumed that in the final segment, the chain is fully converged.We then calculated the R statistic between each segment and the final segment.We used the rank normalized and median folded version of the statistic as implemented in MCMCDiagnosticTools.jl 7 to evaluate how well the samplers converged in the bulk and in the tails of the distribution.We considered it converged once R became less than 1 ± 0.005 for all variables.We evaluated this on all walkers (orvara and orbitize!) and all chains (Octofitter) and took the median.
The most important result of this code comparison is that in the limit of large numbers of samples, all packages produce posterior distributions that largely agree Comparison between packages of the average time until chains converged to a stationary distribution and the rate at which independent posterior samples are generated.The effective sample size (ESS) rate was measured separately using the bulk and tail methods of MCMCDiagnostics.jl.Note that the R and ESS of the slowest variable for a given sampler are used as this is what ultimately limits sampling performance.These results are based on the astrometry presented in Table 1, and are expected to depend strongly on hardware, input data, and choice of priors.
with each other and that include the true orbit.The orbit paths in the plane of the sky are compared in Figure 4 and the parameters are compared in a corner plot in Figure 5.This should serve to further improve confidence in the results of all packages.
The sole exception is that the orbit posterior from Octofitter contains more samples from a higheccentricity, high-semi-major axis branch of the posterior than the two other packages (Figure 5).Nonetheless, running SBC on this model with uncertainties and observation epochs from this dataset (Figure 3) reveals that Octofitter is acceptably calibrated and is not overestimating the spread of the posterior.It is not computationally feasible to perform the same SBC experiment with the other packages, but this likely indicates that the additional samples from the long tail of the posterior are representative of the true posterior and would eventually appear in the outputs of the other packages.
A second result is that for this problem, Octofitter converges to a steady distribution almost immediately using a single computer core, while other packages take many hours to converge (Figure 6).Since assessing the convergence of MCMC chains can be difficult, we suggest that with Octofitter, unsophisticated users are therefore less likely to use insufficiently converged results in their analyses.Of course, in addition to needing the chain to be converged, one must also generate a sufficient number of samples for the problem at hand.Interestingly, despite the much faster convergence, Octofitter is no faster than the other packages at producing independent samples.In most workflows, the bottleneck is ensuring convergence rather than a need to produce tens of thousands of samples, so Octofitter should still offer a decisive improvement in computation time.
The exact results are hardware specific, sensitive to the data, and depend on the choice of priors and parameterization.It is not feasible to fully account for all small differences between software packages, and it is likely that different approaches will perform better or worse depending on the problem at hand.An additional caveat is that the ptemcee sampler used by orvara and orbitize!scales across many cores which reduces the total sampling time.Octofitter supports running multiple chains in parallel but is not configured to use multiple cores to accelerate a single chain.Finally, given the convergence guarantees of Markov chain Monte Carlo, all packages would approach the same distribution given infinite time.This comparison, therefore, is meant only to illustrate the typical efficiency one might expect with Octofitter on similar problems and reasonable run times.

Demonstration with Relative Astrometry, RV, and PMA
We now demonstrate Octofitter on a real system with radial velocity, proper motion anomaly, and relative astrometry data.
The HD 91312 system consists of a 1.6M ⊙ star orbited by a binary companion with a mass of approximately 0.3M ⊙ and separation of roughly 10 AU (Chilcote et al. 2021).The companion was discovered by a direct imaging search targeting accelerating stars (Currie et al. 2021) from the HGCA.We now reproduce the orbital analysis of the discovery paper in order to demonstrate Octofitter's modelling capabilities when applied to a combination of relative astrometry, radial velocity, and proper motion anomaly data.
For this analysis, we used the astrometry data from Chilcote et al. (2021), proper motion anomaly data from the HGCA (Brandt 2021), and limited radial velocity data originally from Borgniet et al. (2019).As the original radial velocity data was not forthcoming, we measured it from the PDF plots submitted to the arXiv alongside the 2021 manuscript.We used similar priors as those described in the discovery paper, namely lognormal priors on primary and companion mass and a uniform prior on the square root of eccentricity.These were chosen to match the analysis of the discovery paper for the purpose of comparing codes and demonstrating Octofitter rather than any physical motivation on our part.
The results of this orbit modelling are presented in Figure 7. Octofitter recovers the orbit of the companion with similar results to those presented by Chilcote et al. (2021).The radial velocity, proper motion anomaly, and relative astrometry are all consistent with the secondary companion having a mass of approximately of 300 M jup , given the choice of priors described in the original discovery paper.

Demonstration with Images
We now present a series of simulations showing how this framework allows us to detect planets using multiepoch direct images.We selected 50 sequences from the Gemini Planet Imager Exoplanet Survey (Nielsen et al. 2020), processed using a forward model matched filter (FMMF, Ruffio et al. 2017), that have a stellar I-band magnitude less than or equal to 6 and have no previously detected point sources.To be clear, we do not search these sequences for real companions but merely use them as realistic noise maps for our simulations.We normalized the contrasts of each sequence to the median at 200 mas separation.This retains the true noise distribution but removes the effects of sequence-to-sequence variation on our results.We consider a hypothetical system at 20 pc with a 1M ⊙ star that is observed once per month for between 1 and 100 months (a little over eight years).Using these sequences and parameters, we generated simulated observations of a planet by injecting a synthetic PSF into the correct positions in each epoch.This input data is presented in Figure 8.
For all models, we adopted a 1 ± 0.1M ⊙ prior on M , a uniform 0 − 30M jup prior on m, a Gaussian prior on parallax, a uniform prior on a, and a uniform circular prior on τ .
We test this model's ability to detect planets in sequences of direct images for the most straightforward case: circular, face-on orbits.We injected the planet into 1, 5, 10, 25, 50, and 100 images with an average SNR ranging between 0 and 5 in each image.Finally, we applied our model to each of these simulated datasets to arrive at a grid of recovered SNR values as a function of number of epochs and SNR per epoch (Figure 9).We find that we are able to recover planet detections with arbitrarily low SNR per epoch despite orbital motion, provided we have a sufficient number of observations.Figure 9 also shows how the recovered SNR compares with the SNR we would expect for combining images without orbital motion in the presence of uncorrelated Gaussian noise.We find that the model detects the injected planets with near perfect √ N scaling when the final, combined SNR is greater than ∼5.For instance, a planet injected into 100 epochs spaced one month apart, at an SNR per epoch of just 0.5 (below the noise), is still robustly detected at a final SNR of 5. We do note that since each image was taken of a different star, this ideal scaling is a best-case scenario.It is possible that repeated observations of the same target could lead to correlated residuals that reduce the final SNR, though the orbital motion of the planet should mitigate this in much the same way as angular differential imaging (Males et al. 2015;Marois et al. 2006).
The left-hand columns of Figure 9 are of particular note.They show that the model fails to detect a point source injected into a single image at an SNR of 4. The failure to detect planets with expected significance below SNR ≈ 5 can be understood by contrasting our detection criteria with the standard used in the field.In these models, we consider the overall SNR marginalized over all locations in the image (or orbits through a sequence), whereas typically, one looks at the maximum  √ N improvement with the number of epochs.We find that the SNR grows as expected unless the final, combined SNR is below ≈ 5.The recovered SNR falls off quickly below this value and levels off at approximately 1.This explains why the recovered SNR of 1 exceeds the injected value of 0.5 in the bottom left corner.
SNR at any given location.For example, in any SNR 5 detection of a planet, numerous other SNR 2 and 3 peaks exist.When looking at the SNR marginalized over all locations, these other peaks serve to reduce the final SNR.This makes the SNR calculated from the marginal photometry posterior a more stringent planet detection test.

Demonstration with Aperture Masking Interferometry and Radial Velocity
We now demonstrate Octofitter using a combination of simulated radial velocity data and aperture masking interferometry (AMI) visibilities.We considered a plausible scenario where a planet has been detected using radial velocity, had its orbit characterized, and is then followed up with a series of three observations with JWST NIRISS AMI (Doyon et al. 2012;Sivaramakrishnan et al. 2023).For this experiment, we connected the mass of the planet to its M s band photometry using Sonora Bobcat models (Marley et al. 2021) for a fixed system age of 10 Myr.The simulated system had a true orbit with a semi-major axis of 2 AU, an eccentricity of 0.1, an inclination of 45 • , a parallax distance of 100 mas, and a mass of 5 M jup .The star had a stellar absolute M s-band magnitude of 2.6.
We simulated the radial velocity by adding Gaussian noise with a standard deviation of 25 m/s to 74 epochs of RV data generated from the above system.The data points were spaced by ∼30 days between 2017 and 2023.
The AMI data was generated using Eq. 3 for a single companion following the given orbit.For this example, we used only the closure phase and did not include the squared visibilities.We again added Gaussian noise to the calculated closure phases, with a closing triangle dependent standard deviation taken from the closure phase uncertainties calculated with AMICAL (Soulain et al. 2020) from NIRISS AMI F480M data of (the presumed single stars) HD 101531 and HD 123991 calibrated against each other.
We applied Octofitter with a joint model of the RV data and three AMI epochs.We also compared this to χ 2 maps from the Fouriever framework (presented in Kammerer et al. 2023) of each individual epoch.The results of this experiment are shown in Figure 10.
Using the AMI data, the Octofitter model is able to constrain the inclination and, therefore, the mass of the planet.Due to the faint signal and presence of noise, the Fouriever results show spurious signals at each epoch.In comparison, Octofitter is able to connect the three epochs by a single higher-significance mode in orbitparameter space.
This example illustrates how joint modelling across epochs can be used to increase the significance of AMI and other similar interferometric observations.

CONCLUSION
We have presented a new code, Octofitter, for modelling exoplanet relative astrometry, radial velocity, and proper motion anomaly, as well as performing nontraditional tasks like detecting moving point sources across images and interferometric observables.
• We demonstrated the simulation-based calibration procedure on a hypothetical orbit fitting task and found that for orbits parameterized with Thiele-Innes constants, Octofitter is acceptably calibrated.
• We compared the results of Octofitter to the popular packages orvara and orbitize!and found that all three arrive at similar posterior distributions.• We showed that Octofitter can converge 1 to 2 orders of magnitude faster than other popular packages, making it computationally feasible to perform a variety of statistical checks like simulationbased calibration for individual datasets.
• We demonstrated a combined fit to relative astrometry, radial velocity, and proper motion anomaly of the HD 91312 system and found a companion mass that agreed with previous results.
• We demonstrated the ability to detect arbitrarily faint companions despite orbital motion using simulations and data from the Gemini Planet Imager Exoplanet Survey.
• We demonstrated a combined model of simulated radial velocity and multi-epoch JWST/NIRISS aperture masking interferometry.
Octofitter is a powerful new tool that will enable the community to broadly apply joint multi-epoch, multiinstrument, and multi-band direct imaging, interferometry, Doppler radial velocity, and astrometric motion.This will allow for the detection and characterization of fainter and lower-mass planets more similar to those in our solar system.
Facilities: Gemini (GPI), Gaia, Hipparcos To properly fit to relative astrometry, proper motion anomaly, and radial velocity data, expressions for the position, velocity, and acceleration of a secondary companion (such as a planet) about a primary mass (such as a star) in Cartesian and celestial coordinates are required.It is well-known that the solution to the Keplerian two-body problem A.2. Celestial Coordinates Now that we know r obs (t), v obs (t), and a obs (t), we can calculate the corresponding expressions in celestial coordinates.Before doing so, a brief discussion of coordinate conventions is required.For this paper, we measure with reference to the celestial north pole.This means that, in reference to the previous section, the positive x-axis points in the direction of positive declination (upwards on the sky), while the positive y-axis points in the direction of positive right ascension (leftwards on the sky).This choice of coordinates is important to note, since it differs from the right/up orientation of the x/y axes commonly seen elsewhere.Furthermore, the positive z-axis points away from the observer.
Keeping this coordinate convention in mind, we can determine position, velocity, and acceleration in celestial coordinates.Normally, an observer a distance d away from an object of size r measures the object to have an angular size of ∆θ = arctan r d . (A20) However, for scenarios relevant to this work, r is on the order of astronomical units and d is on the order of parsecs, meaning the small angle approximation ∆θ ≈ r/d can be used with negligible error.Applying this approximation, we get that the right ascension and declination offsets of the secondary from its primary are where d is the distance to the system, and x obs and y obs are the x-and y-components of Eq.A10.From this, it immediately follows that where v x,obs and v y,obs are the x-and y-components of Eq.A14, and that where a x,obs and a y,obs are the x-and y-components of Eq.A16.The equations for ∆RA and ∆DEC can be used to fit orbital and physical parameters to relative astrometry data, while the equations for ∆ ṘA and ∆ ḊEC can be used to fit orbital and physical parameters to proper motion anomaly data.∆ RA and ∆ D EC are used for calculating model derivatives for higher-order samplers, and may be applied to angular acceleration data in the future.

B. DERIVATION OF THE PARAMETER τ FROM POSITION ANGLE θ
The parameter τ used in orbitize!and Octofitter is the fraction of the orbit completed by a planet after periastron, at some reference epoch t ref .This is a useful parameterization since it lies between 0 and 1 for all orbits; however, if the reference epoch is not similar to the epochs of the observations, it can exhibit pathological behaviour as the other orbital parameters change.One solution is to adjust the reference epoch for each dataset, but a further improvement for many cases is to instead adopt the observed position angle θ at epoch t as an independent variable.θ is an improvement over τ because it is insensitive to changes in the other parameters and has a straightforward interpretation.A derivation of τ from θ is given here.
We begin with the parallax distance to the system ω, the host mass M , the eccentricity e, and the Thiele-Innes constants A, B, F, and G given in Eqs.C29 -C32.From these, we write the transformation matrix  From this, the parameter τ can be calculated using the period as follows.First, calculate the semi-major axis a using Eq.C35.Next, combining Kepler's third law and the value of a gives us the period P = 2π a 3 /(GM ).The mean motion is then n = 2π/P , and the epoch of periastron passage follows from this as t peri = t − MA/n − t ref .Finally, )

Figure 1 .
Figure 1.Conceptual schematic of three different ways to use Octofitter.Other possibilities, like combining images with Doppler or astrometric velocimetry or using a mix of relative astrometry and images without detections to constrain orbits, are also possible..
the Gaia (GAIA-Collaboration et al. 2021) and Hipparcos (van Leeuwen 2007) missions.The Hipparcos-Gaia Catalog of Accelerations (HGCA Brandt, T. D. 2021) cross-calibrates measurements from these two satellites and inflates uncertainties to cover most known systematics.This results in projected velocity measurements of the system's photocentre at the Hipparcos epoch, the Gaia epoch, and between the two.
Figure 2.Posterior density of photometry versus semimajor axis for a simulated detection and non-detection.The solid lines mark the mean of the marginal photometry posterior, and the dashed lines mark ±1σ.We consider a planet detected when the SNR based on the photometry marginalized over all other parameters is greater than some chosen threshold, e.g.5σ.
t a b l e = T a b l e ( C S V .r e a d ( " a s t r o m .c s v " ) ) a s t r o m = A s t r o m e t r y L i k e l i h o o d ( t a b l e ) @ p l a n e t b V i s u a l { K e p O r b i t } b e g i n a ~L o g U n i f o r m ( 2 .5 , 2 5 ) i ~S i n e ( ) e ~t r u n c a t e d ( N o r m a l ( 0 .2, 0 .2 ) , l o w e r = 0 , u p p e r = 1 ) Ω ~U n i f o r m C i r c u l a r ( ) ω ~U n i f o r m C i r c u l a r ( ) τ ~U n i f o r m C i r c u l a r ( 1 .0 ) e n d a s t r o m¦ ¥A similar syntax is used to specify stellar properties: Figure 3.Sample results of running the simulation-based calibration procedure on a model consisting of a single planet parameterized with Thiele-Innes elements (A, B, F , G, etc.) and position angle θ.Each count represents a model fit to a different simulated system.The horizontal band gives a 99% range around the expected value for a perfect sampling procedure.The observation epochs and uncertainties are taken from Table1.

Figure 4 .
Figure 4. Orbit fitting posteriors visualized in the plane of the sky, compared between three packages: orvara, orbitize!, and Octofitter.

Figure 5 .
Figure5.Corner plot of key orbital parameters for the case shown in Figure4compared between three packages: orvara, orbitize!, and Octofitter.Semi-major axis is plotted on a log scale to reveal how the sampler behaviour differs in the long tail.Though we expect all packages to eventually converge to near-identical results, we find that there are small differences in the ω, a, and e marginal posteriors that persist even after running for multiple days.

7
Figure 6.Comparison between packages of the average time until chains converged to a stationary distribution and the rate at which independent posterior samples are generated.The effective sample size (ESS) rate was measured separately using the bulk and tail methods of MCMCDiagnostics.jl.Note that the R and ESS of the slowest variable for a given sampler are used as this is what ultimately limits sampling performance.These results are based on the astrometry presented in Table1, and are expected to depend strongly on hardware, input data, and choice of priors.

Figure 7 .
Figure 7. Sample plot output from Octofitter using data from the HD 91312 system.The top row visualizes the orbital posterior compared with velocity measurements of the host.The horizontal bars in the proper motion panels show the timespans over which the average velocity was measured.The middle row shows the posterior compared with astrometry measurements in the plane of the sky and in separation and position angle over time.The bottom row shows the orbital posterior in physical coordinates to compliment the astrometry plot.The rightmost panel shows a deprojected view of the system where orbits have been rotated face on and to place periastron at the bottom.The conventions used by Octofitter are described in Appendix A.

Figure 8 .
Figure 8. GPI sequences used for simulations in this section.Each image was normalized to have the same average contrast at 200 mas separation from the star (just outside the edge of the mask).The images displayed above contain a simulated planet orbiting CCW at an average of just SNR 1 per epoch, spaced one month apart.Given this data, the model recovers the simulated planet at SNR 7.

Figure 9 .
Figure9.Octofitter's ability to recover planets from simulations of circular, face-on orbits as a function of signal to ratio per epoch and number of epochs.Left: Recovered SNR.Cells above the yellow line would be detected with a 5σ threshold.The cell outlined in cyan corresponds to the data shown in Figure2.Right: The same SNRs relative to an ideal √ N improvement with the number of epochs.We find that the SNR grows as expected unless the final, combined SNR is below ≈ 5.The recovered SNR falls off quickly below this value and levels off at approximately 1.This explains why the recovered SNR of 1 exceeds the injected value of 0.5 in the bottom left corner.

Figure 10 .
Figure 10.Joint modelling of AMI and RV data.A-C: Recovered χ 2 maps at each independent AMI epoch created using Fouriever (presented in Kammerer et al. 2023).D and E: joint radial velocity and AMI posteriors visualized as an RV time series (D) and orbit posterior in the plane of the sky (E).

Figure 11 .
Figure 11.Comparison of θ (blue, second last row) and τ (yellow, last row) paramterizations of the same orbital posterior.The θ marginal posterior is nearly Gaussian, and has simple relationships with all other orbital parameters.By contrast, the τ marginal posterior is less informative and has complex relationships to ω and Ω. Complex structures in a posterior decrease sampling efficiency and could in some cases lead to biased results.
σ_pa, σ_sep, and cor to specify the measurement uncertainties and correlation, respectively.Alternatively, ra (relative), dec (relative), σ_ra, and σ_dec can be substituted if preferred.A sample relative astrometry input table is presented in Table

Table 2 .
Images input sample

Table 4 .
Radial velocity input sample A. and Dec. proper motion at each epoch and σ RA and σ ṽDEC are the associated uncertainties.The input table format for this observation type matches Table 4 of Brandt, T. D. (2021).Compared with Brandt, T. D. et al. (2021) and Brandt, G. M. et al. (