# First VLTI/GRAVITY Observations of HIP 65426 b: Evidence for a Low or Moderate Orbital Eccentricity

Sarah Blunt, W. O. Balmer, J. J. Wang, S. Lacour, S. Petrus, G. Bourdarot, J. Kammerer, N. Pourré, E. Rickman, J. Shangguan, et al.

► To cite this version:

Sarah Blunt, W. O. Balmer, J. J. Wang, S. Lacour, S. Petrus, et al. First VLTI/GRAVITY Observations of HIP 65426 b: Evidence for a Low or Moderate Orbital Eccentricity. The Astronomical Journal, 2023, 166, 10.3847/1538-3881/ad06b7. insu-04479028

HAL Id: insu-04479028

https://insu.hal.science/insu-04479028

Submitted on 10 Mar 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# First VLTI/GRAVITY Observations of HIP 65426 b: Evidence for a Low or Moderate Orbital Eccentricity

Sarah Blunt<sup>1,2</sup>, W. O. Balmer<sup>3,4</sup>, J. J. Wang<sup>1</sup>, S. Lacour<sup>5,6</sup>, S. Petrus<sup>7,8</sup>, G. Bourdarot<sup>9</sup>, J. Kammerer<sup>4</sup>, N. Pourre<sup>10</sup>, E. Rickman<sup>11</sup>, J. Shangguan<sup>9</sup>, T. Winterhalder<sup>6</sup>, R. Abuter<sup>6</sup>, A. Amorim<sup>12,13</sup>, R. Asensio-Torres<sup>14</sup>, M. Benisty<sup>10</sup>, J.-P. Berger<sup>10</sup>, H. Beust<sup>10</sup>, A. Boccaletti<sup>5</sup>, A. Bohn<sup>15</sup>, M. Bonnefoy<sup>10</sup>, H. Bonnet<sup>6,10</sup>, W. Brandner<sup>14</sup>, F. Cantalloube<sup>16</sup>, P. Caselli<sup>9</sup>, B. Charnay<sup>5</sup>, G. Chauvin<sup>10</sup>, A. Chavez<sup>1</sup>, E. Choquet<sup>16</sup>, V. Christiaens<sup>17</sup>, Y. Clénet<sup>5</sup>, V. Coudé du Foresto<sup>5</sup>, A. Cridland<sup>15</sup>, R. Dembet<sup>5</sup>, A. Drescher<sup>9</sup>, G. Duvert<sup>10</sup>, A. Eckart<sup>18,19</sup>, F. Eisenhauer<sup>9</sup>, H. Feuchtgruber<sup>9</sup>, P. Garcia<sup>13,20</sup>, R. Garcia Lopez<sup>14,21</sup>, E. Gendron<sup>5</sup>, R. Genzel<sup>9</sup>, S. Gillessen<sup>9</sup>, J. H. Girard<sup>4</sup>, X. Haubois<sup>22</sup>, G. Heißel<sup>5,23</sup>, Th. Henning<sup>14</sup>, S. Hinkley<sup>24</sup>, S. Hippler<sup>14</sup>, M. Horrobin<sup>18</sup>, M. Houlle<sup>25</sup>, Z. Hubert<sup>10</sup>, L. Jocou<sup>10</sup>, M. Keppler<sup>14</sup>, P. Kervella<sup>26</sup>, L. Kreidberg<sup>14</sup>, A.-M. Lagrange<sup>5,10</sup>, V. Lapeyrère<sup>5</sup>, J.-B. Le Bouquin<sup>10</sup>, P. Léna<sup>5</sup>, D. Lutz<sup>9</sup>, A.-L. Maire<sup>10</sup>, F. Mang<sup>9</sup>, G.-D. Marleau<sup>14,27,28,29</sup>, A. Mérand<sup>6</sup>, P. Mollière<sup>14</sup>, J. D. Monnier<sup>30</sup>, C. Mordasini<sup>29</sup>, D. Mouillet<sup>10</sup>, E. Nasedkin<sup>14</sup>, M. Nowak<sup>31</sup>, T. Ott<sup>9</sup>, G. P. P. L. Otten<sup>32</sup>, C. Paladini<sup>22</sup>, T. Paumard<sup>5</sup>, K. Perraut<sup>10</sup>, G. Perrin<sup>5</sup>, O. Pfuhl<sup>6</sup>, L. Pueyo<sup>4</sup>, J. Rameau<sup>10</sup>, L. Rodet<sup>33</sup>, Z. Rustamkulov<sup>34</sup>, T. Shimizu<sup>9</sup>, D. Sing<sup>3,34</sup>, T. Stolker<sup>15</sup>, C. Straubmeier<sup>18</sup>, E. Sturm<sup>9</sup>, L. J. Tacconi<sup>9</sup>, E. F. van Dishoeck<sup>9,15</sup>, A. Vigan<sup>16</sup>, F. Vincent<sup>5</sup>, K. Ward-Duong<sup>35</sup>, F. Widmann<sup>9</sup>, E. Wieprecht<sup>9</sup>, E. Wiezorrek<sup>9</sup>, J. Woillez<sup>6</sup>, S. Yazici<sup>9</sup>, and A. Young<sup>9</sup>

## The exoGRAVITY collaboration

<sup>1</sup> Center for Interdisciplinary Exploration and Research in Astrophysics (CIERA) and Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208, USA

<sup>2</sup> Department of Astronomy, California Institute of Technology, Pasadena, CA 91125, USA

<sup>3</sup> Department of Physics & Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA

<sup>4</sup> Space Telescope Science Institute, Baltimore, MD 21218, USA

<sup>5</sup> LESIA, Observatoire de Paris, PSL, CNRS, Sorbonne Université, Université de Paris, 5 Place Janssen, F-92195 Meudon, France

<sup>6</sup> European Southern Observatory, Karl-Schwarzschild-Straße 2, D-85748 Garching, Germany

<sup>7</sup> Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de Valparaíso, Av. Gran Bretaña 1111, Valparaíso, Chile

<sup>8</sup> Núcleo Milenio Formación Planetaria—NPF, Universidad de Valparaíso, Av. Gran Bretaña 1111, Valparaíso, Chile

<sup>9</sup> Max Planck Institute for Extraterrestrial Physics, Giessenbachstraße 1, D-85748 Garching, Germany

<sup>10</sup> Université Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, France

<sup>11</sup> European Space Agency (ESA), ESA Office, Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA

<sup>12</sup> Universidade de Lisboa—Faculdade de Ciências, Campo Grande, 1749-016 Lisboa, Portugal

<sup>13</sup> CENTRA—Centro de Astrofísica e Gravitação, IST, Universidade de Lisboa, 1049-001 Lisboa, Portugal

<sup>14</sup> Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

<sup>15</sup> Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands

<sup>16</sup> Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France

<sup>17</sup> School of Physics and Astronomy, Monash University, Clayton, VIC 3800, Melbourne, Australia

<sup>18</sup> Institute of Physics, University of Cologne, Zülpicher Straße 77, D-50937 Cologne, Germany

<sup>19</sup> Max Planck Institute for Radio Astronomy, Auf dem Hügel 69, D-53121 Bonn, Germany

<sup>20</sup> Universidade do Porto, Faculdade de Engenharia, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal

<sup>21</sup> School of Physics, University College Dublin, Belfield, Dublin 4, Ireland

<sup>22</sup> European Southern Observatory, Casilla 19001, Santiago 19, Chile

<sup>23</sup> Advanced Concepts Team, European Space Agency, TEC-SF, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands

<sup>24</sup> University of Exeter, Physics Building, Stocker Road, Exeter EX4 4QL, UK

<sup>25</sup> Université Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange, Nice, France

<sup>26</sup> LESIA, Observatoire de Paris, Université PSL, Place Jules Janssen, 92195F-92195 Meudon, France

<sup>27</sup> Fakultät für Physik, Universität Duisburg-Essen, Lotharstraße 1, D-47057 Duisburg, Germany

<sup>28</sup> Institut für Astronomie und Astrophysik, Universität Tübingen, Auf der Morgenstelle 10, D-72076 Tübingen, Germany

<sup>29</sup> Physikalisches Institut, Universität Bern, Gesellschaftsstr. 6, 3012 Bern, Switzerland

<sup>30</sup> Astronomy Department, University of Michigan, Ann Arbor, MI 48109, USA

<sup>31</sup> Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK

<sup>32</sup> Academia Sinica, Institute of Astronomy and Astrophysics, 11F Astronomy-Mathematics Building, NTU/AS campus, No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan

<sup>33</sup> Center for Astrophysics and Planetary Science, Department of Astronomy, Cornell University, Ithaca, NY 14853, USA

<sup>34</sup> Department of Earth & Planetary Sciences, Johns Hopkins University, Baltimore, MD, USA

<sup>35</sup> Department of Astronomy, Smith College, Northampton, MA 01063, USA

Received 2023 June 2; revised 2023 September 27; accepted 2023 October 4; published 2023 November 21

## Abstract

Giant exoplanets have been directly imaged over orders of magnitude of orbital separations, prompting theoretical and observational investigations of their formation pathways. In this paper, we present new VLTI/GRAVITY astrometric data of HIP 65426 b, a cold, giant exoplanet which is a particular challenge for most formation theories

<figure><img src="image_2.png" /></figure>
at a projected separation of 92 au from its primary. Leveraging GRAVITY's astrometric precision, we present an updated eccentricity posterior that disfavors large eccentricities. The eccentricity posterior is still prior dependent, and we extensively interpret and discuss the limits of the posterior constraints presented here. We also perform updated spectral comparisons with self-consistent forward-modeled spectra, finding a best-fit ExoREM model with solar metallicity and C/O = 0.6. An important caveat is that it is difficult to estimate robust errors on these values, which are subject to interpolation errors as well as potentially missing model physics. Taken together, the orbital and atmospheric constraints paint a preliminary picture of formation inconsistent with scattering after disk dispersal. Further work is needed to validate this interpretation. Analysis code used to perform this work is available on GitHub: https://github.com/sblunt/hip65426.

*Unified Astronomy Thesaurus concepts: Exoplanets (498); Exoplanet formation (492); Astrometry (80); Optical interferometry (1168)*

*Supporting material:* data behind figure

# 1. Introduction

The existence of cold Jupiters (CJs) at large separations (10s of au) is a formation puzzle. Direct imaging surveys have revealed that these planets are intrinsically rare, with an occurrence rate of about 1% (Bowler & Nielsen 2018; Vigan et al. 2021), but given that the nominal core formation timescale at typical CJ separations is much longer than the disk dispersal timescale (Armitage 2020), we would not expect CJs to exist at all, assuming in situ core accretion. Direct gravitational collapse within a disk is an alternative to core accretion, but with its own issues, particularly very fast migration after formation (Nayakshin 2017; Vorobyov & Elbakyan 2018). The CJ population is quite diverse, with masses spanning an order of magnitude and separations spanning several (Bowler 2016), so it is also possible that multiple formation mechanisms are at play.

The eccentricities of CJs are a useful tracer of formation history, both at the population level (Bowler et al. 2020), and for individual systems. Marleau et al. (2019) conducted a detailed investigation of possible formation scenarios via core accretion for the ~14 Myr CJ HIP 65426 b (Table 1), following on the work of Coleman & Nelson (2016a, 2016b). At a projected separation of 92 au (Chauvin et al. 2017), this object is significantly farther from its star than, for example, the famous HR 8799 planets (at 16, 26, 41, and 72 au; Wang et al. 2018), and therefore even more of a challenge for formation via core accretion, as core formation timescale increases with separation in the widely separated regime. Marleau et al. (2019) derived distributions over the initial mass and luminosity of HIP 65426 b under several assumptions of postformation entropy, then coupled these initial conditions with N-body simulations to investigate which models could match the present-day conditions of the planet. They varied the initial number of planets in the system, and included prescriptions for type I and II migration in a protoplanetary disk. Through a suite of such simulations, they found two families of explanations for HIP 65426 b's current separation and luminosity. The first is core formation at close separations, followed by outward scattering and subsequent runaway gas accretion at the present-day location. The slower timescale of type II migration and subsequent disk dispersal allows the planet to remain in place after scattering and damps the postscattering eccentricity. The second scenario is similar, except runaway accretion and subsequent disk dispersal occur before scattering. Under this scenario, most simulations resulted in a high (>0.5) present-day eccentricity for HIP 65426 b. A

third scenario the authors mentioned but did not investigate in detail is the prospect of in situ formation, invoking a more rapid core formation process than typically assumed, such as pebble accretion (Lambrechts & Johansen 2014; Rosenthal & Murray-Clay 2018). This scenario would presumably result in a circular orbit.

In summary, Marleau et al. (2019) propose three statistically distinct formation pathways via core accretion for HIP 65426 b: (1) in situ formation, resulting in a circular orbit, (2) scattering before disk dispersal, resulting in a low-to-moderate eccentricity orbit, and (3) scattering after disk dispersal, resulting in a high-eccentricity orbit. Both scenarios 2 and 3 are often accompanied by an inner giant planet. An important caveat is that these are not hard-and-fast distinctions; an eccentricity of 0, even with no model uncertainty, would not unequivocally rule out two of the three scenarios. However, such a measurement will allow us to assign statistical probabilities to each formation scenario by comparing with the population synthesis outputs of studies like Marleau et al. (2019). Precise measurements of eccentricity will also better constrain the population-level eccentricity distribution of CJs, which is currently still prior dependent (Nagpal et al. 2023), allowing us to model the formation mechanisms responsible for the population as a whole.

HIP 65426 b has been astrometrically monitored since its discovery by Chauvin et al. (2017), and previous papers report an essentially unconstrained eccentricity posterior that reproduces the prior (Chauvin et al. 2017; Cheetham et al. 2019; Carter et al. 2023). In this paper, we update the orbit model of HIP 65426 b using new astrometry from the optical Very Large Telescope Interferometer (VLTI)/GRAVITY (~50× more precise than previous astrometry) as part of the ExoGRAVITY project (Lacour et al. 2020).

Spectral and photometric measurements of HIP 65426 b have also been used to infer the planet's atmospheric properties. Petrus et al. (2021) recovered a bimodal posterior by comparing HIP 65426 b measurements with BT-Settl CIF-IST models, with one small radius (~1 $R_J$) peak and one larger radius peak (~1.2 $R_J$). Their Exo-REM comparisons yielded broader posteriors encompassing both of these possibilities, as Exo-REM predictions were only available for the K band at the time, and they could not compare with all available data. These yielded $T_{\text{eff}} = 1560 \pm 100$ K, $\log g < 4.40$, a slightly super-solar metallicity of $0.05_{-0.22}^{+0.24}$, and an upper limit on the C/O ratio (≤0.55). Carter et al. (2023) subsequently published seven long-wavelength photometric data points of HIP 65426 b using the NIRCAM and MIRI instruments on JWST, and used them to update the empirical bolometric luminosity of the planet to
<table><caption>**Table 1**
Relevant Physical Properties of HIP 65426 A and b</caption><thead><tr><th>Quantity</th><th>Value</th><th>Source</th></tr></thead><tbody><tr><td>Stellar spectral type</td><td>A2V</td><td>Carter et al. (2023)</td></tr><tr><td>Stellar mass</td><td>1.96 ± 0.04 M<sub>☉</sub></td><td>Chauvin et al. (2017)</td></tr><tr><td>Moving group membership</td><td>Lower Centaurus-Crux, 14 ± 4 Myr</td><td>Gagné et al. (2018)</td></tr><tr><td>Stellar parallax</td><td>9.30 ± 0.03 mas</td><td>Gaia DR3 (Gaia Collaboration et al. 2016, 2022)</td></tr><tr><td>Hot-start planetary mass estimate</td><td>7.1 ± 1.2 M<sub>J</sub></td><td>Carter et al. (2023)</td></tr><tr><td>Hot-start planetary radius estimate</td><td>1.44 ± 0.03 R<sub>J</sub></td><td>Carter et al. (2023)</td></tr><tr><td>Hot-start planetary T<sub>eff</sub></td><td>1283<sup>+25</sup><sub>-31</sub> K</td><td>Carter et al. (2023)</td></tr></tbody></table>

<table><caption>**Table 2**
Observing Log</caption><thead><tr><th rowspan="2">Date</th><th colspan="2">UT Time</th><th rowspan="2">Target</th><th rowspan="2">NEXP/NDIT/DIT</th><th rowspan="2">Airmass</th><th rowspan="2">τ<sub>0</sub></th><th rowspan="2">Seeing</th><th rowspan="2">Fiber Pointing<br>ΔR.A./ΔDecl.</th></tr><tr><th>Start</th><th>End</th></tr></thead><tbody><tr><td>2021-01-07</td><td>06:22:05</td><td>06:48:54</td><td>HD 91881</td><td>8/64/1 s</td><td>1.05–1.10</td><td>6.6–8.0 ms</td><td>0"51–0"65</td><td>−1108/710 mas</td></tr><tr><td>2021-01-07</td><td>06:59:10</td><td>07:57:00</td><td>HIP 65426 b</td><td>4/8/100 s</td><td>1.38–1.64</td><td>5.0–6.6 ms</td><td>0"61–1"19</td><td>416/−704 mas</td></tr><tr><td>2021-01-07</td><td>08:03:17</td><td>08:06:18</td><td>HIP 65426 A</td><td>2/64/1 s</td><td>1.35–1.36</td><td>4.6–5.4 ms</td><td>1"05–1"23</td><td>0/0 mas</td></tr><tr><td>2022-01-23</td><td>05:16:56</td><td>05:30:01</td><td>HD 73900</td><td>6/64/1 s</td><td>1.05–1.10</td><td>13.2–20.4 ms</td><td>0"41–0"57</td><td>−825/−455 mas</td></tr><tr><td>2022-01-23</td><td>05:42:05</td><td>05:55:37</td><td>HD 91881</td><td>6/64/1 s</td><td>1.38–1.64</td><td>6.4–13.8 ms</td><td>0"30–0"67</td><td>−1108/710 mas</td></tr><tr><td>2022-01-23</td><td>06:08:31</td><td>06:51:52</td><td>HIP 65426 b</td><td>6/4/100 s</td><td>1.35–1.36</td><td>5.4–13.9 ms</td><td>0"54–0"72</td><td>418/−699 mas</td></tr><tr><td>2023-05-07</td><td>04:18:13</td><td>05:30:00</td><td>HIP 65426 b</td><td>6/4/100 s</td><td>1.16–1.29</td><td>1.6–2.0 ms</td><td>1"34–1"90</td><td>419/−696 mas</td></tr><tr><td>2023-05-07</td><td>05:47:00</td><td>05:55:52</td><td>HD 123227</td><td>4/96/0.3 s</td><td>1.20–1.22</td><td>1.8–2.1 ms</td><td>1"26–1"48</td><td>−404/ 889 mas</td></tr></tbody></table>

Note. NEXP, NDIT, and DIT denote the number of exposures, the number of detector integrations per exposure, and the detector integration time, respectively. τ₀ is the atmospheric coherence time during each exposure. The fiber pointing is the placement of the science fiber relative to the fringe tracking fiber (which is placed on the central star). HD 91881 and HD 73900 are two binary systems used for phase referencing.

$log L/L_{⊙} = (-4.31, -4.14)$. Together with hot-start cooling models, they used this luminosity estimate to infer the planetary mass, radius, and effective temperature (given in Table 1). They found that an independent comparison to the BT-Settl CIFIST grid yielded a good fit, but resulted in an unphysically small planetary radius of $1.06 \pm 0.05 R_J$ (consistent with the smaller radius mode of Petrus et al. 2021).

This paper is structured as follows: in Section 2, we present our new GRAVITY data, including three new astrometric epochs and a new K-band spectrum. In Section 3, we present our updated orbital solution and discuss the significance of our eccentricity measurement in detail. In Section 4, we compare all existing photometric and spectral measurements of HIP 65426 b with self-consistent model spectra to update our understanding of the planet's atmospheric properties. We find results that are consistent with but more precise than previous work, keeping in mind that systematic uncertainties in these fits are likely underestimated. Finally, in Section 5 we discuss the implications of our orbital and atmospheric inferences and call for additional observational and theoretical work.

## 2. Data

### 2.1. GRAVITY Data

We observed HIP 65426 b on 2021 January 7, 2022 January 23, and 2023 May 7 as part of the ExoGRAVITY Large Program (ESO Program ID 1104.C-0651, Lacour et al. 2020). We used the European Southern Observatory (ESO) VLTI's four 8.2 m Unit Telescopes (UTs) and the GRAVITY instrument (Gravity Collaboration et al. 2017). The observations primarily used the “off-axis dual field” mode, in which a

roof mirror is used to split the telescope field into two. The starlight goes to the fringe tracker to correct for atmospheric perturbations (Lacour et al. 2019). The exoplanetary light goes to the spectrograph configured with the medium-resolution grism (R = 500). Phase referencing of the metrology is obtained by swapping on a binary just before the observation. On 2021 January 7 we used the binary system HD 91881, on 2022 January 23 we used two binary systems, HD 73900 and HD 91881, and on 2023 May 07 we used HD 123227. During the night of the 2021 January 7, we also performed “on-axis single-field” observations of the host star to calibrate the spectrum of the planet for that night. The second and third nights of observations do not have an on-axis calibrator, and so we do not consider the spectrum of the planet from those observations. The observing log, presented in Table 2, records the length of the observations, the number of files recorded, and the atmospheric conditions. The placement of the science fiber was based on preliminary orbit predictions fit to the available relative astrometry at the time using the whereistheplanet<fnref n="36" /> software (Wang et al. 2021), and resulted in an efficient coupling of the planet flux into the science fiber.

For the 2021 January 7 epoch, we calculated the complex visibilites of the host and the companion using the Public Release 1.5.0 (2021 July 1³⁷) of the ESO GRAVITY pipeline (Lapeyrère et al. 2014). The observations were phase referenced with the metrology system using observations of the binary calibrator. We decontaminated the flux of the planet due to the host using a custom python pipeline (see Appendix

<footnote n="36">³⁶ http://www.wereistheplanet.com</footnote>

<footnote n="37">³⁷ https://www.eso.org/sci/software/pipelines/gravity/</footnote>
A of Gravity Collaboration et al. 2020), which treats the contamination as a polynomial dependent on time and baseline; a polynomial of fourth order was used for stellar light suppression.

We obtained the astrometric position of the planet relative to the star at each epoch by analyzing the phase of the ratio of coherent fluxes, computing a periodogram power map over the fiber's field of view (Figures 1 and 2). The mean astrometric position is taken to be the minimum $\chi^2$ value of this map. We estimated the uncertainty on each astrometric measurement using the scatter of the mean astrometric values between individual exposures. These new astrometric data points are provided in Table 3.

We then extracted the ratio of the coherent flux between the two sources at the location of the companion, generating a contrast spectrum from 2 to 2.5 µm. We only extract a spectrum for the observations from 2021 January 7, where we observed the star in on-axis mode directly after observing the planet in off-axis mode. We converted this contrast spectrum into a flux calibrated spectrum by multiplying it by a BT-NextGen (Allard et al. 2012) spectrum fit to archive photometry of the host. We used a BT-NextGen model with $T_{\text{eff}} = 8600$ K (Carter et al. 2023) and $\log(g) = 4.2$ (Bochanski et al. 2018), and scaled the model to fit the photometric measurements from Tycho2 (Høg et al. 2000), the Two Micron All Sky Survey (2MASS; Cutri et al. 2003), and the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010) using species (Stolker et al. 2020). The best-fit model is shown in Figure 3. We noted a tension between the absolute flux of the resultant spectrum and published SPHERE K-band photometry (Chauvin et al. 2017) and SINFONI spectrum (Petrus et al. 2021). For our 2021 January 7 observation, $F_{\text{K2}} = 1.40 \times 10^{-16} \text{ W m}^{-2} \text{ µm}^{-1}$ vs the SPHERE value $F_{\text{K2}} = 7.21 \times 10^{-17} \text{ W m}^{-2} \text{ µm}^{-1}$. This is likely due to two factors: the low integration time on the star (<3 minutes) and the degradation/variability of the observing conditions between the observation of the planet and the star. Since SPHERE and SINFONI have more robust measurements of the stellar flux during the course of their observations (from satellite spots for SPHERE, and direct measurements of the stellar flux for SINFONI), we opted to integrate the GRAVITY spectrum over the Paranal/IRDIS_D_K12_2 filter and scale the GRAVITY spectrum by ×0.52 to match the SPHERE K2 photometry. The existing SPHERE photometry matches the existing SINFONI spectrum (described more and compared with our new GRAVITY spectrum in Section 4) without any scaling factor applied.

The fluxes, uncertainties, and interchannel flux covariances are available as a machine-readable FITS file published as data behind the figure.

## 2.2. Literature Data

Literature spectral and photometric data used in this study (all of which are plotted in Figure 14) come from several sources. The VLT/SPHERE IFS spectra at Y-H (0.96-1.64 µm bands), four photometric VLT/SPHERE measurements at the H2 ($\lambda_0 = 1.593$ µm and $\Delta\lambda = 0.11$ µm) and H3 ($\lambda_0 = 1.667$ µm and $\Delta\lambda = 0.12$ µm), at K1 ($\lambda_0 = 2.1025$ µm and $\Delta\lambda = 0.204$ µm) and K2 ($\lambda_0 = 2.2550$ µm and $\Delta\lambda = 0.218$ µm) bands, and two NACO measurements at the Lp ($\lambda_0 = 3.80$ µm, $\Delta\lambda = 0.62$ µm) and Mp ($\lambda_0 = 4.78$ µm, $\Delta\lambda = 0.59$ µm) bands, which were originally published in Cheetham et al. (2019). An additional NACO point obtained

using the NB405 filter was originally published in Stolker et al. (2020). The medium-resolution SINFONI spectrum at K band comes from Petrus et al. (2021). All spectral and photometric data (with the exception of the SINFONI data, which were kindly provided by S. Petrus) were accessed from species<fnref n="38" /> (version 0.5.5).

We also supplemented our new GRAVITY astrometry points with astrometric measurements from the literature (Table 4). These data were also taken verbatim from the papers indicated above. We used all of the available published astrometry in our orbit fits, with the exception of the JWST astrometry published in Carter et al. (2023), and the NB405 NACO astrometry published in Stolker et al. (2020). These points were excluded because they have large error bars relative to the other literature astrometric points (i.e., they contribute little additional information), and there is only a single astrometric point for each of these instrument/filter combinations (i.e., it is hard to tell whether they are affected by systematic errors).

## 3. Orbit Analysis

In order to understand how our new GRAVITY data constrain the orbit of HIP 65426 b, we performed seven different orbit fits, varying the data subset used and the priors applied. These fits use literature data from Table 4³⁹ and new VLTI/GRAVITY data from Table 3. The results are summarized in Table 5 and described in more detail below. We used version 2.2.2 of **orbitize!** (Blunt et al. 2020; Blunt et al. 2023) for all fits, taking advantage of the ability to fit companion RVs introduced in version 2.0.0. **orbitize!** is a Bayesian tool for computing the posteriors over orbital elements of directly imaged exoplanets, and is intended to meet the orbit-fitting needs of the high-contrast imaging community. All orbit fits used the following orbital basis: semimajor axis (a), eccentricity (e), inclination (inc; with inc = 0 for a face-on orbit), argument of periastron ($\omega_p$), position angle of nodes (Ω), and epoch of periastron, defined as the fraction of the orbital period past a specified reference date (see Section 2.1 of Blunt et al. 2020 for a bit more detail). We use MJD = 58,859 as the reference date in this paper. Following the suggestion in Householder & Weiss (2022), we also explicitly specify here that we assume that $\omega_p$ is defined relative to the longitude of ascending node, which in turn is defined as the intersection point between the orbit track and the line of nodes when the object is moving away from the observer. Our coordinate system is illustrated in Figure 1 of Blunt et al. (2020). An interactive tutorial is also available to help users build an intuitive understanding of the coordinate system.<fnref n="40" /> We used the same priors on all parameters as given in Blunt et al. (2020), unless otherwise specified.

The orbit fits presented in this paper include two types of input data: relative astrometry and companion RVs. The ability to fit companion RVs is new in **orbitize!** since the publication of Blunt et al. (2020), so we discuss it in a bit of detail here. Following standard practice (e.g., Chapter 1 of Seager 2010), **orbitize!** by default defines orbital

<footnote n="38">³⁸ https://github.com/tomasstolker/species/blob/main/species/data/companion_data.json</footnote>

<footnote n="39">³⁹ As mentioned in the previous section, there are additional astrometric data points given in Stolker et al. (2020) and Carter et al. (2023), which were excluded because of their large error bars.</footnote>

<footnote n="40">⁴⁰ https://github.com/sblunt/orbitize/blob/main/docs/tutorials/show-me-the-orbit.ipynb</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Detections of HIP 65426 b with VLTI/GRAVITY in epochs 1 and 2. Both periodogram power maps visualize the χ² fit to the interferometric observables assuming a point source, after removing the contribution of the star using a fourth-order polynomial. The outer dashed gray circle indicates the effective fiber field of view (60 mas in diameter), and the red circles indicate the most probable planet position at each epoch. The planet is detected at high confidence in both epochs (periodogram power >500).</figcaption></figure>

parameters (a, e, etc.) as those of the relative orbit, meaning that astrometric and RV measurements of the secondary are both assumed to be relative to the primary.<fnref n="41" /> In order to use the RV measurement of Petrus et al. (2021) in our orbit fits, we therefore subtracted the absolute measurement of the planet’s RV (derived from the planetary spectrum) and an independent

measurement of the star’s RV derived from 78 HARPS spectra of the primary (see Section 4.2 of Petrus et al. 2021 for details), and propagated the uncertainties in both measurements.<fnref n="42" />

<footnote n="41">⁴¹ When RV measurements are available for both the star and the planet, *orbitize!* instead assumes that the planetary and stellar RV measurements are relative to the system barycenter.</footnote>

<footnote n="42">⁴² It is important to note that systematic offsets in the wavelength solutions used for the HARPS and SINFONI spectra (not to mention RV offsets due to astrophysical variability, such as pulsations) would affect the RV value derived here. However, the RV value derived from the SINFONI data is not sufficiently precise to impact our measurements of eccentricity and inclination, the main focus of this paper. Therefore we do not investigate any potential systematics in detail.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2. Detection of HIP 65426 b with VLTI/GRAVITY in epoch 3. See Figure 1.</figcaption></figure>

<table><caption>**Table 3**
New Relative Astrometry of HIP 65426 b Presented in this Paper</caption><thead><tr><th>Date<br>(JD - 2,400,000.5)</th><th>ΔR.A.<br>(mas)</th><th>σ<sub>ΔR.A.</sub><br>(mas)</th><th>ΔDecl.<br>(mas)</th><th>σ<sub>ΔDecl.</sub><br>(mas)</th><th>ρ<sub>Δ R.A., ΔDecl.</sub></th></tr></thead><tbody><tr><td>59,221.312</td><td>415.613</td><td>0.107</td><td>-708.133</td><td>0.073</td><td>-0.095</td></tr><tr><td>59,602.271</td><td>416.269</td><td>0.035</td><td>-705.051</td><td>0.035</td><td>-0.224</td></tr><tr><td>60,071.208</td><td>416.980</td><td>0.281</td><td>-701.373</td><td>0.248</td><td>0.954</td></tr></tbody></table>

**Note.** σ<sub>ΔR.A.</sub> and σ<sub>Δdecl.</sub> denote the uncertainties in astrometric position, and ρ<sub>ΔR.A.,Δdecl.</sub> denotes the correlation between the σ<sub>ΔR.A.</sub> and σ<sub>Δdecl.</sub> measurements.

Relative companion RV measurements like this do not allow us to measure a dynamical mass for the companion, but have the potential to reduce the posterior uncertainties of the relative orbit. However, the uncertainty on the available companion RV measurement is too large to constrain the orbital parameters meaningfully beyond the constraints from the relative astrometry (see Figure 4). It does, however, reduce the 180° degeneracy between Ω and ω, which usually shows up as double-peaked posteriors on those parameters for relative-astrometry-only orbits. A more precise planetary RV would uniquely orient the planet in 3D space (to within the orbital parameter uncertainties). Figure 4 shows the RV predictions for each orbit in the posterior of our accepted fit, together with the measurement from Petrus et al. (2021). The current uncertainties of the companion and stellar RVs of HIP 65426 A and b limit the ability of the relative RV to constrain the orbit fit.

We performed the following orbit fits (summarized in Table 5). All fits included the companion RV as described in the above paragraph:

1. Only including literature data from SPHERE and NACO (i.e., no GRAVITY data).

2. Literature data plus the first epoch of GRAVITY astrometry.

3. Literature data plus the second epoch of GRAVITY astrometry.

4. Literature data plus all GRAVITY astrometry.

5. The first two epochs of GRAVITY data alone (i.e., no literature data, and no third GRAVITY point).

6. Fit 4, except fixing eccentricity to zero.<fnref n="43" />

7. Fit 4, except applying a decreasing prior on eccentricity.

For the final fit, we used the following prior distribution, following Nielsen et al. (2008):<fnref n="44" />

$$p(e) = -2.18e + 2.01. \qquad (1)$$

Fits 1–5 were performed in order to assess the outlier sensitivity of our fits, as well as to understand the relative constraining power of the GRAVITY astrometry and the less precise literature measurements. The final two fits were performed to understand the prior dependence of the inferred eccentricity, as well as the impact of the eccentricity–inclination degeneracy (see, e.g., Ferrer-Chávez et al. 2021).

For all fits, we used the ptemcee implementation of emcee (Foreman-Mackey et al. 2013; Vousden et al. 2016), a parallel-tempered affine-invariant Markov Chain Monte Carlo (MCMC) sampling algorithm. All runs used 20 temperatures and 1000 walkers. After an initial burn-in period of 100,000 steps, each walker was run for 200,000 steps. Every hundredth step was saved.<fnref n="45" /> The chains were examined visually to determine whether they had converged, and the 200,000 steps of each chain post-burn-in were kept as the posterior estimate.

<footnote n="43">⁴³ Note that ω<sub>p</sub> and τ<sub>58849</sub> are undefined for a circular orbit.</footnote>

<footnote n="44">⁴⁴ The linear coefficient is from E. Nielsen (2023, private communication), and the constant value ensures the normalization. Note that this prior has zero probability beyond e = 0.92.</footnote>

<footnote n="45">⁴⁵ In orbitize!, this configuration corresponds to the following variable definitions: num temps=20; num walkers=1000; num -steps=200_000_000; burn steps=100_000; thin=100.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Best-fit spectral energy distribution (SED) model of HIP 65426 A used to convert our transmission functions of photometric bands (in the top panel, where each is plotted as a data point with errors immediately beneath). Middle: MAP BT-NextGen model and photometric data from TYCHO, 2MASS, and WISE. Bottom: MAP model residuals.</figcaption></figure>

<table><caption>**Table 4**

Astrometric Measurements from the Literature Used in the Orbit Fits Presented in this Paper</caption><thead><tr><th>Date<br>[JD − 2,400,000.5]</th><th>Separation (ρ)<br>(mas)</th><th>σ<sub>ρ</sub><br>(mas)</th><th>P.A.<br>(deg)</th><th>σ<sub>P.A.</sub><br>(deg)</th><th>RV<sub>pl</sub><br>(km s<sup>−1</sup>)</th><th>σ<sub>RVpl</sub><br>(km s<sup>−1</sup>)</th><th>Instrument</th><th>Reference</th></tr></thead><tbody><tr><td>57,538.4</td><td>830.4</td><td>4.9</td><td>150.28</td><td>0.22</td><td>...</td><td>...</td><td>SPHERE</td><td>(1)</td></tr><tr><td>57,565.5</td><td>830.1</td><td>3.2</td><td>150.14</td><td>0.17</td><td>...</td><td>...</td><td>SPHERE</td><td>(1)</td></tr><tr><td>57,791.0</td><td>827.6</td><td>1.5</td><td>150.11</td><td>0.15</td><td>...</td><td>...</td><td>SPHERE</td><td>(1)</td></tr><tr><td>57,793.1</td><td>828.8</td><td>1.5</td><td>150.05</td><td>0.16</td><td>...</td><td>...</td><td>SPHERE</td><td>(1)</td></tr><tr><td>57,891.0</td><td>832</td><td>3</td><td>149.52</td><td>0.19</td><td>...</td><td>...</td><td>NACO</td><td>(2)</td></tr><tr><td>57,892.0</td><td>850</td><td>20</td><td>148.5</td><td>1.6</td><td>...</td><td>...</td><td>NACO</td><td>(2)</td></tr><tr><td>58,250.0</td><td>822.9</td><td>2.0</td><td>149.85</td><td>0.15</td><td>...</td><td>...</td><td>SPHERE</td><td>(2)</td></tr><tr><td>58,250.0</td><td>826.4</td><td>2.4</td><td>149.89</td><td>0.16</td><td>...</td><td>...</td><td>SPHERE</td><td>(2)</td></tr><tr><td>58,263.5</td><td>...</td><td>...</td><td>...</td><td>...</td><td>14</td><td>15</td><td>SINFONI/HARPS</td><td>(3)</td></tr></tbody></table>

**Note.** P.A. is the position angle. RV<sub>pl</sub> indicates a measurement of the planet's radial velocity (RV) relative to the primary. Making this measurement involved separately measuring the absolute RV of the planet (from a SINFONI spectrum) and the absolute RV of the star (from a series of HARPS spectra), subtracting these quantities, and propagating the uncertainty. References: (1) Chauvin et al. (2017), (2) Cheetham et al. (2019), and (3) Petrus et al. (2021).

<table><caption>**Table 5**

Marginalized Posterior 68% Credible Intervals for the Free Parameters in Each of the Orbit Fits Described in Section 3</caption><thead><tr><th>Fit name</th><th>a (au)</th><th>e</th><th>inc (deg)</th><th>ω<sub>p</sub> (deg)</th><th>Ω (deg)</th><th>τ<sub>58849</sub></th></tr></thead><tbody><tr><td>Lit. data only</td><td>92.5<sup>+88.7</sup><sub>-32.0</sub></td><td>0.58<sup>+0.29</sup><sub>-0.38</sub></td><td>109.5<sup>+20.5</sup><sub>-11.2</sub></td><td>181.6<sup>+107.6</sup><sub>-128.9</sub></td><td>171.9<sup>+163.9</sup><sub>-45.3</sub></td><td>0.36<sup>+0.30</sup><sub>-0.17</sub></td></tr><tr><td>Lit. + GRAVITY epoch 1</td><td>66.1<sup>+28.3</sup><sub>-12.7</sub></td><td>0.64<sup>+0.26</sup><sub>-0.31</sub></td><td>105.6<sup>+13.9</sup><sub>-5.0</sub></td><td>195.6<sup>+74.0</sup><sub>-95.9</sub></td><td>156.9<sup>+171.0</sup><sub>-8.1</sub></td><td>0.44 ± 0.15</td></tr><tr><td>Lit. + GRAVITY epoch 2</td><td>66.3<sup>+27.2</sup><sub>-13.0</sub></td><td>0.63 ± 0.27</td><td>105.4<sup>+14.9</sup><sub>-4.8</sub></td><td>194.3<sup>+70.5</sup><sub>-101.6</sub></td><td>156.5<sup>+172.2</sup><sub>-7.6</sub></td><td>0.43 ± 0.15</td></tr><tr><td>All data</td><td>62.2<sup>+14.7</sup><sub>-1.2</sub></td><td>0.55 ± 0.19</td><td>109.2<sup>+6.3</sup><sub>-3.9</sub></td><td>189.3<sup>+33.3</sup><sub>-32.3</sub></td><td>155.1<sup>+9.0</sup><sub>-4.9</sub></td><td>0.44 ± 0.10</td></tr><tr><td>GRAVITY epochs 1 + 2 only</td><td>61.4<sup>+9.7</sup><sub>-5.1</sub></td><td>0.55 ± 0.14</td><td>109.6<sup>+4.4</sup><sub>-2.9</sub></td><td>190.3<sup>+25.8</sup><sub>-19.3</sub></td><td>155.4<sup>+7.2</sup><sub>-3.6</sub></td><td>0.44 ± 0.07</td></tr><tr><td>All data, fixed e = 0</td><td>91.2 ± 0.4</td><td>≡0</td><td>103.1 ± 0.5</td><td>...</td><td>153.2<sup>+17.9</sup><sub>-0.3</sub></td><td>...</td></tr><tr><td>All data, decreasing e prior</td><td>74.6<sup>+14.9</sup><sub>-10.0</sub></td><td>0.30 ± 0.17</td><td>105.2<sup>+2.5</sup><sub>-1.9</sub></td><td>187.6<sup>+50.5</sup><sub>-47.2</sub></td><td>153.8<sup>+5.9</sup><sub>-2.5</sub></td><td>0.46 ± 0.14</td></tr></tbody></table>

**Note.** The parallax and total mass were additional free parameters in every case. In all cases, the marginalized parallax 68% credible interval was 9.30 ± 0.03 mas and the marginalized total mass 68% credible interval was 1.96 ± 0.04 $M_{⊙}$.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Planetary RV predictions from the accepted fit (purple histogram), together with the planetary RV measurement from Petrus et al. (2021; dashed black line, with $1\sigma$ range in shaded gray). Takeaways: The planetary RV measurement does not constrain the orbital parameters, beyond breaking the 180° degeneracy for $\Omega$ and $\omega$. A planetary RV of $-3 \text{ km s}^{-1}$ is allowed, given the astrometry alone, but is disfavored (relative to a planetary RV of $3 \text{ km s}^{-1}$) because of the relative RV measurement.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Sky-projected visualization of the posterior of the orbit fit 4 described in the text. Left: 100 orbit tracks projected onto the plane of the sky, colored by elapsed time. The astrometric data are visible as pink points in the bottom left corner of the panel. Middle column: the same 100 posterior orbits (gray) in separation (top) and position angle (bottom) vs. time, together with the astrometric data used for orbit-fitting. Right column: the same 100 posterior orbits (gray), together with earlier (bottom) and later (top) astrometric measurements taken with VLTI/GRAVITY. $1\sigma$ and $2\sigma$ error ellipses are shaded in dark and light pink, respectively. Takeaway: the VLTI/GRAVITY epochs are $\sim 50 \times$ more precise than existing astrometric measurements of HIP 65426, and reduce the posterior uncertainty.</figcaption></figure>

We applied Gaussian priors on the parallax and total mass using the values given in Table 1. Both parameters do not significantly correlate with any other fit parameters, and the marginalized posteriors on these parameters reproduce their priors (visible, for example, in Figure 6).

As a final point, it is clear from scrutinizing Figure 5 that the two NACO points are discrepant from the contemporaneous SPHERE points. We briefly tested whether excluding the NACO points affected our fits by rerunning fit 6 without the two NACO points. Our results were indistinguishable, so we opted to keep the NACO points in all other fits.

### 3.1. The Impact of GRAVITY

Figures 5 and 6 visualize the posterior of fit 4, using all astrometric data and assuming a uniform prior on the eccentricity, and Figure 7 compares each of the orbit fits that use varying data subsets. These plots show the transformative impact of the GRAVITY data on the orbital uncertainties. The semimajor axis, eccentricity, and inclination posteriors all tighten significantly after including just three GRAVITY astrometric epochs. Figure 7 also shows that the results of a preference for moderate eccentricity and near edge-on inclination are robust to outliers, by showing nearly
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Corner plot of the posterior of the accepted orbit fit, using all data and assuming a uniform eccentricity prior. Diagonal panels show marginalized 1D histograms of posterior elements, and off-diagonal panels show 2D covariances between posterior elements. $1\sigma$, $2\sigma$, and $3\sigma$ contours are outlined in the covariance panels, and individual posterior samples outside of the $3\sigma$ boundaries are plotted directly as black dots. Takeaways: The 1D marginalized posterior distributions of the semimajor axis and inclination are well constrained. Strong covariances are apparent, in particular between the eccentricity, inclination, and semimajor axis.</figcaption></figure>

identical posterior distributions regardless of which GRAVITY epoch is used in the fit. It is also apparent from Figure 7 that the majority of the orbital parameter information is coming from the first two GRAVITY epochs, as evidenced by the similarity between the accepted fit and the GRAVITY-only fit. This test highlights the power of GRAVITY precision astrometry, but also warrants a warning: any unquantified systematics in the GRAVITY data could

significantly change the recovered orbital parameters. The orbital constraints reported here rely heavily on the detailed work that the exoGRAVITY team has undertaken to extract accurate astrometry (e.g., Lacour et al. 2019; Gravity Collaboration et al. 2020). HIP 65426 b should continue to be astrometrically monitored by GRAVITY in order to redistribute the impact of the first two (most precise) GRAVITY epochs reported here.
<figure><img src="image_9.png" /></figure>

**Figure 7.** Relative constraining power of the astrometric data for semimajor axis (top), eccentricity (middle), and inclination (bottom). The results of the following fits are shown and compared: (1) only literature astrometry (i.e., no GRAVITY data; gray), (2) literature astrometry and the first epoch of GRAVITY data (dark pink outline), (3) literature astrometry and the second epoch of GRAVITY data (light pink outline), (4) only GRAVITY astrometry (i.e., no literature data), and (5) all astrometric data (i.e., the accepted fit; purple outline). Takeaways: Most of the constraining power of the fit comes from the GRAVITY data, evident by the similarity between the GRAVITY-only fit and the accepted fit. In addition, neither GRAVITY point alone drives the fit, as evidenced by the similarity between fits (2) and (3). In other words, the posterior preference for moderate eccentricities is robust to the possibility that one of the two GRAVITY epochs is an outlier.

## 3.2. *A Moderate Eccentricity?*

All of the orbit fits that include GRAVITY data show a preference for moderate eccentricities ($e \sim 0.5$). In order to understand the significance of this preference, we performed a fit using a linearly decreasing eccentricity prior (fit 6; Equation (1)). Figure 8 shows that the eccentricity posterior is prior dependent, even though both a linearly decreasing and a uniform eccentricity prior result in preferences for nonzero eccentricities. We next performed a series of maximum-likelihood fits, fixing eccentricity to a specific value for each, in order to examine how the changing maximum-likelihood value was influencing the posterior shape. The results of this test are shown in Figure 9. This figure shows that the maximum likelihoods achieved by low ($e < 0.2$) and moderate ($e = 0.5$) orbits are comparable, but that the highest achievable maximum likelihood occurs at moderate eccentricities. This plot, along with Figures 6 and 8 allow us to construct the following explanation of the shape of the eccentricity posterior:

at low ($e < 0.2$) eccentricities, the likelihood is high but the prior volume is low (i.e., there are “fewer” circular orbits that fit the data well, even though those that do fit tend to fit very well). At moderate eccentricities, the prior volume is high and the likelihood is also high, leading to a posterior peak. At high eccentricities ($e > 0.6$), the prior volume is very high but the likelihood is low, leading to a decreasing posterior probability as a function of eccentricity. The “uptick” at very high eccentricities ($e > 0.8$) is caused by a degeneracy between eccentricity and inclination (see Ferrer-Chávez et al. 2021 for a detailed discussion), together with a nonlinear relationship between eccentricity and sky projection. As the eccentricity asymptotically approaches 1, the corresponding inclination must increase more and more rapidly to reproduce the astrometric data (see the covariance between eccentricity and inclination in Figure 6). This causes a characteristic “banana-shaped” covariance, typical of incomplete orbits (see Blunt et al. 2019 for another example in the context of an incomplete
<figure><img src="image_10.png" /><figcaption>**Figure 8.** 1D marginalized eccentricity posteriors for fits with uniform (purple) and linearly decreasing (pink) priors on the eccentricity. The priors themselves are also plotted as lines of the same colors. Takeaways: The eccentricity posterior depends on the choice of prior. However, both the linearly decreasing prior and the uniform prior result in posterior peaks at moderate eccentricity values.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Maximum log(likelihood) as a function of eccentricity. Although the MAP eccentricity is moderate (~0.5), the maximum likelihood at low and moderate eccentricities is comparable. This allows us to understand the shape of the eccentricity posterior (Figure 8): the likelihood is high at lower eccentricities, but the prior volume here is lower. The posterior “drop-off” at higher eccentricities is caused by a real decrease in likelihood. More eccentric orbits are less consistent with the data.</figcaption></figure>

RV orbit). We can explain the maximum a posteriori (MAP) shift to lower eccentricities when using a linearly decreasing prior, i.e., as occurring because the prior volume at moderate

and high eccentricities decreases when we switch to a linearly decreasing eccentricity prior. This interpretation is supported by both the Akaike information criterion (AIC) and Watanabe-Akaike information criterion (WAIC) model selection metrics, a point which is explored further in the Appendix.

## 4. Spectral Analysis

The GRAVITY K-band spectrum ($R = 500$) published in this work overlaps in wavelength coverage almost completely with the higher-resolution SINFONI spectrum ($R \sim 5600$), and therefore does not add significant spectral information to the HIP 65426 b SED. However, it is an independent constraint on the K-band spectrum. In this section we first compare the GRAVITY and SINFONI spectra, finding good agreement once the GRAVITY spectrum has been rescaled as described in Section 2, then globally compare the HIP 65426 b spectral data with self-consistent atmosphere models to update our understanding of this object's atmosphere.

In Figure 10, we plot the GRAVITY and SINFONI spectra. The agreement is excellent. Both spectra agree in magnitude (as expected after the rescaling) and shape.

Following Petrus et al. (2021), we derived atmospheric properties for HIP 65426 b by comparing with two sets of self-consistent atmospheric model grids: BT-SETTL CIFIST 2011c<fnref n="46" /> (Allard et al. 2003, 2007, 2011) and Exo-REM (Allard et al. 2012; Charnay et al. 2018), excellent summaries of which are given in Section 3.1 of Petrus et al. (2021). In the temperature range relevant for HIP 65426 b, the major advantages of Exo-REM include (1) the ability to explore nonsolar metallicities and (2) the success of the Exo-REM

<footnote n="46">⁴⁶ https://phoenix.ens-lyon.fr/Grids/BT-Settl/CIFIST2011c/</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** GRAVITY and SINFONI *K*-band spectra comparison. Top: GRAVITY (gray) and SINFONI (purple) $1\sigma$ flux confidence intervals are shown as filled bands. The SINFONI spectrum was resampled onto the GRAVITY wavelength grid using *spectres* (Carnall 2017). Bottom: the residuals, with propagated uncertainties, are shown relative to the flux = 0 line. Takeaway: the agreement between these two independent data sets is excellent. The GRAVITY fluxes, uncertainties, and interchannel flux covariances are provided in the FITS format as data behind the figure.</figcaption></figure>

(The data used to create this figure are available.)

cloud prescription at reproducing observations of dusty planets (like HIP 65426 b) near the L-T transition (Charnay et al. 2018). Ultimately, however, we are interested in comparing constraints from multiple independent models.

We used *species* (Stolker et al. 2020) to perform comparisons to both sets of models. In all cases, we computed posteriors using the *pyMultinest* (Buchner et al. 2014) Python interface to *multinest* (Feroz & Hobson 2008; Feroz et al. 2009, 2019) with 1000 live points. The results of all atmosphere model fits are shown in Table 6. We also opted not to include rotational broadening as a free parameter in all cases, following an experiment showing that doing so resulted in a broad uniform marginalized posterior over the planet's *vsini* and unchanged other fit parameters. The spectral resolution of the SIFNONI spectra translate to a minimum detectable *vsini* of ~50 km s⁻¹, which translates to a rotation period of ~0.1 day, assuming *i* = 90° and *R* = 1.5 *R*<sub>J</sub>, so the nondetection of rotational broadening is also consistent with our physical expectation.

## 4.1. *BT-SETTL CIFIST*

We performed four variations of *BT-SETTL CIFIST* comparisons to our full spectral data set, by (1) using either the GRAVITY or SINFONI spectrum<fnref n="47" /> in order to compare their relative constraining power and (2) fitting for correlated noise in the SPHERE IFS data, following Wang et al. (2020;

see their Equation (4)). All fits performed in this and the next section are summarized in Table 6. The results of the two fits allowing correlated noise in the SPHERE data are shown in Figure 11. The fit only including GRAVITY data is shown in purple, and the fit only including SINFONI data is shown in pink. The two fits are consistent overall, as we would expect given the consistency of the spectra themselves. Like Petrus et al. (2021; but unlike Carter et al. 2023), we also recover a bimodal posterior in surface gravity and effective temperature, regardless of which of the two *K*-band spectra we use.

Our initial *BT-SETTL CIFIST* fits allowed a free RV offset parameter for the SINFONI data, but this value consistently converged to unphysically large values, perhaps indicating an underlying problem with the model grid. We therefore opted to fix the RV of the SINFONI data to 0 for the fits presented here. Toggling on and off a GP for the SPHERE data does not significantly impact the physical parameters inferred, even though there are clear correlated residuals in the SPHERE data (Figure 12). The SPHERE residuals are visualized in Figure 12, where both modes of the posterior are plotted along with the SPHERE IFS data. Toggling on a GP for this data set allows the model to treat the deviations from both models as correlated noise, perhaps due to imperfect speckle subtraction at the planet location, photometric calibration errors due to telluric absorption, and/or model imperfections.

## 4.2. *Exo-REM*

Because the agreement between the *BT-SETTL CIFIST* model fits is good regardless of whether we use the SINFONI

<footnote n="47">⁴⁷ Although we downsampled the resolution of the SINFONI spectrum in order to compare with the GRAVITY spectrum in Figure 10, we fit the SINFONI spectrum at its native resolution.</footnote>
<table><caption>**Table 6**
68% Credible Intervals of Posterior Fits to Self-consistent Model Spectra Grids</caption><thead><tr><th>fit name</th><th>T<sub>eff</sub> (K)</th><th>log g</th><th>π (mas)</th><th>l<sub>SPHERE</sub> (μm)</th><th>A<sub>SPHERE</sub></th><th>(Fe/H)</th><th>C/O</th><th>R (R<sub>☉</sub>)</th><th>RV<sub>SINFONI</sub></th></tr></thead><tbody><tr><td rowspan="2">BT-Settl (G only, yes GP)</td><td>1637 ± 8</td><td>3.85 ± 0.03</td><td>9.3<sup>+0.03</sup><sub>-0.04</sub></td><td>0.2<sup>+0.2</sup><sub>-0.1</sub></td><td>0.4<sup>+0.3</sup><sub>-0.2</sub></td><td>...</td><td>...</td><td>1.0 ± 0.01</td><td>...</td></tr><tr><td>1477<sup>+8</sup><sub>-7</sub></td><td>3.93 ± 0.07</td><td>9.3 ± 0.03</td><td>0.3<sup>+0.1</sup><sub>-0.0</sub></td><td>0.6 ± 0.1</td><td>...</td><td>...</td><td>1.19 ± 0.03</td><td>...</td></tr><tr><td rowspan="2">BT-Settl (Si only, yes GP)</td><td>1624<sup>+8</sup><sub>-7</sub></td><td>3.89 ± 0.04</td><td>9.3 ± 0.03</td><td>0.4 ± 0.1</td><td>0.5 ± 0.1</td><td>...</td><td>...</td><td>0.99 ± 0.01</td><td>75395<sup>+1</sup><sub>-2</sub></td></tr><tr><td>1469<sup>+5</sup><sub>-4</sub></td><td>3.71<sup>+0.07</sup><sub>-0.06</sub></td><td>9.3 ± 0.03</td><td>0.4 ± 0.1</td><td>0.6 ± 0.1</td><td>...</td><td>...</td><td>1.24 ± 0.02</td><td>75396<sup>+1</sup><sub>-2</sub></td></tr><tr><td rowspan="2">BT-Settl (G only, no GP)</td><td>1639 ± 7</td><td>3.85 ± 0.03</td><td>9.3 ± 0.03</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1.0 ± 0.01</td><td>...</td></tr><tr><td>1487<sup>+7</sup><sub>-6</sub></td><td>3.96<sup>+0.05</sup><sub>-0.06</sub></td><td>9.3 ± 0.03</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1.17<sup>+0.03</sup><sub>-0.02</sub></td><td>...</td></tr><tr><td rowspan="2">BT-Settl (Si only, no GP)</td><td>1630 ± 7</td><td>3.92 ± 0.03</td><td>9.31 ± 0.03</td><td>...</td><td>...</td><td>...</td><td>...</td><td>0.97 ± 0.01</td><td>75389<sup>+1</sup><sub>-2</sub></td></tr><tr><td>1475<sup>+5</sup><sub>-4</sub></td><td>3.81<sup>+0.07</sup><sub>-0.06</sub></td><td>9.3 ± 0.03</td><td>...</td><td>...</td><td>...</td><td>...</td><td>1.21<sup>+0.02</sup><sub>-0.03</sub></td><td>75378 ± 1</td></tr><tr><td>Exo-REM (Si+G, yes GP)</td><td>1337 ± 9</td><td>3.52<sup>+0.03</sup><sub>-0.02</sub></td><td>9.3 ± 0.03</td><td>0.4 ± 0.1</td><td>0.5 ± 0.1</td><td>0.15<sup>+0.08</sup><sub>-0.1</sub></td><td>0.595<sup>+0.008</sup><sub>-0.009</sub></td><td>1.51<sup>+0.03</sup><sub>-0.02</sub></td><td>18 ± 8</td></tr></tbody></table>

**Note.** G is short for GRAVITY (i.e., the GRAVITY K-band spectrum was included in the fit), and Si is short for SINFONI (i.e., the SINFONI K-band spectrum was included in the fit). “yes GP” means that a Gaussian process (GP) regression was enabled for the SPHERE IFS spectrum, the hyperparameter constraints for which are reported in the table.

or GRAVITY spectrum, we only perform Exo-REM fits using both data sets.

The available Exo-REM grid is computed for a grid of metallicities, allowing us estimate the atmospheric metallicity and C/O ratio. However, the grid we used only computes predictions out to 5 µm, so we excluded the two longer-wavelength JWST/MIRI photometric data points for these fits. The results are shown in Figures 13 and 14, and summarized in Table 6. A solar metallicity with a C/O ratio of 0.6 is preferred. We find GP parameters for the SPHERE data set that are consistent with the results from the BT-Settl CIFIST fits (previous section), a good sanity check that the model is similarly treating correlated noise in both cases. We also recover a planetary RV value consistent with that reported in Petrus et al. (2021), which is an excellent sanity check given that their planetary RV was computed by cross-correlation with the SINFONI spectrum alone, and not through SED fitting. Unsurprisingly, our lower-resolution GRAVITY spectrum does not permit a direct RV measurement. The effective temperature recovered by this fit is about 150 K lower than the low-temperature BT-SETTL CIFIST mode, increasing our overall confidence in this low-temperature interpretation. However, the surface gravity value pushes against the low-gravity end of the grid, which casts doubt on the results of this fit.

The constraints we derived in this analysis are significantly more precise than those reported in Petrus et al. (2021), because at the time of that paper’s publication the Exo-REM grid predictions were only available for K band. We were able to use the updated Exo-REM grid to compare with all available spectral information. However, the uncertainties are likely underestimated, in particular because these fits do not account for interpolation errors, an important point that we discuss in more detail in the next section.

From scrutinizing the residuals shown in Figure 14, it is also clear that the NACO and JWST photometry beyond 3 µm is systematically higher than the Exo-REM models. Therefore, in addition to potentially underestimated errors due to unmodeled physics in the atmosphere grid, interpolation errors, and correlated noise in the SPHERE IFS spectrum, the best fit is still not perfect, again pointing to potential inaccuracies in our comparison. Future work could attempt to fit offsets between photometric data from different instruments in order to reduce the discrepancy, and/or inflate the errors in these or other

spectral data points. Altogether, the values and errors derived from the grid comparison should be treated with caution.

# 5. Discussion and Conclusion

## 5.1. Interpreting the Eccentricity Constraints

Because of the strong degeneracy between the eccentricity and inclination, it is most straightforward to report the constraints on these parameters in 2D. The direction of orbital motion on the plane of the sky constrains the inclination to >90°. For fit 4, including all available astrometry and applying a uniform eccentricity prior, the $1\sigma$ upper limit on both parameters is $e = 0.7/\text{inc} = 110°$, and the $2\sigma$ upper limit is $e = 0.8/\text{inc} = 120°$. Applying a nonuniform, linearly decreasing prior on the eccentricity, which previous work shows is appropriate for the CJ population (Bowler et al. 2020; Nagpal et al. 2023), tightens these upper limits even more. These are still tenuous constraints, but they are driven by the likelihood, not just by the prior. This is the first time we are obtaining eccentricity posteriors on this object that do not simply reproduce the prior. It is worth emphasizing the difficulty of measuring the eccentricity of an object almost 100 au from its star; previous studies of this object report that it would take 5–10 years of orbit monitoring before resolving the orbital curvature and constraining HIP 65426 b’s eccentricity. With the precision of VLTI/GRAVITY, we are able to “speed up time” and obtain eccentricity constraints sooner, with fewer measurements. Other directly imaged exoplanets with well-constrained eccentricities are generally much closer to their stars (e.g., β Pic b and c, at 3 and 10 au, respectively).

A potentially useful outcome of this paper is a generalizable prescription for interpreting the eccentricity posterior for incomplete orbits. To have confidence in an eccentricity measurement from a posterior, we suggest showing that the eccentricity posterior is:

1. prior independent, that is, driven by the likelihood and not phenomenologically different under different prior assumptions,

2. inconsistent with a circular orbit, ideally beyond $3\sigma$ for varied prior assumptions, and

3. not dependent on a single data point.

HIP 65426 b does not yet satisfy the first two criteria, so continued orbital monitoring with VLTI/GRAVITY will be
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Results of forward modeling the photometric and spectral data of HIP 65426 b by comparing with the BT-SETTL CIFIST model grid. Posteriors over the free parameters in the fit, as well radius, a derived parameter, are shown. Fits performed using GRAVITY K-band spectra are shown in purple, and fits performed using SINFONI K-band spectra are shown in pink. The GP hyperparameters (defined as in Equation (4) in Wang et al. 2020) to the SPHERE IFS spectral data (length scale and amplitude) are shown as well. Takeaways: As expected, log g correlates strongly with radius and T<sub>eff</sub>. Two families of solutions are apparent at high (1.3 R<sub>j</sub>) and low (0.9 R<sub>j</sub>) radii.</figcaption></figure>

important to measure its eccentricity robustly. Planetary RVs
with sub-kilometer per second precision (Figure 4), which
could be obtained with high-spectral-resolution instruments
like CRIRES, would also further constrain the eccentricity
(e.g., Snellen et al. 2014; Schwarz et al. 2016).

With these caveats in mind, the posteriors presented in this
paper favor a low or moderate eccentricity. Given the results of
Marleau et al. (2019), this is a preliminary hint that

HIP 65426 b did not attain its current position via scattering
after disk dispersal.

5.2. *Interpreting the Atmosphere Constraints*

In this analysis, we have focused on updating the atmo-
spheric model analysis, independent from the evolution-based
analysis of Carter et al. (2023), which constrained physical
<figure><img src="image_14.png" /><figcaption>**Figure 12.** BT-SETTL CIFIST models representing the two posterior peaks shown in Figure 11, together with the SPHERE IFS data. Top: both models, resampled onto the SPHERE IFS wavelength grid using spectres (Carnall 2017), and multiplied by a scalar chosen to minimize the sum of the squared residuals for the SPHERE IFS data alone. The SPHERE IFS data are shown as blue points, with error bars representing their reported statistical uncertainties. Middle: the low-$T_{\text{eff}}$ model (dashed purple line), SPHERE IFS data, and GP $1\sigma$ uncertainties (solid purple band; computed using the MAP GP parameters). Bottom: same as middle, but the high-$T_{\text{eff}}$ model is shown as a solid pink line, and GP uncertainties as a pink band. Takeaway: there are correlated residuals in the SPHERE passband for both of the $T_{\text{eff}}$ modes recovered from comparisons to the BT-Settl CIFIST grid, which we model with a squared-exponential GP.</figcaption></figure>

atmospheric properties using an estimate of $L_{\text{bol}}$. First, it is important to understand the limits of the spectral interpretation approach we have taken in this paper. Self-consistent grid modeling involves interpolating spectra in multiple dimensions during the spectral inversion, while the variations of the synthetic spectra along the grid are not linear in these dimensions (e.g., Czekala et al. 2015; Petrus et al. 2021). Missing or incorrect physics in the model grids, therefore, is not the only source of error. Performing atmospheric retrievals to evaluate the grid comparison results of this study is an important next step.

Keeping this limitation in mind, in this work we repeated the analysis presented in Petrus et al. (2021), which compared the available spectral and photometric data of HIP 65426 b with self-consistent BT-Settl and Exo-REM grids. This work benefited from additional data (in particular, the medium-resolution K-band GRAVITY spectrum; an additional NACO photometric point presented in Stolker et al. 2020; and the JWST photometry presented in Carter et al. 2023), and the expanded capability of Exo-REM to handle data outside of the K band. Like Petrus et al. (2021), we recovered two modes in our BT-Settl posterior fit, regardless of which K-band spectrum we use and whether we included a correlated noise model for the SPHERE IFS data: one at a higher radius of 1.2

$R_J$ and one at a lower radius of 1.0 $R_J$. Both modes are significantly below the hot-start radius of 1.4 $R_J$ derived in Carter et al. (2023), reinforcing the tension that paper originally pointed out.

Like Petrus et al. (2021), we only recover a single posterior mode when comparing with the Exo-REM grid, which is about 150 K cooler than the coolest MAP BT-Settl fit. The GP hyperparameters applicable to the SPHERE data sets are consistent across both grids, which we interpret as the presence of correlated observational noise, but which could also be a systematic problem common to both grids. The Exo-REM fit posterior favors a slightly supersolar metallicity and a C/O ratio of 0.6. However, the log g posterior hits the edge of the grid, so we strongly encourage skepticism of these derived values. The C/O ratio, metallicity, and other atmosphere properties we present here are consistent but more precise than those presented in Petrus et al. (2021), keeping in mind that the systematic errors are likely underestimated.

A planetary metallicity relative to its host star's is generally expected to reflect its formation condition (Öberg et al. 2011; Madhusudhan et al. 2014); in particular, formation via gravitational instability is generally held to produce an object with the same metallicity as its primary. In order to interpret HIP 65426 b's metallicity, it is important to understand the
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Exo-REM posterior fits to all data, showing 2D covariances and 1D marginalized posteriors over the fitted model parameters and derived parameters (radius, luminosity, and mass). Also see Figure 14. Note: parallax is denoted as $\omega$ here, and $\pi$ elsewhere in the text. Takeaways: The Exo-REM derived atmosphere parameters are about 150 K lower than the low-$T_{eff}$ BT-SETTL CIFIST parameters (Figure 11). The surface gravity hits the edge of the available grid. A slightly supersolar-metallicity atmosphere with a C/O ratio of 0.6 is favored, although there are likely systematic errors unaccounted for in this fit.</figcaption></figure>

primary's metallicity. Because HIP 65426 A is a fast rotator with few spectral lines, a direct metallicity measurement is difficult, but if we take the metallicity of other members of its moving group as representative, it should be approximately solar, modulo scatter among individual stars (see Section 4.3 of Petrus et al. 2021).

Taken with a big grain of salt, the supersolar metallicity of HIP 65426 b, expected ~solar metallicity of HIP 65426 A, and the planetary C/O ratio of 0.6, as constrained by comparison

with the Exo-REM grid, are consistent with formation via core accretion past the CO snowline (Öberg et al. 2011). The exact location of this snowline for HD 163296, which has a similar mass to HIP 65426 A, was observed to be 75 au (Qi et al. 2015; Petrus et al. 2021). This (very tentative, with many caveats, see, e.g., Mollière et al. 2022) picture provides a parallel constraint on the picture of how HIP 65426 b attained its current separation by setting an outer limit on the initial formation location before any scattering occurred.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Full best-fit model SED from comparison with the Exo-REM grid, together with all fitted spectral data. Top: transmission functions of the photometric bands. Middle: MAP model spectrum, along with 100 random draws from the posterior (in gray; barely distinguishable from the MAP spectrum). The spectral data are plotted as points with error bars. The horizontal bars of the photometry points indicate their spectral bandpasses. The corresponding band-integrated model predictions are overplotted as empty symbols. Bottom: MAP model residuals. Takeaways: Overall, the model spectrum fits well. Correlated noise is visible in the residuals of the SPHERE data set. The NACO and JWST/NIRCAM points beyond 3 µm are underestimated by the model.</figcaption></figure>

## 5.3. Future Directions

Continued orbital monitoring, both with VLTI/GRAVITY and spectrographs capable of measuring additional planetary RVs will refine the eccentricity measurement of HIP 65426 b over the next few years, allowing us more insight into this specific planet's formation and further constraining the population-level eccentricity distribution of CJs. Uncertainties of ~1 km s⁻¹ or below are needed for relative RV measurements of HIP 65426 b to constrain orbital parameters like a and e, motivating observation with high-resolution spectrographs like CRIRES. Other tracers of its formation condition that will be measurable in the near future include the planet's obliquity (Sepulveda et al. 2023), dynamical mass (hopefully measurable upon the release of Gaia time-series data), and spin (from high-resolution spectra). Further atmospheric characterization work, particularly retrievals that jointly model bulk atmosphere properties and trace chemical species fingerprints (see, e.g., Xuan et al. 2022) represent an important parallel path toward assessing the trustworthiness of the metallicity and C/O ratio, which will be helpful for pinpointing the formation location of HIP 65426 b. Xuan et al. (2022) showed that high-resolution spectra are sensitive to a broader range of atmospheric pressures than lower-resolution spectra, leading to more robust abundance measurements, motivating high-resolution measurements of HIP 65426 b in particular.

It is worth mentioning that the rapid rotation and early spectral type<fnref n="48" /> of the primary star, HIP 65426 A, precludes precise RV measurements, so a dynamical mass measurement will rely entirely on the combination of relative (i.e., from high-contrast imaging) and absolute (i.e., from Gaia) astrometry. However, RV monitoring of the planet, together with continued relative astrometric monitoring, may allow for the detection of unseen inner companions, through Keplerian effects on the host

star (Lacour et al. 2021) and/or planet-planet interactions (Covarrubias et al. 2022).

More theoretical work is also needed in order to interpret these results. In particular, population-level studies à la Marleau et al. (2019) could be conducted for alternate plausible formation pathways, particularly more rapid core formation via pebble accretion and gravitational instability in the protoplanetary disk. It would also be interesting to compare the existing measurements of other formation tracers, for example orbital inclination, metallicity, and C/O ratio, with their corresponding predictions in existing models.

VLTI/GRAVITY is a powerful instrument. This study has mostly focused on its ability to refine orbital eccentricity measurements in order to make dynamical inferences useful for commenting on planet formation. However, the ExoGRAVITY program is actively investigating a number of other scientific questions and observational constraints, particularly precise dynamical mass measurements and bolometric luminosities (e.g., Hinkley et al. 2023).

Planet formation is a complex process, and we will need a diverse set of observational and theoretical tools to unravel its secrets. This paper represents a small step toward better constraints and deeper understanding.

## Acknowledgments

We thank Aldo Sepulveda and Dan Huber for collaborative and informative discussions. S.B. and J.J.W. are supported by NASA Grant 80NSSC23K0280. G.-D.M. acknowledges the support of the DFG priority program SPP 1992 “Exploring the Diversity of Extrasolar Planets” (MA 9185/1) and from the Swiss National Science Foundation under grant 200021_204847 “PlanetsInTime.” Parts of this work have been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation. S.P. acknowledges the support of ANID, - Millennium Science

<footnote n="48">⁴⁸ Hot, early-type stars have few spectral lines, and rapid rotators have broad spectral lines, which both decrease RV precision.</footnote>
Initiative Program—NCN19 171. S.L. acknowledges the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-21-CE31-0017 (project ExoVLTI). This work is based on observations collected at the European Southern Observatory under ESO program 1104.C-0651. It also made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and the Jet Propulsion Laboratory/California Institute of Technology, funded by the National Aeronautics and Space Administration. S.B. wishes to acknowledge her status as a settler on the ancestral lands of the Gabrielino/Tongva people.

Software: tinygp (github.com/dfm/tinygp; Foreman-Mackey et al. 2022), jax (Bradbury et al. 2018), numpy (Harris et al. 2020), pandas (McKinney 2010), matplotlib (Hunter 2007), spectres (Carnall 2017), species (Stolker et al. 2020), corner (Foreman-Mackey 2016), scipy (Virtanen et al. 2020), astropy (Astropy Collaboration et al. 2013, 2018, 2022), and orbitize! (Blunt et al. 2020). Analysis code used to perform this work is available on GitHub with a copy deposited to Zenodo: 10.5281/zenodo.8429161 (Blunt & Balmer 2023).

## Appendix

Here, we further evaluate the differences between orbit fits 4, 6, and 7, which apply different priors on the eccentricity, using a model selection lens. There exist many model selection metrics a statistician can use to pick an “optimal” model. In this section, we briefly motivate the definitions of two such metrics, the AIC and WAIC, drawing heavily from Gelman et al. (2014), and use them to provide a different perspective on the explanation in Section 3.2.

We can choose to define a model’s “goodness” by its ability to predict unseen data points. A “good” model will predict new, unseen data better than a “bad” model. However, a posterior comprises many models, each with its own probability based on existing data, and Bayesian model selection seeks to compare two posteriors. In addition, we do not know a priori what we expect unseen data to look like; this is the purpose of model fitting. The logic Gelman et al. (2014) review is that one can define the “expected (log) predictive density” (elppd) for a new data point as the expected posterior probability of the new data point, weighted by the probability of the new data point itself (typically unknown):

$$ \text{elppd} = \int \log p_{\text{post}}(y_i) f(y_i) dy, \quad (\text{A1}) $$

where $p_{\text{post}}$ is the posterior probability, $y_i$ is an unseen data point, and $f(y)$ is the unknown physical function generating new data. Taking it one step further, we can aim to compute the expected log pointwise predictive density (elppd), which sums the elppd for each new data point over an arbitrary-sized unseen data set. The problem then becomes determining an estimator of this elppd quantity. (See Section 2.3 of Gelman et al. 2014 for more details.)

<table><caption>**Table 7**
Model Comparison Metrics for Orbit Fits Including All Available Astrometric Data and Varying the Eccentricity Prior</caption><thead><tr><th></th><th>Uniform e Prior</th><th>Linear e Prior</th><th>e = 0</th></tr></thead><tbody><tr><td>ΔAIC</td><td>+5.97</td><td>+5.97</td><td>0.00</td></tr><tr><td>ΔWAIC</td><td>+0.22</td><td>+0.04</td><td>0.00</td></tr></tbody></table>

**Note.** All values are computed relative to the $e = 0$ model. The linear e prior is given by Equation (1). For the AIC computation, the $M_{\text{tot}}$ and parallax were not included in the number of free parameters, since they were both highly constrained by their respective priors.

Both the AIC and the WAIC take the approach of defining the elppd of an unseen data set as the sum of two terms: the first term represents the model’s ability to predict existing data points and the second term is an overfitting penalty. They differ in their definitions of each of these terms. Both of these estimators converge to the actual log pointwise predictive density (lppd)<fnref n="49" /> under certain conditions.

The AIC computes the first term as the probability of the existing data given the maximum-likelihood model. The second term is simply the number of free parameters in the model. This definition has the pleasing property of being an unbiased estimator of the lppd for Gaussian posteriors which were computed using a flat prior. It may already be clear that we will choose to argue that the AIC is insufficient for the orbit described in this paper, as the orbital posteriors we seek to compare (e.g., Figure 6) are quite non-Gaussian. In addition, the effects of the prior volume are important for our posteriors, and the maximum-likelihood estimate (which is the only relevant quantity for the AIC) is different from the MAP estimate for both eccentricity priors discussed in Section 3.2. Because the maximum-likelihood estimate is approximately equal for the free-eccentricity and fixed-eccentricity models, while the number of free parameters differs, we can understand why the AIC metric favors the $e = 0$ model (Table 7).

The WAIC uses the whole posterior, not just the maximum-likelihood estimate, to compute both the first and second terms. Because of this, Gelman et al. (2014) call it a “a more fully Bayesian approach for estimating the out-of-sample expectation.” The first term is exactly the average predictive density for all existing data points (Equation 5 in Gelman et al. 2014), and the second term has a few definitions in the literature. In this article, we define the WAIC using Gelman et al. (2014; their Equation 12), which sets the overfitting term to be the sum of the variances of the lppd values of individual existing data points. Under this definition, a model which is overfitting will have (on average) smaller variances in the predicted posterior probabilities of existing data points. Because the WAIC includes information from the whole posterior, implicitly taking into account the effects of the prior volume, we might expect that the WAIC is a better estimator for our thoroughly non-Gaussian posteriors. Indeed, the WAIC shows a much less clear distinction between the various eccentricity models than the AIC (Table 7). The circular model is still preferred, but by a much smaller margin than the AIC suggests (and by an amount which many authors argue is “essentially indistinguishable”). This motivates us to argue that, although the AIC prefers the

<footnote n="49">⁴⁹ The elppd is an estimator of a true underlying value, the lppd. The AIC and WAIC are two different definitions of elppd estimators.</footnote>
fixed *e* = 0 model, the three models are indistinguishable in terms of expected predictive ability. In other words, we do not rule out a moderate eccentricity, despite the AIC preference for a circular orbit, and determine that there is not enough evidence to select conclusively one of the three models as “optimal.” The preference for a low or moderate eccentricity over a high eccentricity (≥0.8), however, is driven by the lower likelihood at higher eccentricities, and is consistent across all three prior options.

## ORCID iDs

Sarah Blunt https://orcid.org/0000-0002-3199-2888

W. O. Balmer https://orcid.org/0000-0001-6396-8439

J. J. Wang https://orcid.org/0000-0003-0774-6502

S. Lacour https://orcid.org/0000-0002-6948-0263

J. Kammerer https://orcid.org/0000-0003-2769-0438

E. Rickman https://orcid.org/0000-0003-4203-9715

J. Shangguan https://orcid.org/0000-0002-4569-9009

A. Amorim https://orcid.org/0000-0003-2990-0726

M. Benisty https://orcid.org/0000-0002-7695-7605

A. Bohn https://orcid.org/0000-0003-1401-9952

M. Bonnefoy https://orcid.org/0000-0001-5579-5339

W. Brandner https://orcid.org/0000-0003-1939-6351

F. Cantalloube https://orcid.org/0000-0002-3968-3780

P. Caselli https://orcid.org/0000-0003-1481-7911

B. Charnay https://orcid.org/0000-0003-0977-6545

V. Christiaens https://orcid.org/0000-0002-0101-8814

Y. Clénet https://orcid.org/0000-0002-8382-2020

A. Eckart https://orcid.org/0000-0001-6049-3132

R. Genzel https://orcid.org/0000-0002-2767-9653

S. Gillessen https://orcid.org/0000-0002-5708-0481

J. H. Girard https://orcid.org/0000-0001-8627-0404

Th. Henning https://orcid.org/0000-0002-1493-300X

S. Hinkley https://orcid.org/0000-0001-8074-2562

S. Hippler https://orcid.org/0000-0002-3912-6108

M. Horrobin https://orcid.org/0000-0001-5451-7847

P. Kervella https://orcid.org/0000-0003-0626-1749

L. Kreidberg https://orcid.org/0000-0003-0514-1147

D. Lutz https://orcid.org/0000-0003-0291-9582

A.-L. Maire https://orcid.org/0000-0002-2591-4138

G.-D. Marleau https://orcid.org/0000-0002-2919-7500

A. Mérand https://orcid.org/0000-0003-2125-0183

P. Mollière https://orcid.org/0000-0003-4096-7067

J. D. Monnier https://orcid.org/0000-0002-3380-3307

C. Mordasini https://orcid.org/0000-0002-1013-2811

T. Ott https://orcid.org/0000-0003-1572-0396

T. Paumard https://orcid.org/0000-0003-0655-0452

G. Perrin https://orcid.org/0000-0003-0680-0167

L. Rodet https://orcid.org/0000-0002-1259-3312

Z. Rustamkulov https://orcid.org/0000-0003-4408-0463

T. Shimizu https://orcid.org/0000-0002-2125-4670

D. Sing https://orcid.org/0000-0001-6050-7645

T. Stolker https://orcid.org/0000-0002-5823-3072

C. Straubmeier https://orcid.org/0000-0002-0671-9302

E. Sturm https://orcid.org/0000-0002-0183-3666

L. J. Tacconi https://orcid.org/0000-0002-1485-9401

A. Vigan https://orcid.org/0000-0002-5902-7828

F. Vincent https://orcid.org/0000-0002-3855-0708

K. Ward-Duong https://orcid.org/0000-0002-4479-8291

J. Woillez https://orcid.org/0000-0002-2958-4738

## References

Allard, F., Allard, N. F., Homeier, D., et al. 2007, A&A, 474, L21

Allard, F., Guillot, T., Ludwig, H.-G., et al. 2003, in IAU Symp. 211, Brown Dwarfs, ed. E. Martín (San Francisco, CA: ASP), 325

Allard, F., Homeier, D., & Freytag, B. 2011, in ASP Conf. Ser. 448, 16th Cambridge Workshop on Cool Stars, Stellar Systems, and the Sun, ed. C. Johns-Krull, M. K. Browning, & A. A. West (San Francisco, CA: ASP), 91

Allard, F., Homeier, D., & Freytag, B. 2012, RSPTA, 370, 2765

Armitage, P. J. 2020, Astrophysics of Planet Formation (2nd edn.; Cambridge: Cambridge Univ. Press)

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33

Blunt, S., & Balmer, W. O. 2023, sblunt/hip65426, v0.1, Zenodo, doi:10.5281/ZENODO.8429161

Blunt, S., Endl, M., Weiss, L. M., et al. 2019, AJ, 158, 181

Blunt, S., Wang, J., Ngo, H., et al. 2023, sblunt/orbitize: tests overwrite previously existing file bugfix, v2.2.2, Zenodo, doi:10.5281/zenodo.8102085

Blunt, S., Wang, J. J., Angelo, I., et al. 2020, AJ, 159, 89

Bochanski, J. J., Faherty, J. K., Gagné, J., et al. 2018, AJ, 155, 149

Bowler, B. P. 2016, PASP, 128, 102001

Bowler, B. P., Blunt, S. C., & Nielsen, E. L. 2020, AJ, 159, 63

Bowler, B. P., & Nielsen, E. L. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte (Berlin: Springer), 155

Bradbury, J., Frostig, R., Hawkins, P., et al. 2018, JAX: composable transformations of Python+NumPy programs v0.3.13., http://github.com/google/jax

Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125

Carnall, A. C. 2017, arXiv:1705.05165

Carter, A. L., Hinkley, S., Kammerer, J., et al. 2023, ApJL, 951, L20

Charnay, B., Bézard, B., Baudino, J. L., et al. 2018, ApJ, 854, 172

Chauvin, G., Desidera, S., Lagrange, A. M., et al. 2017, A&A, 605, L9

Cheetham, A. C., Samland, M., Brems, S. S., et al. 2019, A&A, 622, A80

Coleman, G. A. L., & Nelson, R. P. 2016a, MNRAS, 457, 2480

Coleman, G. A. L., & Nelson, R. P. 2016b, MNRAS, 460, 2779

Covarrubias, S., Blunt, S., & Wang, J. J. 2022, RNAAS, 6, 66

Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, yCat, II/246

Czekala, I., Andrews, S. M., Mandel, K. S., Hogg, D. W., & Green, G. M. 2015, ApJ, 812, 128

Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449

Feroz, F., Hobson, M. P., & Bridges, M. 2009, MNRAS, 398, 1601

Feroz, F., Hobson, M. P., Cameron, E., & Pettitt, A. N. 2019, OJAp, 2, 10

Ferrer-Chávez, R., Wang, J. J., & Blunt, S. 2021, AJ, 161, 241

Foreman-Mackey, D. 2016, JOSS, 1, 24

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Foreman-Mackey, D., Yadav, S. T., et al. 2022, dfm/tinygp: The tiniest of Gaussian Process libraries, v0.2.3, Zenodo, doi:10.5281/ZENODO.7269074

Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv:2208.00211

Gagné, J., Mamajek, E. E., Malo, L., et al. 2018, ApJ, 856, 23

Gelman, A., Hwang, J., & Vehtari, A. 2014, Statistics and Computing, 24, 997

Gravity Collaboration, Abuter, R., Accardo, M., et al. 2017, A&A, 602, A94

Gravity Collaboration, Nowak, M., Lacour, S., et al. 2020, A&A, 633, A110

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357

Hinkley, S., Lacour, S., Marleau, G. D., et al. 2023, A&A, 671, L5

Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27

Householder, A., & Weiss, L. 2022, arXiv:2212.06966

Hunter, J. D. 2007, CSE, 9, 90

Lacour, S., Dembet, R., Abuter, R., et al. 2019, A&A, 624, A99

Lacour, S., Nowak, M., Wang, J. J., et al. 2019, A&A, 624, A99

Lacour, S., Wang, J. J., Nowak, M., et al. 2020, Proc. SPIE, 11446, 1144600

Lacour, S., Wang, J. J., Rodet, L., et al. 2021, A&A, 654, L2

Lambrechts, M., & Johansen, A. 2014, A&A, 572, A107

Lapeyrère, V., Kervella, P., Lacour, S., et al. 2014, Proc. SPIE, 9146, 91462D

Madhusudhan, N., Amin, M. A., & Kennedy, G. M. 2014, ApJL, 794, L12

Marleau, G.-D., Coleman, G. A. L., Leleu, A., & Mordasini, C. 2019, A&A, 624, A20
McKinney, W. 2010, in Proc. of the 9th Python in Science Conf., ed. Stéfan van der Walt & Jarrod Millman (Austin, TX: SciPy), 56

Mollière, P., Molyarova, T., Bitsch, B., et al. 2022, ApJ, 934, 74

Nagpal, V., Blunt, S., Bowler, B. P., et al. 2023, AJ, 165, 32

Nayakshin, S. 2017, MNRAS, 470, 2387

Nielsen, E. L., Close, L. M., Biller, B. A., Masciadri, E., & Lenzen, R. 2008, ApJ, 674, 466

Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJL, 743, L16

Petrus, S., Bonnefoy, M., Chauvin, G., et al. 2021, A&A, 648, A59

Qi, C., Öberg, K. I., Andrews, S. M., et al. 2015, ApJ, 813, 128

Rosenthal, M. M., & Murray-Clay, R. A. 2018, ApJ, 864, 66

Schwarz, H., Ginski, C., de Kok, R. J., et al. 2016, A&A, 593, A74

Seager, S. 2010, Exoplanets (Tucson, AZ: Univ. Arizona Press)

Sepulveda, A., Huber, D., Bedding, T., et al. 2023, AAS Meeting Abstracts, 55, 164.02

Snellen, I. A. G., Brandl, B. R., de Kok, R. J., et al. 2014, Natur, 509, 63

Stolker, T., Quanz, S. P., Todorov, K. O., et al. 2020, A&A, 635, A182

Vigan, A., Fontanive, C., Meyer, M., et al. 2021, A&A, 651, A72

Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261

Vorobyov, E. I., & Elbakyan, V. G. 2018, A&A, 618, A7

Vousden, W. D., Farr, W. M., & Mandel, I. 2016, MNRAS, 455, 1919

Wang, J. J., Ginzburg, S., Ren, B., et al. 2020, AJ, 159, 263

Wang, J. J., Graham, J. R., Dawson, R., et al. 2018, AJ, 156, 192

Wang, J. J., Kulikauskas, M., & Blunt, S., 2021 whereistheplanet: Predicting Positions of Directly Imaged Companions, Astrophysics Source Code Library, ascl:2101.003

Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868

Xuan, J. W., Wang, J., Ruffio, J.-B., et al. 2022, ApJ, 937, 54