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

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 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.


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 (2016b,a).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 widelyseparated regime.Marleau et al. (2019) derived distributions over the initial mass and luminosity of HIP 65426 b under several assumptions of post-formation 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 post-scattering 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-tomoderate 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 priordependent (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. 2022).In this paper, we update the orbit model of HIP 65426 b using new astrometry from the optical interferometer VLTI/GRAVITY (∼50x 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 the HIP 65426 b measurements with BT-Settl CIFIST models, with one small radius (∼1R J ) peak and one larger radius peak (∼1.2RJ ).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 eff =1560 ± 100K, log g < 4.40, a slightly super-solar metallicity of 0.05 +0.24  −0.22 , and an upper limit on the C/O ratio (≤ 0.55).Carter et al. (2022) subsequently published 7 long-wavelength photometric datapoints of HIP 65426 b using the NIRCAM and MIRI instruments on JWST, and used them to update the empirical bolometric luminosity and of the planet to log L/L ⊙ = (-4.31, -4.14).Together with hot-start cooling models, they used this luminosity estimate to infer a 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 ± 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 photometry and spectra measurements of HIP 65426 with self-consistent model spectra to update our understanding of the planet's atmospheric properties, finding 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.Observing log.NEXP, NDIT, and DIT denote the number of exposures, the number of detector integrations per exposure, and the detector integration time, respectively.τ0 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.

GRAVITY Data
We observed HIP 65426 b on the 7th of January, 2021, the 23rd of January, 2022, and the 7th of May, 2023 as part of the ExoGRAVITY Large Program (ESO Program ID 1104.C-0651, Lacour et al. 2020).We used the European Southern Observatory (ESO) Very Large Telescope Interferometer (VLTI)'s four 8.2m 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 star light 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-01-07 we used the binary system HD 91881, on 2022-01-23 we used two binary systems: HD 73900 and HD 91881, and on 2023-05-07 we used HD 123227.During the night of the 2021-01-07, we also performed an observation "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 whereistheplanet1 software (Wang et al. 2021), and resulted in an efficient coupling of the planet flux into the science fiber.
For the 2021-01-07 epoch, we calculated the complex visibilites of the host and the companion using the Public Release 1.5.0 (1 July 20212 ) of the ESO GRAVITY pipeline (Lapeyrère et al. 2014).The observations were phasereferenced 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 A, 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 analysing 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 χ 2 value of this map.We estimated the uncertainty on each astrometric measurement using the scatter of mean astrometric values between individual exposures.These new astrometric datapoints 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 microns.We only extract a spectrum for the observations from 2021-01-07, 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 eff = 8600 K (Carter et al. 2022) 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), 2MASS (Cutri et al. 2003), and 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-01-07 observation, F K2 = 1.40 × 10 −16 W/m 2 /µm vs the SPHERE value F K2 = 7.21 × 10 −17 W/m 2 /µm.This is likely due to two factors: the low integration time on the star (< 3 minutes) and the degradation/variability of 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.

Literature Data
Table 4. Astrometric measurements from the literature used in the orbit fits presented in this paper.P.A. is the position angle.Here, RV pl indicates a measurement of the planet's radial velocity relative to the primary.Making this measurement involved separately measuring the absolute radial velocity of the planet (from a SINFONI spectrum) and the absolute radial velocity of the star (from a series of HARPS spectra), subtracting these quantities, and propagating the uncertainty.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. (2022), 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).

ORBIT ANALYSIS
In order to understand how our new GRAVITY data constrains 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 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 (ω p ), position angle of nodes (Ω), and epoch of periastron, defined as fraction of the orbital period past a specified reference date (see Blunt et al. (2020) section 2.1 for a bit more detail).We use mjd=58859 as the reference date in this paper.Following the suggestion in Householder & Weiss (2022), we also explicitly specify here that we assume that ω 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 system5 .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 radial velocities.The ability to fit companion radial velocities 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 parameters (a, e, etc.) as those of the relative orbit, meaning that astrometric and radial velocity measurements of the secondary are both assumed to be relative to the primary. 6In 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. 7Relative companion RV measurements like this do not allow us to measure a dynamical mass for the companion, but have the potential to reduce posterior uncertainties of the relative orbit.However, the uncertainty on the available companion RV measurement is too large to meaningfully constrain the orbital parameters 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 6 When radial velocity measurements are available for both the star and the planet, orbitize!instead assumes that planetary and stellar RV measurements are relative to the system barycenter. 7It 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.
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 zero8 .7. Fit 4, except applying a decreasing prior on eccentricity.
(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 eccentricityinclination degeneracy (see, e.g., Ferrer-Chávez et al. 2021).
For all fits, we used the ptemcee implementation of emcee (Vousden et al. 2016;Foreman-Mackey et al. 2013), a parallel-tempered affine-invariant 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 100th step was saved. 10he 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.
We applied Gaussian priors on 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.

The Impact of GRAVITY
Figures 5 and 6 visualize the posterior of fit #4, using all astrometric data and assuming a uniform prior on eccentricity, and Figure 7 compares each of the orbit fits that use varying data subsets.These plots show the transformative impact of 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 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 evident 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.

A Moderate Eccentricity?
All of the orbit fits that include GRAVITY data show a preference for moderate eccentricities (e ∼ 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 non-zero 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 achieveable 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, likelihood is high, but 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, prior volume is high, and likelihood is also high, leading to a posterior peak.At high eccentricities (e > 0.6), prior volume is very high, but 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 radial velocity orbit).We can explain the maximum a posteriori (MAP) shift to lower eccentricities when using a linearly decreasing prior, therefore, 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.The GRAVITY K-band spectrum (R=500) published in this work overlaps in wavelength coverage almost completely with the higher-resolution SINFONI spectrum (R∼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  4) only GRAVITY astrometry (i.e.no literature data), and (5) all astrometric data (i.e. the accepted fit; purple outline).Takeaway: 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.Table 5. Marginalized posterior 68% credible intervals for the free parameters in each of the orbit fits described in Section 3. Parallax and total mass were additional free parameters in every.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⊙. 153.8 +5.9

Fit
−2.5 0.46 ± 0.14 Table 6.68% credible intervals of posterior fits to self-consistent model spectra grids.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 regression was enabled for the SPHERE IFS spectrum, the hyperparameter constraints for which are reported in the table.
fit name Although the maximum a posteriori 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.
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 2011c11 (Allard et al. 2003, Allard et al. 2007, Allard et al. 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) ability to explore non-solar metallicities, and 2) the success of the Exo-REM 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, Feroz et al. 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 mininum detectable vsini of ∼50km s −1 , which translates to a rotation period of ∼0.1d, assuming i=90 and R=1.5R J , so the non-detection of rotational broadening is also consistent with our physical expectation.

BT-SETTL CIFIST
We performed four variations of BT-SETTL CIFIST comparisons to our full spectral dataset, by 1) using either the GRAVITY or SINFONI spectrum12 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. 2022), 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 dataset 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.

Exo-REM
Because the agreement between the BT-SETTL CIFIST model fits is good regardless of whether we use the SINFONI or GRAVITY spectrum, we only perform Exo-REM fits using both datasets.
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 longerwavelength JWST/MIRI photometry 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 dataset 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 radial velocity 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 150K 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 un-modeled 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 datapoints.Altogether, the values and errors derived from grid comparison should be treated with caution.Because of the strong degeneracy between eccentricity and inclination, it is most straightforward to report constraints on these parameters in two dimensions.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-σ upper limit on both parameters is e = 0.7/inc=110 • , and the 2-σ upper limit is e = 0.8/inc=120 • .Applying a non-uniform, linearly decreasing prior on eccentricity, which previous work shows is appropriate for the cold Jupiter 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 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).
A potentially useful outcome of this paper is a generalizable prescription for interpreting eccentricity posterior for incomplete orbits.To have confidence in an eccentricity measurement from a posterior, we suggest showing that the eccentricity posterior is:  4) 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 eff .Two families of solutions are apparent at high (1.3RJ ) and low (0.9 RJ ) radii.
• prior-independent, that is, driven by the likelihood and not phenomenologically different under different prior assumptions, • inconsistent with a circular orbit, ideally beyond 3-σ for varied prior assumptions, and   11).The surface gravity hits the edge of the available grid.A slightly super-solarmetallicity atmosphere with a C/O ratio of 0.6 is favored, although there are likely systematic errors unaccounted for in this fit.
• 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 important to robustly measure its eccentricity.Planetary RVs with sub-km s −1 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.

Interpreting the Atmosphere Constraints
In this analysis, we have focused on updating the atmospheric model analysis, independent from the evolution-based analysis of Carter et al. (2022), which constrained physical atmospheric properties using an estimate of L 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 the 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, andthe JWST photometery presented in Carter et al. 2022), and the expanded capability of Exo-REM to handle data outside of 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. (2022), 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 150K cooler than the coolest MAP BT-Settl fit.The Gaussian process hyperparameters applicable to the SPHERE datasets 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 super-solar 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( Öberg et al. , 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 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 super-solar 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.

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 cold Jupiters.Uncertainties of ∼ km s −1 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 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 timeseries 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 13 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, radial velocity 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 luminiosities (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.
Table 7. Model comparison metrics for orbit-fits including all available astrometric data and varying the eccentricity prior.All values are computed relative to the e = 0 model.The linear e prior is given by Equation (1).For the AIC computation, Mtot and parallax were not included in the number of free parameters, since they were both highly constrained by their respective priors.
Uniform e prior linear e prior e = 0 ∆AIC +5.97 +5.97 0.00 ∆WAIC +0.22 +0.04 0.00 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.

Figure 1 .
Figure 1.Detections of HIP 65426 b with VLTI/GRAVITY in epochs 1 and 2. Both periodogram power maps visualize the χ 2 fit to the interferometric observables assuming a point source, after removing the contribution of the star using a 4th-order polynomial.The outer dashed grey 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).

Figure 3 .
Figure 3. Best-fit SED model of HIP 65426 A used to convert our Top: transmission functions of photometric bands (each 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.

Figure 4 .
Figure 4. Planetary radial velocity predictions from the accepted fit (purple histogram), together with the planetary RV measurement from Petrus et al. (2021) (dashed black line, with 1-σ range shaded grey).Takeaway: The planetary RV measurement does not constrain the orbital parameters, beyond breaking the 180 • degeneracy for Ω and ω.A planetary RV of −3 km s −1 is allowed, given the astrometry alone, but is disfavored (relative to a planetary RV of 3 km s −1 ) because of the relative RV measurement.

Figure 5 .
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 (grey) 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 (grey), together with earlier (bottom) and later (top) astrometric measurements taken with VLTI/GRAVITY.1-and 2-σ error ellipses are shaded in dark and light pink, respectively.Takeaway: The VLTI/GRAVITY epochs are ∼50x more precise than existing astrometric measurements of HIP 65426, and reduce the posterior uncertainty.

Figure 6 .
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-diagonals show 2D covariances between posterior elements.1, 2, and 3-σ contours are outlined in the covariance panels, and individual posterior samples outside of the 3-σ boundaries are plotted directly as black dots.Takeaway: the 1D marginalized posterior distributions of semimajor axis and inclination are well constrained.Strong covariances are apparent, in particular between eccentricity, inclination, and semimajor axis.

Figure 7 .
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; grey), (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).Takeaway: 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.

Figure 8 .
Figure 8. 1D marginalized eccentricity posteriors for fits with uniform (purple) and linearly decreasing (pink) priors on eccentricity.The priors themselves are also plotted as lines of the same colors.Takeaway: 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.

Figure 9 .
Figure9.Maximum log(likelihood) as a function of eccentricity.Although the maximum a posteriori 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 (Figure8); 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.

Figure 10 .
Figure 10.GRAVITY and SINFONI K-band spectra comparison.Top: GRAVITY (grey) and SINFONI (purple) 1σ 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 datasets is excellent.

Figure 11 .
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 Wang et al. (2020) Equation4) 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 eff .Two families of solutions are apparent at high (1.3RJ ) and low (0.9 RJ ) radii.

Figure 12 .
Figure12.BT-SETTL CIFIST models representing the two posterior peaks shown in Figure11, 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 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 eff model (dashed purple line), SPHERE IFS data, and Gaussian Process 1−σ uncertainties (solid purple band; computed using the MAP GP parameters).Bottom: same as middle, but the high-T 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 eff modes recovered from comparisons to the BT-Settl CIFIST grid, which we model with a squared-exponential Gaussian process.

Figure 13 .
Figure 13.Exo-REM posterior fits to all data, showing 2D covariances and 1D marginalized posteriors over fitted model parameters and derived parameters (radius, luminosity, and mass).Also see Figure14.Note: parallax is denoted ϖ here, and π elsewhere in the text.Takeaways: the Exo-REM derived atmosphere parameters are about 150K lower than the low-T eff BT-SETTL CIFIST parameters (Figure11).The surface gravity hits the edge of the available grid.A slightly super-solarmetallicity atmosphere with a C/O ratio of 0.6 is favored, although there are likely systematic errors unaccounted for in this fit.

FFigure 14 .
Figure 14.Full best-fit model SED from comparison with the Exo-REM grid, together with all fitted spectral data.Top: transmission functions of photometric bands.Middle: MAP model spectrum, along with 100 random draws from the posterior (in grey; barely distinguishable from MAP spectrum).The spectral data are plotted as points with error bars.The horizontal bars of the photometry points indicate their spectral bandpass.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 dataset.The NACO and JWST/NIRCAM points beyond 3 µm are underestimated by the model.

Table 1 .
Relevant physical properties of HIP 65426 A and b