ethraid: A Simple Method for Characterizing Long-period Companions Using Doppler, Astrometric, and Imaging Constraints

We present ethraid, an open-source Python package designed to measure the mass (m c ) and separation (a) of a bound companion from measurements covering a fraction of the orbital period. ethraid constrains m c and a by jointly modeling radial velocity, astrometric, and/or direct imaging data in a Bayesian framework. Partial orbit data sets, especially those with highly limited phase coverage, are represented well by a few method-specific summary statistics. By modeling these statistics rather than the original data, ethraid optimizes computational efficiency with minimal reduction in accuracy. ethraid uses importance sampling to efficiently explore the often broad posteriors that arise from partial orbits. The core computations of ethraid are implemented in Cython for speed. We validate ethraid's performance by using it to constrain the masses and separations of the planetary companions to HD 117207 and TOI-1694. We designed ethraid to be both fast and simple, as well as to give broad, “quick look” constraints on companion parameters using minimal data. ethraid is pip installable and available on Zenodo and GitHub.


INTRODUCTION
Fitting the orbits of exoplanets is one of the best tools for uncovering their origins, evolutionary history, and relationship to neighboring planets.Keplerian fits to RV observations produced the first detections of hot Jupiters orbiting Sun-like stars (Mayor & Queloz 1995;Butler et al. 1997;Cochran et al. 1997), and higherprecision measurements have facilitated mass determinations of smaller planets (Pinamonti et al. 2018;Luque et al. 2019).
Planets with long periods, from years to decades, present multiple observational challenges.In particular, their orbits require consistent and prolonged monitoring to characterize.Long-term RV surveys (see, e.g., Rosenthal et al. 2021) confront this problem directly by compiling observational baselines over multiple decades.Though fruitful, these efforts are costly in terms of both human time and telescope time.
Furthermore, these surveys are observationally inefficient in two ways.First, the precise mass and orbital parameters of each companion are often not of immediate interest.For example, studies of giant planet occurrence rates (Fulton et al. 2021) or the correlation between gia) ethraid (2024) b) https://github.com/jvanzand/ethraidants and other planet classes (Bryan et al. 2019;Rosenthal et al. 2022) are principally concerned with differentiating between giant planets, brown dwarfs, and stars, which may be accomplished with less than a full orbit.Although precise constraints may be useful, some science cases benefit from greater statistical breadth over precision measurements of specific orbital parameters.
Second, once legacy surveys end, their final catalogues may still include targets with too little phase coverage to characterize (Rosenthal et al. 2021).Without further observation, these time series represent significant investments of telescope time with minimal scientific return.Both of these limitations motivate tools for extracting planetary information from partial orbits.
Many existing Bayesian orbit fitters are well-suited to performing precise fits using RVs (radvel, Fulton et al. 2018; The Joker, Price-Whelan et al. 2017), astrometry (orbitize!, Blunt et al. 2020;OFTI, Blunt et al. 2017), or both (orvara, Brandt et al. 2021).Some of these codes are designed to be robust against specific orbit fitting challenges (e.g., non-Gaussian posteriors, low phase coverage, and high-dimensional parameter spaces).However, their performance tends to suffer as these limitations are taken to the extreme.
In this work, we present ethraid, an open source Python package designed specifically to constrain companion mass and separation given partial phase cov-erage.ethraid can simultaneously model summary statistics of up to three independent data types: linear/quadratic RV trends, astrometric trends from the Hipparcos-Gaia Catalog of Accelerations (Brandt 2018(Brandt , 2021)), and direct imaging contrast curves.In the remainder of this paper, we refer to these summary statistics as the "data." This paper takes the following structure.In Section 2 we describe ethraid's fitting algorithm, including parameter sampling, forward modeling and likelihood calculation, and marginalization to derive posteriors.Section 3 gives the mathematical details of our forward model for each of RVs, astrometry, and direct imaging data.We discuss ethraid's performance in Section 4 and basic usage in Section 5. We test the strengths and weaknesses of this code in Section 6.Finally, we give a brief list of future improvements planned for ethraid in Section 7 before concluding in Section 8.

ETHRAID'S FITTING ALGORITHM
Given some RV, astrometric, or imaging data, ethraid samples a space of orbital models, assesses each model's probability of having produced the signal based on both its a priori probability and the likelihood of the measured signal in light of the model, and uses the resulting posterior surface to calculate confidence intervals for the inferred companion's mass and semi-major axis.ethraid's approach to this problem consists of three steps: 1. Random sampling of orbital parameters from prior distributions 2. Forward modeling and likelihood calculation

Marginalization over orbital parameters
These steps are summarized in Figure 1.Before detailing them below, we note an important assumption that ethraid makes, namely, that any measured signals were produced by a single bound companion.Consequently, RV and astrometric trends produced by multiple companions, stellar activity, or instrumental systematics will not be properly interpreted and may produce misleading results (see Section 6).

Parameter Sampling from Priors
The algorithm used to sample observations (in our case, model parameters) from a distribution is a key component of any Monte Carlo method.ethraid uses importance sampling (Kloek et al. 1978), an approach in which observations θ are sampled according to an importance function q(θ).The result is that the histogram of a set of N samples, {θ n } N n=1 , converges to q in the limit of large N : where "∼" denotes "has the distribution of."However, because our goal is to calculate the posterior probability p(θ), we must weight each observation: where the weights w n are given by Following the earlier partial orbit fitters OFTI and The Joker, we sample orbital elements directly from their respective priors, π(θ).This choice is motivated by the fact that partial orbit data is generally uninformative, resulting in prior-dominated posteriors.Conveniently, choosing q = π(θ) simplifies the weight of each model to the model likelihood: Thus, ethraid approximates the posterior surface p(θ) by creating a histogram of orbital models, each sampled from the appropriate prior distributions and weighted by its likelihood.The first step of this procedure is to sample six orbital parameters according to their respective priors: semi-major axis a, eccentricity e, inclination i, argument of periastron of the companion (not the host star, as is sometimes the case in RV codes; see Householder & Weiss 2022) ω, the mean anomaly at a reference epoch M 0 , and companion mass m c .As all data products we model are insensitive to the longitude of the ascending node (Ω), we fix it to 0 arbitrarily.The prior PDFs for these parameters are: where U denotes a uniform distribution and E is a userselected eccentricity distribution.Options for the latter include -'zero': a constant distribution of e = 0 -'uniform': a uniform distribution between 0 and 0.99 -'kipping': the piecewise beta distribution for short-and long-period exoplanets presented by Kipping (2013) -'piecewise': a combination of the Kipping (2013) distribution for planets ( mc MJ ≤ 13), the beta distribution derived by Bowler et al. (2020) for brown dwarfs (13 < mc MJ ≤ 80), and a uniform distribution between 0.1 and 0.8 for stars (80 < mc MJ ) found by Raghavan et al. (2010).
While the last of these is the most physically motivated in general, the others give users the option to include prior knowledge about a given system.For example, a dynamically cool system with multiple known planets on circular orbits may be unlikely to host an eccentric outer giant.In this case, a user might opt for the 'kipping' distribution to favor lower-eccentricity models.
Our importance sampling strategy differs from MCMC sampling (e.g.radvel), where each model draw in a chain is dependent on the previous step.Drawing samples independently is advantageous for characterizing highly non-gaussian or multi-modal posterior surfaces, which arise often in partial orbit analysis.Our method is similar to the rejection sampling technique used by OFTI and The Joker, which also draws models independently from a proposal distribution.However, rather than accepting or rejecting models based on their posterior probability, ethraid assigns each model a weight according to its likelihood.In practice, models with a high chance of rejection under rejection sampling are assigned vanishingly small weights under importance sampling, yielding similar performance.

Forward Modeling and Likelihood Calculation
After sampling a set of parameters, ethraid evaluates the model likelihood, that is, the probability of observing the measured data given the model under consideration.ethraid accomplishes this by generating synthetic RV, astrometric, and/or imaging data sets, and comparing them to the observed data.This procedure is repeated for a number of model orbits specified by the user (typically 10 6 -10 8 ).The likelihood calculations for each data type are detailed in Section 3.
In the case that multiple data sets are provided, ethraid records the likelihood of each data set conditioned on each model.Because the RVs, astrometry, and imaging are independent, ethraid calculates the combined likelihood of all provided data sets as the product of the likelihoods of each individual data set.

Marginalization
ethraid can generate visualizations of the 1D and 2D marginalized posteriors for a and m c .The marginalization step displays the desired posterior by summing the likelihoods of all models in the same interval of parameter space.For a 1D posterior, these are intervals in the desired parameter, for example, ∆m c,i .For a 2D posterior, they are interval pairs (∆a i , ∆m c,i ).

FORWARD MODEL
In this section we detail ethraid's likelihood calculations for RV, astrometry, and imaging data, including the forward modeling procedure which produces an analytic counterpart to each measured quantity.

RV Constraints
If a body on a long-period orbit is left unmodeled, it will present as a gradual increase or decrease in the RV residuals.For an orbital period much longer than the observing baseline, the variation is almost purely linear; for a period only a few times the baseline, there may also be significant quadratic curvature.The RV data therefore comprises the linear and quadratic coefficients to a second order polynomial fit to the RVs, denoted γ (m/s/day) and γ (m/s/day 2 ), respectively.
For each set of sampled parameters, our goal is to forward model γ and γ, whose analytic forms are given by differentiating the stellar radial velocity γ: where K is the RV semi-amplitude: and ν is the true anomaly, related to the eccentric anomaly E by (Murray & Dermott 2010) The derivatives of ν are We obtain E by numerically solving Kepler's equation: where M is the mean anomaly, calculated by advancing the sampled initial mean anomaly M 0 some duration t past a reference epoch via where P is the orbital period calculated from Kepler's third law.We choose our reference epoch to be 1989.5 (see Section 3.2).We evaluate the model log-likelihood by calculating γ and γ at the epoch of the RV observations and comparing them to their measured counterparts via where γdata and γdata are the measured slope and curvature of the RV time series, σ γdata and σ γdata are their respective uncertainties, and C RV is a constant.From Equations 5 one may derive (see Appendix B) that for small companion masses, surfaces of constant γ and γ follow m c ∝ a 2 and m c ∝ a 7/2 , respectively.Thus, for any sampled companion mass, there is a physical separation at which such a companion would produce the observed RV trend and curvature.In other words, RV constraints are generally consistent with a large range of (a, m c ) pairs.We seek to break this degeneracy by incorporating an independent data set: astrometry.

Astrometry Constraints
A star with no orbiting companions will have a constant proper motion.Meanwhile, a massive companion will induce a change in its host star's proper motion vector over some time interval, called the proper motion anomaly (PM a ; Kervella et al. 2019).Our goal is to model the magnitude of the PM a vector.
Our forward modeling procedure for astrometry is similar to the RV case, but is tailored to the data provided by the Hipparcos-Gaia Catalog of Accelerations (HGCA; Brandt 2021).This catalog aligned the reference frames of Hipparcos (1989.85-1993.21, Hip 1997) and Gaia EDR3 (25 July 2014-28 May 2017; Lindegren et al. 2021) to measure the accelerations of over 115,000 stars.We choose 1989.5 as the reference epoch t 0 for our RV, astrometric, and direct imaging models.
The HGCA reports three proper motions for each target: the epoch proper motions of Hipparcos (⃗ µ H ) and Gaia (⃗ µ G ) near 1991.25 and 2016.0,respectively, and the average proper motion given by the difference in position between these epochs divided by their ∼25-year baseline, ⃗ µ HG .The position-derived proper motion is generally the most precise proper motion measurement owing to the long baseline between the two missions, while ⃗ µ G is usually the next-most precise.We therefore quantify PM a as ∆⃗ µ, the difference between ⃗ µ G and ⃗ µ HG , and model its norm: We likewise obtain σ ∆µ by propagating the measurement uncertainties provided in the HGCA.The astrometric picture is summarized in Figure 2.
Our derivations of ⃗ µ G and ⃗ µ HG differ in one key respect from the procedure used to construct the HGCA.Brandt (2021) obtained the best-fit positions and proper motions by recalibrating astrometric fits to the full Hipparcos or Gaia time series.Because simulating and fitting a model to a full time series would be computationally expensive, we instead calculate the analytic average of position and proper motion during each mission by integrating over the host star's orbit.We take these average quantities to approximate the stellar position and proper motion at the mission midpoints.To minimize positional uncertainty, Brandt evaluated the fitted positions at target-specific "characteristic epochs" rather than the mission midpoints.This disparity introduces an error on the order of 0.01 mas/yr to our model, which we consider negligible because it is below both the median (0.047 mas/yr) and minimum (0.016 mas/yr) ∆µ uncertainty in the HGCA.
Below we derive expressions for the stellar position and velocity due to a companion, and use them to model ∆µ.Calculating ⃗ µ HG requires the average stellar position over each mission, while ⃗ µ G requires the average stellar proper motion over only the Gaia mission.
In a two-body system composed of a host star and a bound companion, the star's position in the orbital plane relative to the system barycenter is (Murray & Dermott 2010, eq. 2.41) where a ⋆ = a mc m⋆+mc is the semi-major axis of the star's orbit and the positive X-axis is in the direction of the companion's periastron.We use uppercase letters to refer to coordinates in the orbital plane, and lowercase letters for the sky plane.The average position over a time interval [t 1 , t 2 ] can be calculated analytically in the usual way: for the X-component, it is with a similar formula for ⟨Y ⟩.We may reformulate this integral as: We have X(E) from Equations 13 and dt dE = P (1−e cos E) 2π can be obtained from Equation 9. Carrying out the above integral and the corresponding one for Y yields To calculate average proper motion over the same time interval [t 1 , t 2 ], we simply evaluate Equations 13 at the interval boundaries and divide their difference by the elapsed time: Equations 16 and 17 allow us to define the average stellar position and proper motion relative to the system barycenter in the orbital plane: To express these vectors in the observer frame, that is, with the x-and y-axes in the sky plane and the zaxis toward the observer, we apply the rotation matrix R (Murray & Dermott 2010), which is given in its full form in Appendix C.

⃗ r obs
We project these 3-dimensional vectors onto the sky by ignoring the z-component, and convert to angular coordinates using the system's distance from Earth.The anomalous position ⃗ ρ anom and proper motion ⃗ µ anom over either Hipparcos or Gaia are given by and the absolute position and proper motion at the mission midpoint t mid are where the subscript '0' refers to the barycentric position and proper motion at the reference epoch t 0 and α * ≡ α cos δ, where δ is the Dec at the epoch at which α is evaluated.1 We use the above results to estimate 1) {α H,abs * , δ H,abs }, the stellar RA/Dec at the Hipparcos midpoint t H,mid ; 2) {α G,abs * , δ G,abs }, the stellar RA/Dec at the Gaia midpoint t G,mid ; and 3) {µ α,G,abs * , µ δ,G,abs }, the stellar proper motions in the RA/Dec directions at t G,mid .We provide here the procedure for calculating the RA component of the position-derived proper motion µ HG,abs .The calculation for the Dec component is analogous.
Combining ⃗ µ HG,abs with ⃗ µ G,abs from Equations 21, we have our final result: The modeled ∆µ is thus a function only of the anomalous proper motions, i.e., those induced by the companion.As in Section 3.1, we calculate the model loglikelihood as RV and astrometric trends are often insufficient to constrain a companion's properties, either because one data set is not available (e.g., any star that is not in the original Hipparcos catalog does not have an HGCA acceleration) or because the constraints they do provide rule out the same regions of parameter space, leaving the same a-m c degeneracy described in Section 3.2.In these cases, we turn to direct imaging to place an upper limit on the companion mass.
Our data is the measured contrast curve, a table of angular separations and corresponding magnitude differences indicating the dimmest detectable companion at a given projected separation from the host star.For simplicity, we treat the contrast curve as a step function with no associated uncertainties.That is, at a given angular separation, any companion dimmer than the listed limit is completely undetectable, whereas any companion brighter than the limit would be detected in all cases.Our goal is to calculate a model companion's angular separation and magnitude difference to compare to the contrast curve.We begin with the angular separation.

Angular Separation Calculation
Similar to Equations 13, we calculate the hostcompanion separation via where we now use the full semi-major axis a instead of a ⋆ .Rotating the separation vector into the sky plane with the rotation matrix R and dividing by the distance to the system d, we obtain the projected angular separation ρ: This approach is fully consistent with Sections 3.1 and 3.2 in that it incorporates all orbital parameters into the modeled quantity.
The angular separation calculation above accounts for the fact that each model companion's projected separation is a function not only of a, but also of e, i, ω, and E.
However, this detailed procedure is expensive, requiring nearly 300% and 70% the time of the RV and astrometry calculations, respectively.We detail an approximate calculation of ρ in Appendix A which significantly decreases run time at the cost of a slight reduction in accuracy.

Contrast Calculation
After modeling a companion's angular separation, we seek to convert its mass m c to a ∆mag contrast in the appropriate band.We model contrast as dependent on m c and the stellar V mag only.We linearly interpolated magnitude-mass relations for both stellar (Pecaut & Mamajek 2013, Table 5) and brown dwarf (Baraffe et al. 2003, Table 4) companions in multiple bandpasses to define a function F (m c ) which converts a model companion mass to ∆mag in whichever of the following bands most closely matches the band of the given contrast curve: {V , R, I, J, H, K, L ′ }.Pecaut & Mamajek (2013) do not include magnitude relations for the K or L ′ bands, though they do for K s and W 1. We approximated that K s ≈ K (2.15 µm ≈ 2.2 µm) and W 1 ≈ L ′ (3.37 µm ≈ 3.78 µm) to concatenate the two tables.The brown dwarf mass-magnitude relations of Baraffe et al. (2003) correspond to mature systems (5 Gyr).
We also performed a linear interpolation of the measured contrast curve, yielding another function f data (ρ) which converts measured angular separation to measured ∆mag contrast.The imaging 'data' is then the function f data , together with the assumption that at any separation ρ ′ , no companion with ∆mag < f data (ρ ′ ) was detected.
We therefore have that for a given model separation ρ and contrast ∆mag = F (m c ), both derived from the initial model parameters θ, the log-likelihood that no companions were detected with ∆mag < f data (ρ) is given by We visualize ethraid's imaging likelihood calculation in Figure 3.The independence of the RV, astrometric, and imaging measurements allows us to find the log-likelihood of all three data sets in light of a model orbit by simply summing Equations 11, 24, and 27.

PERFORMANCE
Four processes contribute meaningfully to the computational cost of an orbital fit.These include sampling orbital parameters from their prior PDFs and, for each sampled orbit, calculating the likelihood of each of the RV, astrometry, and/or direct imaging data sets.Therefore, the total time required to run an orbital fit is predetermined by the number of orbital models sampled (a value input by the user) and the number of data sets provided.
The time required to sample one model orbit and calculate the likelihood of all three data sets conditioned on that model is ∼5.75 µs using a single core of a 2.6 GHz Intel Core i7 processor.ethraid takes another 0.25 µs to multiply the three likelihoods to produce the combined likelihood, for a total of 6 µs per model.The average fraction of total run time required for each of these steps is: sampling -6.8%, RVs -11.0%, astrometry -46.3%, imaging -31.8%, combined likelihood -4.1%.If the approximate angular separation calculation is used for the imaging posterior instead of the detailed calculation (see Section 3.3), the share of time required by the imaging calculations drops to < 1%, reducing the overall run time by ∼30%.
ethraid implements the core calculations of Section 3 in Cython, which provides a ∼5x speedup over the standard Python equivalent functions, allowing users to produce informative posterior plots with 10 6 − 10 7 sampled models on one core in 6 − 60 seconds.
Two potential improvements to the current approach would substantially increase ethraid's performance.First, the likelihood calculations are trivially parallelizable.Calculating the likelihoods for all data sets concurrently could decrease the total run time by 50%, and further parallelizing within data sets would give greater speedups.Second, the likelihood and sampled parameter lists tend to overload the working memory when the number of model orbits is large (≳ 10 8 ), causing significant slowdowns.An alternative would be to store these arrays in temporary cached files, or to perform all necessary calculations associated with a given model or likelihood and then overwrite it immediately.

USAGE, FITTING, SAVING, AND PLOTTING
ethraid is installable through pip, and basic usage instructions are available on GitHub.
The typical use case of ethraid is for a system with an unexplained RV trend, though it may be run on systems with any combination of constraints from RVs, astrometry, or imaging.
The majority of parameters required to run an ethraid fit are stored in the corresponding system's .pyconfiguration file.An ethraid configuration file contains all of the RV, astrometry, and direct imaging data needed to carry out the forward modeling algorithm laid out in Sections 2 and 3 (with the exception of the optional contrast curve for direct imaging constraints), as Beginning with a set of sampled parameters, ethraid models angular separation (either exactly or approximately) using the orbital elements, and contrast using the companion mass, together with a linear interpolation function derived from Table 5 of Pecaut & Mamajek (2013) and Table 4 of Baraffe et al. (2003), denoted by 'PM13' and 'B03,' respectively.ethraid then determines whether the model companion's brightness would have exceeded the minimum detectability threshold at the model separation in the real imaging data (gold line) and assigns a likelihood of 0 for detectable companions and 1 for non-detectable companions.The example model in the diagram falls above the threshold, so we assign it a likelihood of 0.
well as parameters for saving and plotting the calculated posteriors.We provide descriptions of these inputs in Table 1. 2 5.1.Command Line Interface ethraid supports three different commands from the command line: run, plot, and lims.run samples the six-dimensional parameter space and calculates a posterior for each data set provided, plus an additional combined posterior conditioned on all provided data sets.Using the save parameter in the configuration file, the user has the option to save the 'raw' unbinned posteriors, which they may reshape as desired, or the 'processed' posteriors, which are smaller and load more quickly but cannot be reshaped.All arrays are saved as hdf5 3 files in the specified output directory.
After running a fit, the user may access and plot the saved arrays using the plot command.plot loads a 2 ethraid's configuration files are analogous to those used by orvara, and we model this  The lims command simply prints the 2σ confidence intervals (i.e., 2.5th and 97.5th percentiles) of both the m c and a posterior PDFs to the terminal.It is intended to give a quick quantitative summary of the orbit fit.Users may load the raw or processed posterior arrays directly if they wish to calculate more detailed statistics.We caution that the bare numerical results returned by lims will not necessarily give an accurate summary of the posterior distribution, particularly in cases where the posterior is under-sampled.We encourage users to use both the one-and two-dimensional posterior plots to guide their interpretation of the confidence intervals returned by this command.Finally, the all command runs run, plot, and lims sequentially.We give more information on the command line interface in Table 2.

VALIDATION
We tested ethraid on four systems to judge its strengths and weaknesses in estimating companion parameters.We chose the first two systems, HD 117207 and TOI-1694, to illustrate ethraid's capacity to characterize single companions with fractional orbital coverage.Both systems host RV-characterized giant planets with well-constrained orbital parameters.We find that ethraid's parameter estimates are consistent with these 'ground truth' values, regardless of the planet's phase during the interval of RV observation.We chose the other two systems, HD 114729 and HD 12661, to demonstrate ethraid's failure modes.Both of these systems host two long-period giant companions, violating one of ethraid's key assumptions that RV/astrometric signatures originate from a single companion.ethraid's failure to accurately characterize the companions in either system illustrates the perils of applying this code blindly.
We designed our validation algorithm to mimic real applications of ethraid.For systems exhibiting periodic RV variability, we identified and isolated a subset of the RVs which we then subdivided into 'slices' of equal duration.We chose the number of slices so that the RV variability was approximately linear/quadratic over each slice, emulating a real data set with limited coverage of the companion's phase.We used radvel to fit a second order polynomial to each slice and recorded the fitted trend and curvature terms.We then applied ethraid to each fitted trend/curvature pair and determined whether the results were consistent with the known companion parameters.For all fits, we used the 'piecewise' eccentricity prior and simulated 10 8 model orbits.We included astrometric trends from the Hipparcos-Gaia Catalog of Accelerations (HGCA; Brandt 2021) in our analysis when available.
6.1.Case Study: HD 117207's 7.5-year super-Jupiter is identified using 1.8-year baselines HD 117207 is a chromospherically quiet (log R ′ HK =-5.06)G8 dwarf hosting a super-Jupiter companion, HD 117207 b, with minimum mass and separation initially reported as 2.06 M J and 3.78 AU (Marcy et al. 2005).Rosenthal et al. (2021), using RVs spanning over 20 years, refined these values to m c sin i = 1.87 ± 0.075 M J and a = 3.744 0.059 0.060 AU.We identified a roughly 7.5-year section of the full time series and divided it into four slices.We used radvel (Fulton et al. 2018) to fit linear and quadratic trend terms to these slices.The results of this fitting procedure are shown in the top five panels of Figure 4.
We supplied the linear and quadratic trends from each slice as inputs to a partial orbit fit using ethraid.We also retrieved this star's astrometric trend of 0.13 ± 0.03 mas/yr from the HGCA and included it in our model.ethraid performed each joint RV/astrometric fit in about six minutes.We provide the posteriors of our partial orbit analysis of each slice in Figure 4 and compare these results to the known parameters of HD 117207 b, shown in each panel as a gold star.In all four cases, the 68% confidence interval derived by ethraid encompasses the true planet parameters.
We also validated our results for this system against orvara.We used the same HGCA astrometry and the raw RVs for each slice to emulate our ethraid fits.We ran 10 temperatures and 20 walkers over 10 5 steps, which took about seven minutes (comparable to our ethraid runs).However, these fits did not converge and gave substantially broader mass and separation estimates than we obtained with ethraid.After increasing to 10 6 steps, we obtained parameter estimates comparable to those from ethraid in about 80 minutes.We conclude that ethraid's data compression approach allows it to match the performance of other robust fitters with minimal loss of information in the regime of highly limited orbital coverage and with only RVs and absolute astrometry.

Case Study: TOI-1694's companion mass is bounded by direct imaging
TOI-1694 is an early K-dwarf with low chromospheric activity (log R ′ HK =-5.0) hosting a Jupiter analog with a minimum mass of 1.05 ± 0.05 M J at a separation of 0.98 ± 0.01 AU (Van Zandt et al. 2023).
We divided the full two-year time series into three slices and fit linear and quadratic trend terms to them using radvel (Fulton et al. 2018).The results of this fitting procedure are shown in the top four panels of Figure 5.
We supplied the linear and quadratic trends from each slice as inputs to a partial orbit fit using ethraid.This star is not listed in the HGCA, and therefore has no available astrometric trend.However, we included direct imaging of this system (Mistry et al. 2023) obtained in the I-band (832 nm) with the 'Alopeke speckle imager (Scott et al. 2021) in our analysis.We provide the posteriors from our partial orbit analysis of each slice in Figure 5.Although in this case the RV trend and imaging cannot disambiguate between planets, brown dwarfs, and low-mass stars, they are sufficient to rule out high-  '1d' and/or '2d'; whether to save 1d or 2d posterior plots No brackets or quotation marks needed; separate multiple arguments by a space mass stars as well as planetary and brown dwarf companions beyond ∼10 AU.Additionally, detection of curvature in some slices strongly favors models with shorter periods.The RV posteriors are consistent with TOI-1694 b's true parameters in all cases.6.3.Case Study: HD 114729's two companions conspire to produce misleading parameter estimates HD 114729 is an inactive (log R ′ HK =-5.02), metalpoor ([Fe/H]=-0.22)solar analog (G0 V) hosting a sub-Jupiter companion with a minimum mass and semimajor axis of m c sin i = 0.892 ± 0.053 M J and a = 2.094 ± 0.022 AU (Rosenthal et al. 2021).We refer to this planet as HD 114729 Ab.This system also hosts an M dwarf companion: Mugrauer et al. (2005) directly imaged HD 114729 and found a bound companion at a projected separation of 282 AU.Using the stellar cooling models of Baraffe et al. (1998), as well as Butler et al. (2003)'s chromospherically derived age of 6 Gyr, Mugrauer et al. (2005) estimated HD 114729 B's mass to be 0.253 M ⊙ .
We isolated a subset of the RV baseline covering approximately four years, divided it into four equal slices, and fit a second order polynomial to the RVs in each slice.The time series slices and polynomial fits are shown in the top of Figure 6.
We used ethraid to derive m c -a posteriors for each slice, incorporating both the fitted trend and curvature as well as HGCA astrometry (0.11±0.03 mas/yr), which is the same for all slices.In all cases, ethraid's derived posterior distributions are discrepant with the true mass and semi-major axis of HD 114729 Ab at > 2σ (Figure 6, bottom).This disagreement is driven by the astrometric data.The RV variability in the one-year slices is dominated by HD 114729 Ab because of its relatively short ∼3-year orbital period, and indeed the RV-only posteriors agree more closely with that planet's true parameters.Meanwhile, the astrometric variability reported in the HGCA is measured over 25 years, and thus includes a strong contribution from HD 114729 B's 4400-year orbit.We also plot this companion's parameters in Figure 6 to demonstrate its consistency with the astrometry data.
Our analysis of HD 114729 illustrates how ethraid's fits may be misleading in the case of multiple companions.It is very valuable to provide direct imaging constraints whenever available to mitigate this failure mode.6.4.Case Study: HD 12661's two planets are misinterpreted as a single companion HD 12661 is an inactive (log R ′ HK =-5.06), metal-rich ([Fe/H]=0.293)G6 V star with two gas giant compan-ions (Fischer et al. 2001(Fischer et al. , 2003)).HD 12661 b has a minimum mass of m c sin i = 2.283 +0.062 −0.063 M J and a separation of a = 0.824 ± 0.011 AU, while HD 12661 c has m c sin i = 1.855 ± 0.054 M J and a = 2.86 +0.038 −0.039 AU (Rosenthal et al. 2021).
We chose a ∼2.5-year subset of the RVs and divided it into three slices, each covering the same time interval.We fit linear and quadratic terms to each slice.The RV subset and slices are shown in the upper panels of Figure 7.
We used the trend and curvature derived for each slice in an ethraid partial orbit fit.We also included the astrometric trend of 0.23 ± 0.05 mas/yr.The orbit posteriors are shown in the lower panels of Figure 7, along with the the orbital parameters of HD 12661 b and c. ethraid's RV posteriors are consistent with planet b's parameters, which is expected because with a larger mass than planet c and a period of 264 days, this planet dominates the RV variability over the 200day slices.Furthermore, this consistency is not sensitive to the phase or curvature of the RV slices.The slices in Figure 7 exhibit linear, quadratic, and higher order curvature, but are all nevertheless consistent with HD 12661 b's parameters in the lower panels.
The astrometry posterior, on the other hand, is consistent with neither planet's parameters.We modeled the combined astrometric signature of both planets and found that planet c dominates this signal due to its wider separation.Furthermore, we found that if planet c's inclination is within 25 degrees of 0 (i.e., a "face-on" orbit), then the resulting underestimate of planet c's mass (a factor of ≳ 2.4) may account for the disparity between this system's measured and expected astrometric signatures.A third, more distant companion could also explain this signal, though we did not explore this possibility in detail.
Like HD 114729, HD 12661 shows the limitations of ethraid's approach to treating RV/astrometric trends.While planet b dominated the RV trend/curvature because of its short period, planet c had a much greater influence on the astrometric signal.ethraid's assumption that both signals were due to one companion produced constraints inconsistent with either planet.We caution that unless the presence of multiple giant planets can be ruled out, ethraid's single companion assumption may produce false parameter estimates.2023): a = 0.98 AU, mc sin i = 1.05 MJ.Without astrometry, the trend and imaging data provide too little information to constrain the planet's mass and separation.However, they are able to rule out massive stars as well as the (a, mc) pairs in the white regions.Note also that the posteriors in panels e) and f ) favor shorter-period models due to the higher RV curvature associated with their trends.
2. Dynamic likelihood storage to mitigate RAM overload.
3. Analytical removal of the RV and astrometric signatures of known stellar companions, improving ethraid's predictive power in systems like HD 114729 (Section 6.3).
4. Improved companion mass prior to reflect the relative prevalences of planets, brown dwarfs, and low-mass stellar companions.
5. Inclusion of brown dwarf cooling models for systems of ages other than 5 Gyr.
We compiled the above improvements based primarily on our own use of ethraid, but we welcome contributions from the exoplanet community to make ethraid a useful and accessible tool.

CONCLUSIONS
In this work we presented an open-source Python package for constraining the masses and separations of long-period companions using RVs, astrometry,  RV signature due to its short period, as shown by its consistency with the green RV posteriors in each slice.Meanwhile, planet c is responsible for the majority of the astrometric signal.The inconsistency between these planets' measured parameters and the combined posterior (red) shows that ethraid misinterpreted these separate contributions as originating from a single object of higher mass.and direct imaging.ethraid uses simple data inputs to calculate RV and astrometric posteriors, requiring only two values and their errors for the first and one value and its error for the second.Due to the limited information content of this data, ethraid is designed to broadly constrain companion parameters rather than to produce tight orbital fits, giving a 'first look' at the nature of a companion and helping to disambiguate between multiple competing hypotheses (e.g., a close-in planet vs. a distant M dwarf).ethraid is pip-installable and may be downloaded from GitHub.
We demonstrated ethraid's strengths and weaknesses by testing it on four systems hosting one or more known companions.HD 114729 (Butler et al. 2003) hosts a giant planet and a stellar companion, and HD 12661 (Fischer et al. 2001(Fischer et al. , 2003) ) hosts two gas giant planets.In both cases, ethraid's assumption of a single companion led to spurious parameter predictions.HD 117207 and TOI-1694 each host only one distant giant planet.ethraid consistently recovered HD 117207 b's parameters, and although a lack of astrometric data prevented a recovery of TOI-1694 b, our combined RV-direct imaging analysis ruled out high-mass stellar companions and high-separation (≳ 10 AU) planets and brown dwarfs.ethraid leverages up to three exoplanet detection techniques to constrain companion properties given minimal observational data.It is intended to help users evaluate the merit of potential companion models, as opposed to giving precise parameter constraints.We caution that due to ethraid's implicit assumption that any input signals are caused by a single companion, applying it to multi-companion systems may give nonphysical results.Users may run ethraid on a system using one configuration file.We plan to maintain and improve ethraid so that it may be of use to the community.
Software: radvel (Fulton et al. 2018) Figure 8. Left: Panel d) from Figure 5, showing the companion models consistent with a subset of the RVs of TOI-1694, as well as the imaging data.Note that despite falling above the gray 95% contour, some high-mass models are consistent with the imaging data because their orbital geometries result in small projected angular separations from the host star.Right: Parameter constraints for the same RV and imaging data, but using the approximate imaging calculation.In the approximate case, the imaging contour represents a hard bound, above which companions are ruled out with certainty.This causes the 1σ upper bound on mass (dark red) to be underestimated by ∼ 30% at 16 AU, while being roughly unaffected at 8 AU.In the long-period regime, using the small angle approximation for θHG is appropriate.Right: A similar diagram for short periods (≲ 2 AU).Arrows t1 and t2 show the position of the host star at the beginning and end of the Gaia EDR3 data collection window.We chose t2's position arbitrarily to illustrate that for short periods, the angle between the two position vectors is not small, but rather can take any value on [0, 2π].We do not show the analogous vectors for ⃗ µHG because THG is nearly nine times the duration of TG; this larger denominator drives ⃗ µHG to nearly zero for short-period orbits.
We begin with the long-period regime, that is, separations corresponding to orbital periods much longer than T HG , the 24.75-year interval between the Hipparcos and Gaia characteristic epochs.From Equation B11 we see that ∆µ depends only on the anomalous proper motions, so we may neglect center-of-mass proper motion in this derivation.
Because the companion orbit is circular and face-on, the proper motion vector has the same magnitude at all epochs, that is, d⋆P .Furthermore, in the long-period regime, the average proper motion approximates the instantaneous value: |⃗ µ HG | ∼ |⃗ µ G |.The proper motion components at the Gaia and HG epochs are thus where θ HG is the position angle of the host star at the HG epoch (see Figure 10, left panel), given by From these proper motions, we may compute ∆µ: where we approximated that sin θ HG ∼ θ HG and cos θ HG ∼ 1 − θ 2 HG 2 .We further approximate that the x-component of ∆⃗ µ dominates in the long-period regime, yielding: Thus, m c ∝ a 2 in the long-period regime.This behavior is shown for large a values (a ≳ 20 AU) in the left panel of Figure 11.The short-period portion of the astrometry posterior surface displays two distinct features: an overall negative slope and a series of spike-like formations extending to high mass.We derive the general linear character of the astrometry posterior in the short-period regime using the same picture we used to understand the long-period behavior.First examining ⃗ µ G , we have from Equation B10: where we have applied our assumptions of a face-on circular orbit and converted velocity to proper motion.Unlike in the long-period regime, T G is comparable to or multiple times the duration of orbits in the short-period regime.As a result, M (t 1 ) and M (t 2 ) may differ by any angle from 0 to 2π.To continue, we set M (t 1 ) = 0 without loss of generality, and average over the possible values of M (t 2 ): where ⟨ ⟩ denotes the average over the interval [0, 2π].We may apply this same procedure to ⃗ µ HG , but note that in this case, the denominator of the right-hand side of Equation B17 is ∼ 25 years, as opposed to ≲ 3 years for ⃗ µ G .Thus, ⃗ µ G dominates in the short-period regime, and we neglect ⃗ µ HG .See the right panel of Figure 10 for a pictorial representation of this scenario.Finally, we use the relation between a and a ⋆ to find that m c ∝ a −1 for short periods: We show this relation in Figure 11.
The spikes represent harmonics of T G , the 34-month Gaia EDR3 mission baseline.For periods shorter than the ∼ 25-year Hipparcos-Gaia baseline, averaging the orbital position over that interval produces a cancellation effect, reducing the magnitude of ⃗ µ HG so that ⃗ µ G dominates ∆µ in the short-period regime.For orbital periods which divide T G , i.e.T G itself, T G 2 , etc., ⃗ µ G → 0 as well because the host star returns to its original position over the course of the Gaia mission, resulting in an average measured velocity of zero.The spikes in the right panel of Figure 11 demonstrate that model orbits must have higher companion mass near the harmonics to produce a constant value of ∆µ.We observe a similar effect at harmonics of T HG (left panel of Figure 11), though these spikes are less pronounced because both ⃗ µ HG and ⃗ µ G contribute to ∆µ in that regime.
C. APPENDIX: ROTATION MATRIX R is the rotation matrix P 3 P 2 P 1 , described in Equations 2.119-2.121 of Murray & Dermott (2010).It is not given explicitly in that text.We exploit a symmetry in Sections 3.2 and 3.3 which warrants explanation here.In those sections, we apply R to the stellar position and velocity vectors, despite the fact that i, ω, and Ω correspond to the companion orbit.Ω and i are the same for the star and the companion, but the stellar argument of periapsis is offset from that of the companion by π radians: ω ⋆ = ω c + π.This means that rotating the stellar position and velocity vectors using ω c is not correct in general.We demonstrate below that for our calculations, this asymmetry vanishes.
First, we relabel the matrix R above as R c to denote its use of the sampled ω c .We also define the matrix R ⋆ , which results from the substitution ω c → ω ⋆ : where we have abbreviated the matrix elements for clarity.We summarize the rotation steps of Section 3.2 as The situation is analogous for ⃗ v obs , as well as for ⃗ ρ in Section 3.3.Thus, the use of R c instead of R ⋆ introduces a single negative sign to each rotated vector.Because all vectors are equally affected, and in all cases we ultimately work with the vector moduli, this negative sign has no effect on our likelihood calculation.

Figure 1 .
Figure 1.ethraid's fitting algorithm.The steps are the same for RVs (left), astrometry (center), and direct imaging (right): ethraid samples a set of model orbital parameters, forward models the data-specific parameters, evaluates the model likelihood, and finally marginalizes many such likelihoods to derive the posterior PDF.Arrows represent inputs used in a given step: blue arrows show inputs from the previous step, red arrows show inputs provided by the user, and purple arrows show inputs that are fixed in the code.The bottom row shows three of the plots available to the user through the plot command.Left: The 2D mc-a posteriors for each of the RV (green), astrometry (blue), and imaging (gray line) data sets, plus the combined posterior (red) conditioned on all three.The dark/light regions of each posterior show 68/95% confidence intervals.For imaging, the gray contour marks the 95% confidence boundary, with models below/above the contour ruled in/out by the data.Right: CDFs of the marginalized a and mc posteriors.The vertical red lines indicate 95% confidence intervals.ethraid also produces corresponding PDFs by default.

Figure 2 .
Figure 2. Our model of astrometric proper motion anomaly.Red ellipses show the stellar orbit associated with a nonluminous bound companion (not shown).Green arrows show the average proper motion vectors fitted to the data from Hipparcos (left) and Gaia (right).Gold stars show the average positions during Gaia and Hipparcos, and the blue arrow shows ⃗ µHG, the average proper motion between them, calculated as the difference in average positions divided by the difference between the mission midpoints.The proper motion anomaly ∆µ, shown in red, quantifies the difference between the average Gaia proper motion and the positionderived proper motion.Figure adapted from Kervella et al. (2019).

Figure 3 .
Figure3.ethraid's likelihood calculation for direct imaging.Beginning with a set of sampled parameters, ethraid models angular separation (either exactly or approximately) using the orbital elements, and contrast using the companion mass, together with a linear interpolation function derived from Table5ofPecaut & Mamajek (2013) and Table4ofBaraffe et al. (2003), denoted by 'PM13' and 'B03,' respectively.ethraid then determines whether the model companion's brightness would have exceeded the minimum detectability threshold at the model separation in the real imaging data (gold line) and assigns a likelihood of 0 for detectable companions and 1 for non-detectable companions.The example model in the diagram falls above the threshold, so we assign it a likelihood of 0.
3 https://www.h5py.org/set of posterior arrays from the specified file path and optionally produces and saves three figures.The first includes two one-dimensional marginalized posterior PDFs side by side, one for a and the other for m c .The second is analogous to the first, but shows CDFs instead of PDFs.The third is a two-dimensional plot in m c -a space displaying the posteriors conditioned on each data set individually as well as the combined posterior.We show examples of the one-dimensional CDFs and two dimensional PDFs in in Figure 1.

Figure 4 .
Figure 4.The top panel shows a 7.5-year subset of HD 117207's RV time series, divided into four slices based on the epoch in which the RVs were measured.Panels a)-d) show the RVs from each slice, along with the linear/quadratic trend fit to those RVs using radvel.We show the measurement error of each RV with horizontal black lines.We used ethraid to model each trend/curvature pair, along with HGCA astrometry, in panels e)-h), respectively.The gold star in each panel shows HD 117207 b's semi-major axis and mc sin i values measured by long-baseline RV surveys: a = 3.744 AU, mc sin i = 1.87 MJ, where we have approximated a conversion from mc sin i to mc by dividing by the median sin i value of 0.866.The overlap between ethraid's predicted parameters and the true values in all four panels demonstrates ethraid's reliability in estimating companion parameters over a range of orbital phases.
7. FUTURE IMPROVEMENTSDevelopment of ethraid is ongoing, and we plan to implement a number of improvements in future versions, including:1.Parallelization of likelihood calculations.

Figure 5 .
Figure 5. Same as Figure 4 for TOI-1694's two year baseline.The gray lines in panels d)-f ) show posteriors derived from direct imaging.The gold stars show TOI-1694 b's median semi-major axis and adjusted mc sin i values measured by Van Zandt et al. (2023): a = 0.98 AU, mc sin i = 1.05 MJ.Without astrometry, the trend and imaging data provide too little information to constrain the planet's mass and separation.However, they are able to rule out massive stars as well as the (a, mc) pairs in the white regions.Note also that the posteriors in panels e) and f ) favor shorter-period models due to the higher RV curvature associated with their trends.

Figure 6 .
Figure 6.Same as Figure 4 for HD 114729.The gold star at the lower left of each panel shows the semi-major axis and adjusted mc sin i of HD 114729 Ab measured by long-baseline RV surveys: 2.094 AU, mc sin i = 0.892 MJ.The upper right star shows the same for HD 114729 B, with a projected separation and mass of 282 AU and M = 0.253 M⊙.We estimated HD 114729 B's true separation by dividing the projected separation of 282 AU by 0.79, as described in Section 3.3.Agreement between the planet parameters and the RV posteriors indicates that ethraid is correctly recovering HD 114729 Ab's RV signature.Meanwhile, the consistency between HD 114729 B's parameters and the astrometric posterior suggests that HD 114729 B produced the astrometric trend.Note that in panels a) and d), the fitted trend/curv values had precisions < 2σ, leading to broad RV posteriors in panels e) and h).

Figure 7 .
Figure 7. Same as Figure 4 for HD 12661.For this system, we include RVs from both the Keck/HIRES and APF/Levy instruments.The gold stars show the a and adjusted mc sin i values of HD 12661 b (left) and c (right).Planet b dominates theRV signature due to its short period, as shown by its consistency with the green RV posteriors in each slice.Meanwhile, planet c is responsible for the majority of the astrometric signal.The inconsistency between these planets' measured parameters and the combined posterior (red) shows that ethraid misinterpreted these separate contributions as originating from a single object of higher mass.

Figure 10 .
Figure 10.Left: A diagram of the circular, face-on orbit we use to understand the behavior of the astrometry posterior.The black circle shows the orbit of the host star about the system barycenter due to a companion at wide separation (≳ 20 AU).Solid and dotted arrows show the average stellar position and proper motion vectors, respectively, in each of the Hipparcos, HG, and Gaia EDR3 epochs.θHG gives the position angle at the HG epoch relative to the position angle during the Hipparcos epoch.In the long-period regime, using the small angle approximation for θHG is appropriate.Right: A similar diagram for short periods (≲ 2 AU).Arrows t1 and t2 show the position of the host star at the beginning and end of the Gaia EDR3 data collection window.We chose t2's position arbitrarily to illustrate that for short periods, the angle between the two position vectors is not small, but rather can take any value on [0, 2π].We do not show the analogous vectors for ⃗ µHG because THG is nearly nine times the duration of TG; this larger denominator drives ⃗ µHG to nearly zero for short-period orbits.

Figure 11 .
Figure 11.Left: The posterior surface calculated from astrometric data.The red line follows the relationship mc ∝ a −1 , while the green line follows mc ∝ a 2 .Right: The short-period regime of the surface at left.Iso-period lines are shown in red: the right-most line traces P = TG ∼ 2.83 yr, the next follows P = T G2 , then P = T G 3 , and so on.High-likelihood models near these periods must have large companion masses to counter the small net displacement of the star.The same effect can be seen at harmonics of the ∼ 25-yr Hipparcos-Gaia baseline (∼ 4 − 16 AU), though it is less pronounced because ⃗ µG remains nonzero in this regime even when ⃗ µHG vanishes.

Table 1 .
Configuration file parameters gammaddot Float -Quadratic RV curvature (m/s/day/day); setting to None excludes curvature gammaddot_err Float -Error on gammaddot (m/s/day/day); setting to None excludes curvature

Table 2 .
Command line interface parameters