The Orbital Eccentricities of Directly Imaged Companions Using Observable-Based Priors: Implications for Population-level Distributions

The eccentricity of a sub-stellar companion is an important tracer of its formation history. Directly imaged companions often present poorly constrained eccentricities. A recently developed prior framework for orbit fitting called ''observable-based priors'' has the advantage of improving biases in derived orbit parameters for objects with minimal phase coverage, which is the case for the majority of directly imaged companions. We use observable-based priors to fit the orbits of 21 exoplanets and brown dwarfs in an effort to obtain the eccentricity distributions with minimized biases. We present the objects' individual posteriors compared to their previously derived distributions, showing in many cases a shift toward lower eccentricities. We analyze the companions' eccentricity distribution at a population level, and compare this to the distributions obtained with the traditional uniform priors. We fit a Beta distribution to our posteriors using observable-based priors, obtaining shape parameters $\alpha = 1.09^{+0.30}_{-0.22}$ and $\beta = 1.42^{+0.33}_{-0.25}$. This represents an approximately flat distribution of eccentricities. The derived $\alpha$ and $\beta$ parameters are consistent with the values obtained using uniform priors, though uniform priors lead to a tail at high eccentricities. We find that separating the population into high and low mass companions yields different distributions depending on the classification of intermediate mass objects. We also determine via simulation that the minimal orbit coverage needed to give meaningful posteriors under the assumptions made for directly imaged planets is $\approx$ 15% of the inferred period of the orbit.


INTRODUCTION
The field of exoplanet direct imaging has advanced significantly in recent years, with the detection capabilities of ground based telescopes improving due to technological developments in adaptive optics (AO). Direct imaging also allows for the measurement of the relative astrometry of planets over time, which can in turn allow astronomers to better constrain the orbital parameters of these objects and the overall orbital architecture of their systems. By better understanding the architecture of these extrasolar systems, we may obtain an understanding of formation pathways of planetary systems in our Galaxy. Thus, orbit fitting of directly imaged exoplanets has the potential to shed light on the early formation history and dynamical evolution of planetary systems with widely-separated, gas giant planets.
In particular, the eccentricity of a planet can tell us much about how it was formed. In the theory of protoplanetary disk formation, disk fluid elements can initially be in an eccentric orbit, but soon lose energy through collisions and settle into the minimum energy orbit, which is a circular orbit (Lodato 2008). Therefore, planets forming in a protoplanetary disk via core or pebble accretion should have lower eccentricities, but can develop higher eccentricities through processes such as planet-planet scattering or disk migration , Dong & Dawson 2016, Papaloizou & Terquem 2006, Chatterjee et al. 2008, Johansen & Lambrechts 2017. Various protoplanetary disk processes can also either damp eccentricities (e.g. Kley & Nelson 2012;Bitsch & Kley 2010) or excite eccentricities (e.g. Moorhead & Adams 2008;Goldreich & Sari 2003). Formation of planets via gravitational instability can potentially form planets with higher eccentricities. Generally, gravitational instability planet formation occurs further from the host star, causing eccentricity damping timescales to be longer at such large separations (e.g. Mayer et al. 2004). These various processes may lead to a variety of eccentricities for individual planetary systems, so it is difficult to tell which mechanism is operating for individual systems. However, the eccentricity distribution of many planetary systems (i.e., the parent distribution) can give us strong constraints on which formation mechanisms are at work. Planets that slowly form from pro-toplanetary disks in unperturbed orbits should have low eccentricities, while planets that undergo planetary migration or outward scattering can have a range of varying eccentricities.
Efforts to obtain an eccentricity distribution at a population-level for sub-stellar companions have already been made for a number of exoplanet and brown dwarf populations. Hogg et al. 2010 presents the method of hierarchical Bayesian modeling to obtain parameters for the underlying parent eccentricity distribution from the posteriors of exoplanet parameters given from orbital fits. The planet population simulated by Hogg et al. 2010 was inferred from radial velocity measurements of hot Jupiters. Kipping 2013 has constrained eccentricity distributions of short-period exoplanets using data obtained from 396 companions using the radial velocity (RV) method. Bowler et al. 2020 used data from 27 long period, directly imaged companions to obtain an eccentricity distribution of these objects. They find that extrasolar companions present different eccentricity distributions if they are low mass or high mass. They conclude that "exoplanets", or lower mass companions, have lower eccentricities (so they likely formed from a disk) and "brown dwarfs", or higher mass companions, have higher eccentricities (consistent with binary star formation). This conclusion, if true, is groundbreaking -it suggests that these companions can be observationally distinguished via their eccentricities. In this work, we revisit this analysis using a different set of priors and modified sample of directly imaged substellar companions.
In order to determine the orbital eccentricities of exoplanets, the standard procedure is to fit astrometric data of the companion relative to the star at different points in time using Bayesian statistics. Within the Bayesian framework, prior probability distributions (priors) for each orbit parameter must be assumed to ultimately infer posterior probability distributions (posteriors) for each parameter. The most common priors assumed for orbital parameters are model-based priors (i.e., they assume a prior distribution in the 6 model parameters of a Keplerian orbit). They generally assume distributions that are uniform (for the eccentricity, argument of periastron and period of ascending nodes), log-uniform (semi-major axis) and sine (inclination). It is also common practice to narrow down the parameter space with physically motivated priors, such as stability constraints, coplanarity and stellar rotation rates (e.g. Pearce et al. 2015;Wang et al. 2018;Wang et al. 2020a; Thompson et al. 2023;Yimiao Zhang et al. 2023). However, for orbits where less than 40% of the orbital arc is covered, which is the majority of directly imaged companions, this standard method of model-based priors (which herein we will refer to as "uniform" priors, as is done by O'Neil et al. 2019) has been shown to introduce biases in the resulting posterior distributions, generating inaccurate parameters and confidence intervals (Martinez et al. 2017, Lucy 2014. The eccentricity parameter is generally affected by this issue, with many objects presenting a bias towards high eccentricities when fit with uniform priors. O'Neil et al. 2019 presented a new approach to priors for the orbit fitting of long period resolved companions, which is based on uniformity in the observable parameters rather than in the model parameters. This approach is called "observable-based" priors. This method has been shown to reduce this bias in orbital parameters where the orbital arc of the object is less than 40% covered. In this work, we aim to obtain the eccentricity distribution of extrasolar companions at a population level, using observable-based priors to fit the orbits of 21 directly imaged companions. Given the majority of our sub-stellar companion sample is under-sampled (i.e., their current orbital coverage spans a small amount of their orbital arc), we also perform a series of simulations to assess the minimum orbital coverage needed in order to get meaningful eccentricity posteriors for directly imaged companions.
Our analysis is outlined as follows. We introduce our sample and present the astrometric and radial velocity data used for our orbit fits in sections 2.1 and 2.2. In section 3.1, we present the concept of observable-based priors. We present our results in section 4, describing the results of our orbit fitting of the 21 companions (section 4.1) and simulations for the minimal orbital coverage needed (section 4.2). We discuss the implications of our work in 5 and our conclusions in 6.

Sample Selection and Astrometric Data
The companions chosen are a sub-sample of objects from Bowler et al. 2020. The criteria for choosing these objects is outlined in detail in Bowler et al. 2020; we summarize them here: the objects must be at a projected separation of 5-100 AU from the host star at the time of discovery, and the hosts must be stars (M > 75 M Jup ), such that any companions that are clearly part of binary systems are excluded from the sample. For the systems with more than one detected planet, such as HR 8799, PDS 70, HD 206893 and β Pictoris, we chose only one planet as the "representative" of the system's eccentricity in this analysis. We made this choice so that multiple planets from the same system would not bias our resulting eccentricity distributions, since it has been shown that planets with multiple systems present lower eccentricities due to stability requirements ), and therefore the eccentricities within a single system are correlated. We also do not include any companions that have less than 4 epochs of observation, as our software requires this value as a minimum amount of astrometry for orbit fitting given the number of free parameters. With these requirements, 21 companions are in the sample used in this work. The sample is presented on Table 1.

Radial Velocity Data
Relative radial velocities (RV), when available, contribute a great deal of information to visual orbit fits. Since RV measurements are challenging to obtain for directly imaged companions, only six objects in this study have radial velocities available in the literature. The RVs used are shown in Table 3.
Using RVs for these six objects allows us to eliminate the degeneracy in Ω and therefore better constrain the orbital plane of the object. Additionally, using RV data for HD 1160 b changed the eccentricity distribution of the companion when using uniform priors (see Section 4.1.2).

Orbit Fitting
We use observable-based priors for fitting the relative astrometry to orbits for all 21 directly imaged companions. These priors are implemented in the orbit fitting software Efit5 (Meyer et al. 2012). Efit5 uses MULTI-NEST (Feroz & Hobson 2008;Feroz et al. 2009), a multimodal (or nested) sampling algorithm, to perform a Bayesian analysis on the data. For all of these fits, we use 3,000 live points in the nested sampling algorithm. We also include relative radial velocities (RVs) when available. In order to compare these results to the uniform prior approach, we also fit for the orbits of these companions with uniform priors. For our orbit fitting with uniform priors, we use both the Efit5 package and the orbitize! package ) with the Markov-Chain Monte Carlo (MCMC) approach. We verify that both fitting methods with uniform priors give consistent results. For the uniform priors with orbitize, we fit for 50 million possible orbits, using 1,000 walkers, 20 temperatures, 20 threads, 10,000 burn steps and a thin factor of 10, as is done in Bowler et al. 2020. Figure 1. The weighted eccentricity distribution derived via our orbital fits for the sample of 21 companions using observablebased priors. The x-axis represents the eccentricity values and the y-axis represents the normalized incidence of eccentricities. Here, the fits include both astrometric and radial velocity data points .
After fitting the orbits of the 21 companions using observable-based priors, we use their posteriors to fit for the eccentricity distribution of our entire population of 21 companions. We use a Beta distribution as the model for our parent eccentricity distribution, as is done in other works (e.g. Bowler et al. 2020;Kipping 2013). In order to recover the parent distribution from the sample's posteriors, we use Maximum Likelihood Estimate fitter in a bootstrapping manner. The details of our recovery technique are presented in Sections 3.2 and 3.3. We then split the sample into subsamples to assess possible distinctions between the population distributions.
The purpose of observable-based priors is to improve orbital estimates for orbits where the data covers a low percentage of the orbital arc (< 40%). Here, we briefly summarize the formulation of observable-based priors. A detailed formulation is outlined in O'Neil et al. 2019.
In its general form, observable-priors assume that all regions of observable parameter space that can be observed are equally likely, thus emphasizing uniformity in observables rather than in model parameters. In the case of orbit fitting, our fit starts with measured observables D from the astrometry: or, for astrometry and radial velocities: where x and y are the object's positions ( right ascension (R.A.) and declination (Dec)) in the plane of the sky relative to the position of the primary (x o and y o ) and v z is the velocity relative to the star. These measured observables are linearly related to the orbital observables (which describe the position and motion in the orbital plane) by the Thiele-Innes constants (e.g Hartkopf et al. 1989, Wright & Howard 2009). Due to this linear relationship, a uniform distribution in the measured observables would imply a uniform distribution in the orbital observables. The orbital observables, denoted here as X, Y, V x and V y , are also connected to the model parameters according to the following equations (e.g. Hilditch 2001; Ghez et al. 2003): where G is the gravitational constant, E is the eccentric anomaly, e is the eccentricity, a is the semimajor axis and M is the mass of the system. By transforming between measured observables and orbital observables, and then between orbital observables and model parameters, we can transform between measured observables and model parameters. This allows us to express a distribution that is uniform in the measured observables in terms of model parameters. This fitting method reduces biases in orbital parameters by a factor of two in orbits with low phase coverage (O'Neil et al. 2019). For the eccentricity parameter, observable-based priors reduce the bias towards artificially high eccentricities usually found in fits with flat priors. Flat priors lead to a biased region in periastron passage (T o ) parameter space, where the T o tends to artificially coincide with the observation epochs ). This bias is mitigated in observable-based priors because they suppress this biased region of the parameter space when sampling it. Thus, observable-based priors present a suitable fitting method for the orbital parameters of our object sample, which all have phase coverage <40%, and provide a different view on the eccentricity distribution of the companions when compared to the standard method for orbit fitting.

Recovering the Population Eccentricity Distribution
The individual eccentricity posterior distributions obtained with the observable-based prior method and the data described in the previous sections are presented on Figure 1. From these individual distributions, we now seek to determine the eccentricity distribution of the population. We adopt a Beta distribution for our underlying parent distribution for exoplanet eccentricity, similar to Bowler et al. 2020, Hogg et al. 2010, Kipping 2013and Van Eylen et al. 2019. The Beta distribution is a convenient choice for this study because its values, like the eccentricities, range from 0 to 1. It also only requires two positive parameters, α and β, which allow it to present a variety of shapes, such as linear, Gaussian and uniform. The Beta distribution is represented by: Here, Γ is the the Gamma function. We use the posteriors in our sample and the beta function fitter present on the SciPy package (Virtanen et al. 2020), which uses the Maximum Likelihood Estimate (MLE) method, for fitting our posteriors to the distribution. We perform our fitting using a bootstrapping method by repeating the following steps: (1) sampling an eccentricity point from each posterior chain and then (2) fitting a Beta distribution to those points using the SciPy Beta fitting functionality.

Beta Distribution Recovery Simulation
In order to validate the choice of using the Maximum Likelihood Estimate fitter for obtaining our estimated α and β distributions, we develop a procedure similar to the recovery analysis presented in Bowler et al. 2020. The purpose of this exercise is to test our bootstrapping method. We obtain points from the α = 0.867, β = 3.03 distribution presented by Kipping 2013 for warm Jupiters and points from the α = 0.95, β = 1.30 distribution presented by Bowler et al. 2020. We vary the number of points taken, taking N = 10, 20, 50 and 100. N represents the number of companions that contribute to our fitting sample. Our resulting recovered distributions are presented in Tables 4 (Kipping) and 5 (Bowler). In the plots with Beta distribution PDFs (Figure 2), we plot the "true" input distribution in a darker blue/red shade along with 100 randomly sampled distributions from within the confidence intervals in a lighter blue/red shade.  Table 4 continued

RESULTS
In this section, we show our results for the eccentricity distributions of individual objects and different populations of directly imaged companions. We find that some objects presented significant eccentricity shifts when using observable-based priors instead of uniform priors (4.1.2). For our full sample of 21 objects, we obtain shape parameters α = 1.09 +0.30 −0.22 and β = 1.42 +0.33 −0.25 (4.1.1). We test whether leaving one object out of the sample changes our results (4.1.3), and find that our sample does not have any true outliers. We also test splitting the sample by companion mass (low and high), and find that the consistency between the distributions changes depending on where intermediate mass objects are placed (4.1.4). Finally, we simulate how much orbital coverage is needed to obtain reliable posteriors in our distributions, both for uniform and observable priors (4.2). We obtain 15% of the period of the orbit as a result.

Population Eccentricity Distribution
We use the posteriors from 21 companions to investigate changes in the inferred eccentricity of individual systems due to observable-based priors, as well as to obtain an inferred eccentricity parent distribution. Following the procedure from section 3.3, we randomly sample the posteriors to obtain a possible range of α and β shape parameters of the parent Beta distribution.
Our result for the eccentricity distribution of the entire sample is presented in Figure 3. We obtain a beta distribution with α = 1.09 +0.30 −0.22 and β = 1.42 +0.33 −0.25 . This shape presents a near-uniform distribution for the eccentricities of substellar companions, with a slight tendency towards lower eccentricities.
In order to determine if two derived population distributions are consistent with each other, we estimate the consistency of their respective parameters α, and β. The procedure defined here is a general method, which we use to determine consistency of two Beta probability distributions throughout this work. We define two random variables as the differences δα = α 2 − α 1 and δβ = β 2 − β 1 where α 1 , β 1 defines the first population, and α 2 , β 2 defines the second population. If the two populations are consistent, the distributions of both δα and δβ should be consistent with zero. The distributions of δα and δβ are calculated by sampling the respective distributions of α 1 , α 2 , β 1 , and β 2 in a bootstrapping manner from the resulting fits (i.e., sampling the α 1 , β 1 and α 2 , β 2 pairs from the fit chains), and calculating their differences at every iteration. We then calculate the p-value of the sample: by integrating the 2D distribution of (δα,δβ) inside the contour of constant density that goes through the origin. The integration is performed using a Gaussian kernel density estimator. An example of this process along with more details on how P is calculated is shown in Appendix B. We alternatively calculate P using the Kolmogorov-Smirnov Test (KS Test) using in the package ndtest in Python (github.com/syrte/ndtest). For this test, the discrete values of x and y for the PDF must be used, so we sample 100 pairs of α and β from the chains of two distributions, obtain 1 value from each of the pairs and compare the two using the KS Test. We obtain fully consistent results between the kernel density estimate and the KS test, and thus report P from the kernel density herein.
For the distributions presented in Bowler et al. 2020 and Kipping 2013, we do not have the posterior chains of α and β pairs so we randomly sample from uniform distributions with the α and β ranges (including the uncertainties) given by these works. When comparing our Beta distribution to the distributions found by Bowler et al. 2020 andKipping 2013, we find that our distribution is consistent with Bowler et al. 2020's (P of 0.77) but not with Kipping 2013's (P of 0.000) (Figure 4).
The reasons for the inconsistency with Kipping 2013 could involve the fact that Kipping 2013's distribution is for short-period exoplanets (therefore not including the long period objects that populate our sample), the lower N in our sample (21 vs. 396 objects) and the difference in detection method for the objects. Kipping 2013's sample used RVs as the primary detection method of the companions, while we used direct imaging. Kipping 2013 also has much smaller confidence intervals than both our shape parameters and Bowler et al. 2020's shape parameters. This, again, is likely due to a small confidence interval in the individual eccentricities as well as the limited sample size (< 10%) of the direct imaging sample when compared to the RV sample.
As mentioned in section 3.1, the traditional prior framework often presents a bias towards higher eccentricities. This is why, when comparing our distribution to the one obtained by Bowler et al. 2020, we note that our parent distribution presents a reduced curved peak near high eccentricities, favoring slightly lower eccentricities. This could indicate that the bias towards higher eccentricities in individual orbits can also provide a bias towards higher eccentricities in the underlying parent distribution, particularly when the sample size is small. The observable-based prior aims to avoid such a bias in individual orbit fits; however, it is not completely capable of eliminating it. We also note that our loweccentricity end presents a smaller peak (near e = 0), likely due to the fact that we did not include more than one planet from multi-planet systems (e.g. HR 8799 b,d, and e are not included) into our sample. Including those objects could cause an artificial increase in the PDF at low eccentricities since their eccentricities are correlated by being in the same system and being stable.

Individual Objects with Significant Eccentricity
Shifts Some objects present a significant shift in their eccentricity distributions when comparing orbit fits from uniform priors and observable-based priors. The orbits with significant eccentricity shifts in our sample are HD 49197 b, HR 2562 b, HIP 65426 b, HD 1160 b and PZ Tel b. These 5 objects compose more than 20% of our companion sample.
HD 49197 b had its eccentricity posterior change from 0.62 +0.33 −0.43 to 0.28 +0.29 −0.20 , (weighted median of the posterior distribution (50th percentile)), with the uncertainties encompassing the 68th percentile values, favoring lower eccentricity solutions when fit with observable- Figure 4. Comparison of the obtained beta distribution for the entire sample with the distributions of Bowler et al. 2020 andKipping 2013. We note that the distribution obtained by Kipping 2013 is for short period warm Jupiters, whereas Bowler et al. 2020's distribution is for long period extrasolar planets and brown dwarfs. Our sample closely resembles the one used by the latter's work, but with a lower high eccentricity incidence -likely due to the change from uniform to observable priors. based priors. HR 2562 b had its eccentricity posterior change from 0.49 +0.44 −0.36 to 0.66 +0.16 −0.33 , changing from a bimodal distribution with eccentricity modes at ≈ 0.2 and 0.9 to a fit that spans a large amount of eccentricity values but disfavors eccentricities > 0.9. HIP 65426 b had its posteriors change from 0.57 +0.29 −0.38 to 0.26 +0.41 −0.20 , favoring lower eccentricity values with the observable-prior fit.
One object that is of particular interest is HD 1160 b, one of two known companions orbiting HD 1160. When fitting the orbit of HD 1160 b using uniform priors, there is a strong preference for higher eccentricities (≥ 0.9; see Figure 14). With observable-based priors, however, the eccentricity distribution strongly shifts towards lower values. We validate the lower-eccentricity solution in two ways: (1) by checking consistency with an updated orbit fit that incorporates the two new RV measurements reported in Section 2.2, and (2) by assessing the stability of the system, given its age. As discussed below, both approaches validate the lower-eccentricity solution, indicating that the observable-based prior, in this case, is doing exactly what it was designed to do: minimize the biases in the posteriors when in the priordominated regime.
1. We fit HD 1160 b's orbit using uniform priors with orbitize!. The main difference from the previous fit is that we include two radial velociy (RV) data points, presented in Section 2.2. The addition of two RV data points was sufficient to shift the eccentricity distribution, yielding the distribution presented in Figure 15. This distribution resembles the observable-based priors distribution (even without RVs) much more closely (see Figure 5).
2. To analyze whether the system would allow for such a high eccentricity for one of the companions, we ran stability simulations for the age of the system. Using the WHFast integrator on REBOUND (Rein & Tamayo 2015), we integrate the system for 120 Myr, which is the oldest estimate for the system's age given by Garcia et al. 2017. We assess stability using the Mean Exponential Growth of Nearby Orbits (MEGNO) (Cincotta et al. 2003). 1160 c from orbitize! to assess the stability of the system. We bin the b and c posteriors into 10 bins, with each bin spanning a 0.1 spacing in eccentricity (e.g., 0.1 < e b < 0.2, 0.3 < e c < 0.4). For each bin, we draw 1,000 random orbit combinations from the posteriors and assess whether they remain stable for the age of the system. This yields 100,000 possible orbit combinations from the objects' posteriors, with 1,000 per bin (and a total of 100 bins). We then assess how many orbits in each bin remain stable using the MEGNO parameter. We do this for two different posteriors for HD 1160: with and without the RVs.
We present a contour plot of the 100,000 possible orbit combinations of HD 1160 b and c, with the contours representing what percent of the random draws remained stable for the age of the system. The contour is shown in Figure 6. The figure illustrates that HD 1160 b is unlikely to have eccentricity values > 0.9, which was a value favored by the uniform prior fit without radial velocity measurements. When we test this simulation set but with HD 1160 b's orbital posteriors without RVs, we obtain an equivalent plot, showing that the majority of solutions for HD 1160 b and c disfavor very eccentric orbits for both objects.
For PZ Tel b, a new astrometry point was sufficient to also change the eccentricity distribution of the object. The 2018 epoch added to our fits is shown in Table  2. Prior to 2018, works that presented orbit fits for  . Stability contour for the HD 1160 system obtained by binning the companions' posteriors in eccentricity bins of 0.1 width. The colors represent percent of the 1,000 random draws for each bin that remain stable for the age of the system. HD 1160 b and c orbits are shown to disfavor very eccentric orbits by the stability constraint. this object indicated that it had a very high eccentricity (e.g., Ginski et al. 2014, Beust et al. 2016, Maire et al. 2016b). The new astrometry point for PZ Tel b, obtained in 2018, changes the orbital period coverage of the companion from 3.9 +1.9 −1.3 % to 6.2 +2.9 −2.1 %, calculated from the weighted median of the posterior distribution (50th percentile), with the uncertainties encompassing the 68th percentile of our orbit fits for the period and

Leave-One-Out Cross Validation
Given the small sample size of directly imaged companions, there is a possibility that a single object with a narrow eccentricity posterior may noticeably impact the resulting population distribution. To test for this, we analyze the impact of removing one single object from the sample for each of the 21 objects. We do that by re-running the beta distribution fit on the entire sample, but excluding one object at a time. We characterize a large impact as a change in expected α and β parameters for the underlying distribution of eccentricity that falls outside of the 1-σ confidence level for the parameter values shown in Section 4.1.1.
We find that removing any single object from our population does not yield δα and δβ to be outside of the 68% confidence interval using P defined in the previous Section. We therefore conclude that our sample does not have obvious outliers. The object which produces the most different distribution is 1RXS0342+1216 b, which is a high mass companion with eccentricity of 0.94 +0.04 −0.11 (given from our orbit fit's 68th percentile, with the central value being the weighted median of the posteriors). Even in this case, we estimate that the two distributions have parameters that are consistent with a P of 0.50.

Separating the Sample by Companion Mass
We also separate the sample in two populations based on the mass of the companion. The goal is to identify any possible differences in the eccentricity as a function of mass. Such a distinction could be correlated with the masses expected for gas giant planets (companions with mass ≤ 15 M J ) and brown dwarf companions (companions with mass > 15 M J ). If differences in eccentricity distributions are found between the populations, eccentricities can potentially be used to constrain formation mechanisms. The complexity here is the use of a mass boundary for distinguishing between populations. If the definition of a planet versus a brown dwarf is related to formation, it is unlikely that 15 M J is a meaningful boundary (e.g., Bodenheimer et al. 2013;Mordasini et al. 2012). Additionally, most of the masses assumed for these companions are from evolutionary models and rely on uncertain properties such as system age (e.g., Carson et al. 2013b;Hinkley et al. 2013). Even with these caveats, it is still beneficial to look for differences between more and less massive companions in the context of distinction between these populations.
When classifying these objects into the low-mass and high-mass categories, the intermediate mass objects at the boundary, namely β Pic b, HR 2562 b and κ And b, could have a strong impact on the final result. Thus, we re-run our "planet" vs. "brown dwarf" population placing them in either of the populations to analyze whether the eccentricity distribution of each population changes. The mass estimates are plotted against the eccentricity estimates (68th percentile) in Figure 8. The mass estimates and respective references are shown in Table  1.
Including all of the boundary objects in the planet sample makes it comprised of 9 companions. When none of them are in the sample, the planet population is comprised of only 6 objects. This low sample size can cause an outlier to significantly skew the eccentricity of the entire population.
When switching β Pic b from the planet to brown dwarf population, we find no significant change in our P values that measure the consistency between the two populations. This is not the case for κ And b and HR 2562 b, in particular the former. Table 6 represents our Consistency Parameter (defined as the P value) for different combinations of these objects' classifications. κ And b, when considered a brown dwarf, produces distributions of planets and brown dwarfs that are different from each other. We illustrate two of these extreme cases (highest P and lowest P) in Figure 9. The "extreme" cases yield, respectively, planet population parameters  Table 6 that a small portion of objects can shift whether the two populations have similar distributions or not. This is likely due to the small sample size both for the entire population but especially for the sub-samples separating high and low mass objects.
We also test the mass separation distributions completely omitting the three boundary objects, and find planet and brown dwarf populations that are similar to Figure 9, bottom panel. However, the low sample size of the population (only 6 planets, none of which have eccentricities above ≈ 0.7 in their posteriors' 68th percentile) is significantly affected by a single object with higher eccentricity, such that our population distribution derivation is prone to higher uncertainties and to being skewed by any outliers in the sample. Given these results and the uncertainties and caveats around masses, we conclude that our current sample and its components' individual eccentricity distributions do not provide enough constraints to affirm that different formation pathways are underway for sub-stellar companions above and below ∼15 M Jup .

Separating the Sample by Companion Separation
We also test splitting our companion sample by separation from the host star. The interest in performing this comparison comes from the possible distinction between core accretion and gravitational instability companions, since gravitational instability models tend to favor formation further away from the host star (e.g. Mayer et al. 2004), while core/pebble accretion companions are most efficiently formed closer to the host star (e.g. Dodson-Robinson et al. 2009). The 68th percentile eccentricity posteriors are plotted against the separation of the companions, in AU, in Figure 10.
We also test splitting the population into "Close" and "Wide" separation companions, with the cutoff of 30 AU. This cutoff allows the companion sample to be split ≈ in half between the 5-30 AU and 30-100 AU bins. It is also close to the 35 AU threshold set by Dodson-Robinson et al. 2009, beyond which the critical mass obtained by core accretion objects cannot occur and where gravitational instability becomes a more efficient companion formation mechanism. We find that the two populations are consistent with each other, with a P of 0.85. Their distributions are shown in Figure 11.

Minimal Orbital Coverage Simulations
Given the major shift in the eccentricity distribution of PZ Tel b given the addition of new data (Section Note-Impact of the classification of intermediate mass objects on the consistency between the planet and brown-dwarf populations. The consistency is measured as P measuring the similarity of the estimated parameters α and β that define the Beta distribution. A P >0.68 would mean that the two distributions are consistent, while <0.01 means that the two distribution are likely different. For this table, β Pic b is considered a planet, although considering it a brown dwarf produces no significant changes in P. 4.1.2), we simulate how much orbital period coverage we need in order to obtain a meaningful posterior distribution for eccentricity. The main goal of this analysis is to stipulate whether we can obtain meaningful results on these objects with the currently available astrometry. For context given the sample we have been exploring, we examine the best estimate for the orbital coverage as a fraction of the period (in %), plotted against the best estimate for the eccentricity of each source, with error bars encompassing the possible values presented in each individual eccentricity distribution. We obtain the orbital coverage by dividing the current astrometric coverage for each object by its inferred orbital period. The results are presented in Figure 12. The average orbital period coverage for the sources in our sample is 7.4%, calculated from weighted median values of our sample's orbital period fits and the astrometric coverage, in years, for each object in the sample. In order to assess the minimum required orbital phase coverage, we simulate data to indicate how an object's orbital configuration affects the amount of coverage needed to obtain reliable posteriors. We perform this simulation in hopes of finding a trend for how much information must be obtained for each system such that a meaningful posterior can be extracted -and, in turn, a meaningful eccentricity underlying parent distribution can be obtained. First, we simulated orbits with properties that are representative of the average object in our sample. We define a period of 200 years, and inclination of 60 degrees and periastron passage occurring in 2150. The total system mass is fixed at 0.68 M . We vary the eccentricity of a given orbit to see how its value affects the posteriors for a given orbital period coverage. We calculate astrometry points based on these orbital parameters, starting in the year 2021. We generate astrometry points in order to simulate increasing orbital phase coverage. The astrometry "time step" represents 1% of the period of the orbit. The astrometry encompasses uncertainties on the order of milliarcsecond precision, as is typical for high contrast imaging instruments (Konopacky et al. 2014). To simulate instrument differences and other systematic factors, the simulated astrometry points also have noise randomly sampled from a Gaussian distribution centered at 0 and with the width of the astrometric uncertainty. This simulated astrometry is used to run orbit fits analogous to those in the rest of our analysis. We run our observable-based prior and uniform prior orbit fitter 100 times for each orbital eccentricity value, increasing the astrometry and hence the phase coverage with each successive run. In order to consider a simulation "successful", the real input eccentricity, To, inclination and period must be within the 68% confidence interval given by the fit.
We find that the orbital coverage of 15% was the minimum value needed to obtain 68 successful posteriors (out of 100 trials) for all eccentricity value variations, both for observable and uniform priors. This encompasses 30 years of observations from 2021 to 2051 for this specific orbit. The orbital coverage value is defined differently from the Lucy 2014 40% minimum orbital Figure 10. The semi-major axis plotted against the eccentricity of the companions. Semi-major axis values are calculated from the 68th percentile of our orbital fits. Error bars on the eccentricity represent the 68th percentile of the orbit fit. The blue dots represent "planets", or objects under 15 MJ , while red dots represent "brown dwarfs", or objects with masses above 15 MJ . Green dots are "boundary" objects: intermediate mass objects that could be in either distribution. In orange, we plot the separation cutoff where we split our sample between "Close" and "Wide" separations. Figure 11. The "Close" vs "Wide" separation companion populations, with the split happening at 30 AU. The two distributions are consistent with each other, with a P of 0.85. Figure 12. The inferred orbital coverage of each object in the sample, calculated as a fraction of the astrometric coverage and the period fit, plotted against the eccentricity of the object. This is from the real data used in our sample. Most of the sources have less than 10% orbital period coverage, with 4 objects presenting less than 1% period coverage.
phase coverage, since here it is based off of the period of the orbit and not the orbital arc. Thus the 15% period coverage is dependent on an orbital arc where the planet does not go through periastron passage during observation, which should be the case for the majority of directly imaged planets, as the planet spends the least amount of time of its orbit at periastron. The phase coverage obtained by this specific orbit spans ≈ 3 -15% of the orbit's phase (depending on eccentricity, which was varied in the interval from 0.1 to 0.9 in steps of 0.1), which is less than the value obtained by Lucy 2014. Another difference is that in Lucy 2014 the total system mass was a free parameter in the fit, while for our case (as is standard for most direct imaging orbit fits) the mass was kept fixed. Having the total system mass as a free parameter is likely a reason as to why Lucy 2014's phase coverage requirement is larger than the ones obtained with our simulations. We also test the needed orbital coverage for 100-year and 400-year period orbits and find no dependence on period as long as the orbital period coverage remains the same (at least 15% of period). Note that this is a rule of thumb only -some orbits may require more coverage than this, while others may require less. The 15% coverage number found here represents a minimum point at which the orbit posteriors for directly imaged planets should realistically be used to derive population distributions. Our results are presented in Figure 13.
For the median period of objects in our sample (≈ 251 years), we would need 38 years of data to begin to have reliable posteriors. Looking at specific example cases, HD 19467 B has an estimated period of 426.14 +115.89 −98.13 years (weighted median of the fit, with uncertainties spanning the 68th percentile of orbital fit) but our astrometry only covers 7 years (or 1.6 +0.5 −0.4 %), calculated from the 68th percentile of the fit, of the estimated orbital period -which means that our parameter posteriors for this object may not be reliable. This result is also exemplified in Section 4.1.2, where PZ Tel b's orbital coverage increase from 3.9 +1.9 −1.3 completely shifted the eccentricity distribution of the companion. However, given that we still have less than 15% coverage for PZ Tel b, we might expect additional changes in the posteriors as more data is collected.

DISCUSSION
In this work, we use observable-based priors to fit for the orbital parameters of 21 directly imaged substellar companions with updated astrometry and radial velocity measurements from the literature. Revisiting the analysis done by Bowler et al. 2020 with a different set of priors, the goal is to obtain the eccentricity distribution of the entire population and analyze if there is any difference in the eccentricity distributions of planet and brown dwarf populations. Differences in these populations' eccentricity distributions could potentially be signs of different formation mechanisms: with planets mainly forming from protoplanetary disks and thus presenting low eccentricities, and brown dwarfs forming similarly to binary star systems (for instance via disk fragmentation) and therefore spanning all values of eccentricities .

Choice of Priors
The vast majority of exoplanets and brown dwarfs in our directly imaged population have long periods, due to their fairly large distance from the host star. Consequently, their orbital arcs covered by current astrometry span a small fraction of their orbital periods, yielding orbital fits that are severely undersampled. Such undersampling of data can cause all results of orbit fitting to be in a prior-dominated regime. Choosing uniform priors for all orbital parameters results in equal weighting between all possible phases of the orbit. This is demonstrably not the case from basic orbit physics, which dictates that the planet presents the highest orbital speed during periastron passage (T o ). Thus the probability of catching these long period companions right at T o is low, since it spends the least amount of time there throughout its orbit. In the case of well-sampled data covering Figure 13. The bar plot of how many successful tries (i.e., true parameters are within 1 σ of posterior distribution) we obtain as a function of eccentricity for observable and uniform priors. For 15% of the orbital period in coverage, all of the eccentricity values from 0.1 to 0.9 have over 68 of the 100 trials considered successful. a large fraction of the orbit, the information contained in the data itself is sufficient to overcome the choice of prior. However, in the case of minimally sampled data, which corresponds typically to linear motion of the planet without measureable higher order motion such as acceleration, landing in the prior-dominated regime can result in unintended effects. One such example is that astrometric systematics could be interpreted by the fitter as rapid acceleration, hence leading to a best fit that is near periastron. This known bias caused by using uniform priors warranted a different approach to priors when it came to orbit fitting. Observable-based priors (O'Neil et al. 2019) provided such an approach, where the uniformity was in the orbital observables rather than in the orbital parameters. These priors decrease the bias towards high eccentricities and periastron passage during observations, effectively by down-weighting the likelihood of an orbit fit that has T o near the time of observations. While it is the case that proper approaches to prior definition in Bayesian statistics is a matter of debate, it is important to recognize that in the priordominated regime, as we are here, prior choices have measurable consequences in the posteriors and mitigation of biases is desirable.
Differences in eccentricity posteriors of some of the objects in our sample are indeed noticed when it comes to comparing these two different priors -in particular for HR 2562 b, HD 1160 b, HD 49197 b, HIP 65426 b and PZ Tel b (for corner plots of these objects in comparison to uniform priors, see Appendix A). In some cases, the long tail of high eccentricities completely disappears (e.g. HD 1160 b) when using observable priors. This is particularly important given our simulation showing that the system only remains stable if HD 1160 b has an eccentricity of < 0.9. Indeed, when radial velocity data is added to HD 1160 b's orbital fit with uniform priors, this is validated -we see a significant shift in its eccentricity posteriors to lower eccentricities, more consistent with what observable priors obtained with an astrometry-only solution. Stability arguments have been successfully used in other cases to help inform the orbital parameters of directly imaged planets, and represent a powerful means of constraint beyond astrometry alone (e.g., Wang et al. 2018).

Population Results
Given the major shift in many of the posteriors using observable-based priors, we expected some change in the population distribution from the combined dataset. While the specific values of α and β are different from Bowler et al. 2020, the overall result is very similar -a nearly uniform distribution, but with a lower incidence of high eccentricity orbits. Since we did not include all planets from the same multi-planet system in the sample (but rather chose one as the system's "representative") due to the likely preference that multiple planet systems have for lower eccentricities due to stability requirements (Wright & Howard 2009), it is possible that inclusion of these objects (such as HR 8799 b,d,e) could have yielded a population eccentricity distribution with a stronger preference for lower eccentricities. Overall, since the distributions are statistically consistent with each other, we can conclude that current data yields an eccentricity distribution that is largely flat with eccentricity.

Planets and Brown Dwarfs
Given the fact that the sample size is small and some of the posteriors are prior-dominated, care must be taken when interpreting parameter-based subdivisions of the sample. This is exemplified when we separate the sample by mass into "planets" and "brown dwarfs" (i.e. low and high mass objects), as we have 3 objects (HR 2562 b, κ And b and β Pic b) that fall into a "boundary" or "intermediate" mass classification -objects whose mass estimates allow for placement in both classifications. In such cases, we tested placing them in either population to see if results would change. Indeed, placing κ And b and HR 2562 b into the planet population yields completely different results from placing them into the brown dwarf population. Given these results, we conclude the uncertainty in both model-derived masses and individual eccentricity distributions are too large to allow us to distinguish between different populations as a function of mass. This means that at this time, we cannot use eccentricities to say that different formation pathways are underway for sub-stellar companions above and below ∼15 M Jup . We note again that the distinction between populations at this mass boundary is not particularly meaningful from a formation perspective, but rather is a useful "breakpoint" to explore these possibilities.

Minimal Orbital Coverage
Since we are mostly in the prior-dominated regime and there were major changes in posteriors with only a small addition of data, it is important to have some handle on how much data is needed to begin to make meaningful population distributions for directly imaged planets. Our approach here, to look at when we have enough orbital coverage for the eccentricity distribution to cover the true value in simulated data, is one possible way to answer that question. Our finding of 15% coverage is valid for both observable-based and model-based, or uniform, priors. The results converge at that point, suggesting that ultimately the distinction between these two prior sets really only matters for interpreting results in the prior-dominated regime. However, this is not the only possible set of simulations that should be conducted to get a firm handle on the necessary orbital phase coverage. For instance, we did not consider the role of some of the orbital parameters, such as inclination, in defining the needed coverage. Edge-on orbits are often more quickly constrained than face-on, suggesting that this parameter is important to consider. Furthermore, we did not explore in depth the role that true time from periastron plays in the defining of meaningful posteriors. Instead, we take 15% as a useful guideline for a general orbit where the companion is not close to periastron, which should encompass the majority of real systems.
As noted above, the average percent coverage of our sample is 7.4%, with most objects having astrometry spanning less than 10% of the orbital period. Given these results, we conclude again that the undersampling of the orbital period for each individual object in our sample coupled with our small sample size does not allow us to affirm that planets and brown dwarfs have different formation pathways. Until the orbits of these objects are well sampled enough to provide meaningful posteriors, the eccentricity of the population of directly imaged substellar companions remains essentially unconstrained.

CONCLUSION
The main findings of this study are: • We derive new orbital parameter posteriors for a set of 21 directly imaged substellar companions using observable-based priors.
• Several companions have resulting eccentricity distributions that change significantly from previous results. The inclusion of radial velocity (RV) data points or new astrometry in the orbit fitting process shifts the resulting posteriors for both observable-based and uniform prior fits.
• We derive a population-level eccentricity distribution for the 21 companions and obtain shape parameters α = 1.09 +0.30 −0.22 and β = 1.42 +0.33 −0.25 . These values are consistent with Bowler et al. 2020's parameters obtained using uniform priors, but with a lower incidence of high eccentricity objects.
• Separating the population into "giant planets" and "brown dwarfs" produces different results depending on where intermediate mass objects are placed. This result implies that our current sample size and large uncertainties may not be sufficient to determine whether these objects do in fact present distinct eccentricity populations.
• From our orbital coverage simulations, we find that one generally needs 15% orbital period coverage to obtain a reliable posteriors on eccentricity, period and T o posterior .
Following the conclusions of this work, the addition of new astrometry or radial velocity points for directly imaged companions can help further constrain the eccentricity posteriors of directly imaged companions to reliable intervals. This can in turn allow for a more robust eccentricity distributions at a population level.

ACKNOWLEDGEMENTS
We thank Christopher Theissen and Sarah Blunt for helpful discussions on the analysis and results presented in this work. We also thank an anonymous referee for their comments which helped improve this manuscript.
Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The W. M. Keck Observatory was made possible by the financial support of the W. M. Keck Foundation. The authors wish to acknowledge the significant cultural role that the summit of Maunakea has always had within the indigenous Hawaiian community. The author(s) are extremely fortunate to conduct observations from this mountain. Portions of this work were conducted at the University of California, San Diego, which was built on the unceded territory of the Kumeyaay Nation, whose people continue to maintain their political sovereignty and cultural traditions as vital members of the San Diego community. C