VLTI/GRAVITY Provides Evidence the Young, Substellar Companion HD 136164 Ab Formed Like a “Failed Star”

Young, low-mass brown dwarfs orbiting early-type stars, with low mass ratios (q ≲ 0.01), appear to be intrinsically rare and present a formation dilemma: could a handful of these objects be the highest-mass outcomes of “planetary” formation channels (bottom up within a protoplanetary disk), or are they more representative of the lowest-mass “failed binaries” (formed via disk fragmentation or core fragmentation)? Additionally, their orbits can yield model-independent dynamical masses, and when paired with wide wavelength coverage and accurate system age estimates, can constrain evolutionary models in a regime where the models have a wide dispersion depending on the initial conditions. We present new interferometric observations of the 16 Myr substellar companion HD 136164 Ab (HIP 75056 Ab) made with the Very Large Telescope Interferometer (VLTI)/GRAVITY and an updated orbit fit including proper motion measurements from the Hipparcos–Gaia Catalog of Accelerations. We estimate a dynamical mass of 35 ± 10 M J (q ∼ 0.02), making HD 136164 Ab the youngest substellar companion with a dynamical mass estimate. The new mass and newly constrained orbital eccentricity (e = 0.44 ± 0.03) and separation (22.5 ± 1 au) could indicate that the companion formed via the low-mass tail of the initial mass function. Our atmospheric fit to a SPHINX M-dwarf model grid suggests a subsolar C/O ratio of 0.45 and 3 × solar metallicity, which could indicate formation in a circumstellar disk via disk fragmentation. Either way, the revised mass estimate likely excludes bottom-up formation via core accretion in a circumstellar disk. HD 136164 Ab joins a select group of young substellar objects with dynamical mass estimates; epoch astrometry from future Gaia data releases will constrain the dynamical mass of this crucial object further.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.lowest-mass "failed binaries" (formed via disk fragmentation or core fragmentation)?Additionally, their orbits can yield model-independent dynamical masses, and when paired with wide wavelength coverage and accurate system age estimates, can constrain evolutionary models in a regime where the models have a wide dispersion depending on the initial conditions.We present new interferometric observations of the 16 Myr substellar companion HD 136164 Ab (HIP 75056 Ab) made with the Very Large Telescope Interferometer (VLTI)/GRAVITY and an updated orbit fit including proper motion measurements from the Hipparcos-Gaia Catalog of Accelerations.We estimate a dynamical mass of 35 ± 10 M J (q ∼ 0.02), making HD 136164 Ab the youngest substellar companion with a dynamical mass estimate.The new mass and newly constrained orbital eccentricity (e = 0.44 ± 0.03) and separation (22.5 ± 1 au) could indicate that the companion formed via the low-mass tail of the initial mass function.Our atmospheric fit to a SPHINX M-dwarf model grid suggests a subsolar C/O ratio of 0.45 and 3 × solar metallicity, which could indicate formation in a circumstellar disk via disk fragmentation.Either way, the revised mass estimate likely excludes bottom-up formation via core accretion in a circumstellar disk.HD 136164 Ab joins a select group of young substellar objects with dynamical mass estimates; epoch astrometry from future Gaia data releases will constrain the dynamical mass of this crucial object further.

Introduction
The origins of low-mass brown dwarf companions to stars are uncertain.These objects could form from the low-mass tail of the initial mass function (IMF; see, e.g., Chabrier 2003) via fragmentation of a molecular cloud, and be captured into binary orbits (e.g., Padoan & Nordlund 2004;Boyd & Whitworth 2005;Bate 2009Bate , 2012)), or they could form via the gravitational instability of a circumstellar disk (Boss 1997;Stamatellos et al. 2007;Stamatellos & Whitworth 2009;Kratter et al. 2010;Forgan & Rice 2013;Forgan et al. 2015Forgan et al. , 2018)).Some even argue that a handful of these objects might form via core accretion like planets, and undergo runaway accretion (D'Angelo et al. 2010;Mollière & Mordasini 2012;Bodenheimer et al. 2013) to reach masses greater than the deuterium burning limit (;13 M J ; Saumon et al. 1996;Chabrier & Baraffe 2000).Schlaufman (2018) argues to the contrary, showing that planetary mass companions with masses below 4 M J orbit stars with on average higher metallicity (associated with a larger amount of material for core-accretion formation), while planetary mass companions with masses above this cutoff do not necessarily orbit metal-rich stars.This would appear to indicate there is population-level evidence for distinct formation channels for exoplanets and brown dwarf companions with masses above and below 4 M J .
There are a number of curious systems that blur these observational rules of thumb that enforce the boundaries between formation channels, for instance the HD 206893 system, where two nearly coplanar low-mass companions reside within a debris disk (Delorme et al. 2017;Milli et al. 2017;Grandjean et al. 2019;Kammerer et al. 2021;Ward-Duong et al. 2021;Hinkley et al. 2023).For newly discovered low-mass companions, especially low-mass companions to early-type stars (which have low mass ratios), this confusion is inevitable.However, further orbital characterization can provide clues to their formation (Bowler et al. 2020), and spectral characterization can provide a further line of evidence.
Additionally, evolutionary and spectral models of substellar objects (brown dwarfs and massive giant planets), like their higher-mass stellar counterparts (e.g., Rodet et al. 2018;Dieterich et al. 2021), need to be tested and refined based on model-independent measurements of their masses derived from orbital characterization (Dupuy & Liu 2017;Brandt et al. 2019;Fontanive et al. 2019;Rickman et al. 2020;Vrijmoet et al. 2020;Brandt et al. 2021;Bonavita et al. 2022;Franson et al. 2022;Rickman et al. 2022;Vrijmoet et al. 2022;Balmer et al. 2023;Franson & Bowler 2023).The difficulty is that substellar objects amenable to these measurements, unlike stellar binaries, appear to be intrinsically rare (e.g., Dieterich et al. 2012;Fontanive et al. 2018;Nielsen et al. 2019;Duchêne et al. 2023), and generally substellar companions are orders of magnitude fainter than their hosts, necessitating a greater observational expense than stellar binary orbits.So far, the amount of substellar objects with dynamically derived masses is low, about 20 in total, and the number of these with accurate age estimates younger than a few hundred megayears is of order unity (see Figure 4 in Franson & Bowler 2023).The youngest brown dwarf with a dynamical mass, PZ Tel b ( - + M 27 9 25 J , 24 ± 3 Myr; Franson & Bowler 2023), has been predicted by evolutionary models to have a higher mass on average than its dynamical mass.Older, cold brown dwarfs like HD 19467 B have been found to be underluminous compared to evolutionary models (Brandt et al. 2021;Greenbaum et al. 2023).It remains to be seen to what extent there are systematic discrepancies between dynamical and evolutionary models with age, as more dynamical masses need to be estimated for each age category, particularly for young ages.
The sensitivity of direct imaging to higher intrinsic luminosities has revealed many faint, substellar companions to stars in young stellar associations, like the greater Scorpius-Centarus (Sco-Cen) association (Janson et al. 2012;Hinkley et al. 2015;Cheetham et al. 2018;Bohn et al. 2022).These discoveries present an excellent opportunity to determine and study the dynamical masses of young brown dwarfs.
The system has a low Gaia renormalized unit weight error (RUWE; 1.053), but has a strong proper-motion anomaly between Hipparcos and Gaia, as noted in Kervella et al. (2019Kervella et al. ( , 2022)), and a χ 2 = 15.9 in the Hipparcos-Gaia Catalog of Accelerations (HGCA; Brandt 2021), indicative of the presence of a perturbing companion.Kouwenhoven et al. (2007) and Wagner et al. (2020) identify a nearby low-mass star (HD 136164 AB, M B ∼ 0.3 M e , Ks mag = 11.2) at a separation of 5 2 (∼650 au), which cannot reproduce the observed astrometric signal (Kervella et al. 2022).Wagner et al. (2020) first identified the substellar companion HD 136164 Ab as part of their coronagraphic survey of A-type stars in Sco-Cen (Wagner et al. 2022).At a separation of 125-150 mas, the companion is consistent with inducing the astrometric signal observed between Hipparcos and Gaia; at a separation of ∼30 au the acceleration is consistent with a substellar mass (Kervella et al. 2022).The near-infrared magnitudes are consistent with the companion being substellar (Wagner et al. 2020, give 25 ± 5 M J based on evolutionary models from Baraffe et al. 2003), andWagner et al. (2020) report an estimate of the companion's spectral type (M6-L2) consistent with the observed absorption due to water in lowresolution Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE) spectrophotometry.Only considering two epochs of astrometry, they partially constrain the separation of the companion's orbit (∼30 ± 15 au), but were unable to estimate the eccentricity reliably.
With a low mass ratio and close separation, Wagner et al. (2020) posit that the companion could have formed via disk fragmentation (e.g., Boss 1997;Kratter et al. 2010;Forgan et al. 2018), or even the very high end of a "planetary" coreaccretion formation channel (e.g., Pollack et al. 1996;Mordasini et al. 2012) although the possibility (e.g., Emsenhuber et al. 2021) and occurrence (e.g., Schlaufman 2018) of such high-mass core-accretion planets is still heavily debated.Placing the companion's orbit in context could provide insight into this question, as recent population-level studies of substellar companions indicate two populations in eccentricity (Bowler et al. 2020;Doó et al. 2023;Nagpal et al. 2023).It appears that low-mass objects (that likely formed within a protoplanetary disk, whose gas damped their initial eccentricity) exhibit a distribution of eccentricities that is on average lower than the distribution for higher-mass objects (that likely formed via core fragmentation along with the stellar population, and formed within or became captured into binary orbits).Coupling a measurement of the mass and eccentricity of an object in a Keplerian orbit could therefore inform the interpretation of its formation history.If an object forms bottom up and has a high-eccentricity orbit, it appears to be an outlier; if the object forms like a "failed star," its eccentricity is less interpretable.A low mass/low mass ratio and low eccentricity could favor a "planetary" formation channel via either core accretion or disk instability, whereas a moderate mass and lower eccentricity could likewise favor a diskinstability formation, and a higher mass and/or eccentricity could favor a core fragmentation interpretation, i.e., the "failedbinary" interpretation.
Here, we present new interferometric measurements of the companion from the Very Large Telescope Interferometer (VLTI)/GRAVITY, including orbital coverage through periastron passage.GRAVITY has been used in the past few years to observe, for the first time with long-baseline optical interferometry, a handful of exoplanets and substellar companions: HR 8799 e (GRAVITY Collaboration et al. 2019), β Pic b and c, (Gravity Collaboration et al. 2020;Lagrange et al. 2020;Nowak et al. 2020;Lacour et al. 2021), HD 206893 B andc (Kammerer et al. 2021;Hinkley et al. 2023), PDS 70 b and c (Wang et al. 2021), HD 79246 B (Balmer et al. 2023), andHIP 65426 b (Blunt et al. 2023).For high-contrast targets, GRAVITY consistently achieves 50-100 μ as precision astrometry, while providing valuable spectral information about carbon-bearing molecules in the K band, which can be used to constrain atmospheric modeling (Gravity Collaboration et al. 2020;Kammerer et al. 2021;Blunt et al. 2023), in particular atmospheric retrievals (Mollière et al. 2020;Balmer et al. 2023).We use our new measurements of HD 136164 Ab along with the proper-motion measurements of the host star from the HGCA to derive a new dynamical mass estimate.We detect absorption due to carbon monoxide in the companion's K-band spectrum.The well-determined properties of this young brown dwarf companion make it an excellent prospect for future observations, and additional tests of evolutionary models.

VLTI/GRAVITY
We observed HD 136164 A and Ab four times between 2022 February and 2023 May using the GRAVITY instrument (Gravity Collaboration et al. 2017) at the European Southern Observatory (ESO)'s VLTI.We used the four 8.2 m unit telescopes (UTs) in dual-field on-axis fringe tracking mode (Lacour et al. 2019), where the fringe tracking fiber was placed at the location of the host star HD 136165 A and used to track the fringe pattern, and the science fiber was placed at the predicted location of the companion and integrated to detect the shifted fringe pattern from the companion.Observations made on 2022 February 19 and 2023 May 10 were taken as target visibility or bad weather backups to program 1104.C-0651 (PI: Lacour) in visitor or designated visitor mode (dVM).Observations made on 2022 April 21 and 2022 June 14 were taken as part of program 109.237J.001(PI: Balmer) in service mode (SM).Table 1 records our observing log.
During the first observation of the companion with GRAVITY on 2022 February 19, its position was relatively uncertain.The fiber placement was made based on a preliminary orbit fit to the two previous SPHERE observations (Wagner et al. 2020) taken 2 yr before the initial GRAVITY observation, and therefore the first GRAVITY observation suffered from fiber injection losses on the companion (with a theoretical coupling efficiency of γ = 0.88 < 1, as opposed to γ > 0.95 for subsequent observations).
We determined the complex visibilites of the host and the companion, which were phase referenced with the metrology system, using public release 1.5.0 (2021 July 133 ) of the ESO GRAVITY pipeline (Lapeyrere et al. 2014).We proceeded to decontaminate the science fiber flux due to the halo from the host star by simultaneously fitting the stellar contamination as a low-order polynomial, and the companion as a point source.

This pipeline is described in detail in Appendix A of Gravity Collaboration et al. (2020).
We determined the relative astrometry of the companion by analyzing the phase of the ratio of coherent fluxes.We calculate a 100 × 100 χ 2 periodogram power map over the fiber's field of view (Figure 1), and take the minimum of the χ 2 map as the preliminary companion position.We then calculate another 100 × 100 χ 2 grid with a range restricted to ±10 mas around the initial χ 2 grid minimum, and we estimate the uncertainty and covariance of the measurement by computing the rms of the χ 2 minima for each exposure.Because HD 136164 Ab is relatively bright we only observed three exposures per epoch.In order to estimate our errors, we divide each exposure into four subsections, each four detector integrations long, for a total of 12 χ 2 grid computations, except for 2023 May 10, where we only have two exposures with three detector integrations each, for six χ 2 grid computations.We find the typical precision from this technique is on the order of ∼100 μas, greater than 16.5 μas (the theoretical limit of VLTI/GRAVIY), due to systematic phase errors. 34Once the astrometry is determined, we extract the ratio of the coherent flux between the two sources at the estimated position of the companion, calculating the "contrast spectrum" of the companion.
Observations from 2022 February 19, 2022 April 21, and 2022 May 10 yield robust detections of the companion.Due to a metrological glitch during the fifth DIT in the observing sequence on 2022 June 14, the baselines for UT1 are unrecoverable for about half of the observing sequence, and we neglect these baselines in our extraction.We nevertheless recover the companion's signal in these data, albeit with larger astrometric uncertainties (of the same magnitude as the original SPHERE observations).Table 2 records our new GRAVITY relative astrometry of HD 136164 Ab around HD 136164 A. Figure 1 illustrates the detections of HD 136164 Ab within the fiber's field of view.
We transformed our contrast spectrum of the companion into a flux calibrated spectrum using a synthetic spectrum of the host star.We scaled a BT-NextGen (Allard et al. 2011) spectrum with T eff = 8100 K and et al. 2018) to archival photometry from Gaia (Gaia Collaboration et al. 2023), Tycho2 (Høg et al. 2000), and the Two Micron All Sky Survey (2MASS; Cutri et al. 2003) using species (Stolker et al. 2020).We broadened the spectrum with specutils (Astropy Collaboration et al. 2013Collaboration et al. , 2018Collaboration et al. , 2022) ) using a Gaussian kernel with an FWHM of 2 Å to match the rotational broadening of the star, and then sampled the broadened spectrum on the SPHERE and GRAVITY wavelength grids using spectres (Carnall 2017).The spectrum and stellar photometry are shown in Figure 2.  1 is presented chronologically, left to right.The dashed gray circle indicates the fiber's field of view.The origin is the placement of the science fiber on sky for a given observation, a prediction based on the previous available orbit fit, and the units are in ΔR.A. and Δdecl., e.g., the displacement of the fiber with respect to the host star HD 136164 A. The strongest peak in the χ 2 map (indicated with a blue arrow) reveals the position of the companion, with characteristic interferometric side lobes whose shape and distribution depend on the u-v plane coverage.Note.NEXP, NDIT, and DIT denote the number of exposures, the number of detector integrations per exposure, and the detector integration time, respectively, and τ 0 denotes the atmospheric coherence time.The fiber pointing is the placement of the science fiber relative to the fringe tracking fiber (which is placed on the central star).γ is the coupling efficiency at the position of the companion (see Table 2).

VLT/SPHERE
HD 136164 Ab was first detected using SPHERE (Beuzit et al. 2019) at ESO's VLT.We adopt the companion's 2015 June 19 and 2019 June 29 astrometry from Wagner et al. (2020) as is, although we note that, when compared with our updated orbital fits including our GRAVITY observations, the observation on 2019 June 29 appears systematically biased inwards in separation, likely because the off-axis point-spread function becomes asymmetric due to the transmission function of the coronagraph at close separations.Because our GRAVITY observations drive the orbital fits, and the 2015 June 19 astrometry is unaffected by this systematic error, we do not correct the SPHERE astrometry in this work.
We rereduced the IFS data from 2015 June 19 in order to extract the spectrum of the companion.The data were taken in the IRDIFS-EXT mode where the IFS (Claudi et al. 2008) and IRDIS (Dohlen et al. 2008) observe in parallel.The IFS covers wavelength bands Y, J, and H1 between 0.95 and 1.65 μm and the dual-band imager IRDIS covers K1 and K2 at 2.0 μm and 2.1 μm, respectively (Vigan et al. 2010).Because our GRAVITY data cover the K band, we do not rereduce the IRDIS data.The observations were obtained with the N_ALC_YJH_S coronagraph, which is optimized for close inner working angles.The raw data from the ESO archive were preprocessed with vlt-sphere, an open-source pipeline (Vigan 2020), and then we performed starlight subtraction using pyKLIP (Wang et al. 2015), an open-source implementation of Karhunen-Loève image processing (KLIP; Soummer et al. 2012).We used the forward modeling capabilities of KLIP (Pueyo 2016) to extract the spectrum of the companion using pyKLIP (Wang et al. 2016;Greenbaum et al. 2018), and applied the theoretical coronagraph transmission function to the resulting spectrum.The contrast spectrum was flux calibrated as in Section 2.1.

Orbital Analysis
In order to measure the orbit of HD 136164 A and b, we use orbitize!35(Blunt et al. 2020).We use orbitize!ʼsparallel-tempered (Vousden et al. 2016) affine-invariant (Foreman-Mackey et al. 2013) Markov Chain Monte Carlo (MCMC) algorithm.This algorithm fits for the six-parameter visual orbit (Green 1985), the system parallax, the proper motions and reference epoch astrometry, and the masses of the star and companion.We used 1000 walkers at 20 temperatures to sample our posterior distribution of orbits, with 15,000 steps of burn in discarded, before 5000 steps were recorded (for a total of 50,000 orbits recorded in our posterior distribution).We place a physically motivated normally distributed prior ( )   M 1.7, 0.1 (Bochanski et al. 2018;Kervella et al. 2022), and a log-uniform prior on the companion mass of  Note.We report the median and 68% confidence interval (CI) on each parameter derived from the posterior visualized in Figure 7. Propagating errors between both mass estimates yields q = 0.018 ± 0.005.a Next periastron passage after τ ref = 2020.0.
( )  M log 1200 J , a Gaussian prior on the system parallax based on the Gaia DR3 measurement ( )  8.202, 0.040 mas, and otherwise use default priors on all orbital elements as described in Blunt et al. (2020).We fit the model to relative astrometry from SPHERE (Wagner et al. 2020), GRAVITY (this work), and absolute astrometry of the host star from the HGCA (Brandt 2021).
Table 3 records the median and 1σ confidence intervals on our fit orbital parameters.Figure 3 plots the best-fit orbits to the observations of the system, and Figure 7 illustrates the posterior distribution of the best-fit orbits from our orbitize!MCMC.We note that the second SPHERE epoch deviates by a few milliarcseconds from our orbit fits, likely because the coronagraphic transmission at the location of the companion distorts the photocenter, biasing the astrometry in separation.We find a highly eccentric orbit, with e = 0.44 ± 0.03, and derive a companion mass of 35 ± 10 M J .Our derived stellar mass (1.87 ± 0.07 M e ) deviates by 1σ from our prior on the stellar mass, resulting in a mass ratio of q = 0.02.Our four epochs of relative astrometry at <100 μas precision, combined with the long-duration baseline information from the absolute astrometry, has strongly constrained the visual orbit of the companion.

Spectral Analysis
We compared the spectrum of the companion to two selfconsistent, cloudless, chemical equilibrium model grids, ATMO (Phillips et al. 2020) and SPHINX (Iyer et al. 2023).We chose these model grids because they use modern opacity sources and are computed for temperature ranges (2300-3000 K) appropriate for the spectral type of HD 136164 Ab.In particular, the SPHINX grid accounts for spectral broadening of key molecules, a range of metallicities and C/O ratios, and has been benchmarked to observations of M dwarfs with known bulk properties.However, it is computed for a lower resolution (R = 250) than our GRAVITY data (R = 500), so when fitting the grid to our data, we convolve the GRAVITY spectrum to the resolution of the grid using spectres (Carnall 2017).
We used the species36 package to fit our spectra to the model grid, linearly interpolating spectra between grid points.Using version 0.71 of species, we ingested the SPHINX grid from Iyer et al. (2022), which assumes a mixing length parameter of one.We initialized pyMultiNest37 (Feroz & Hobson 2008;Feroz et al. 2009;Buchner et al. 2014) via species to sample the interpolated grid with 500 live points.We measured the posterior distributions on the grid parameters, namely effective temperature (T eff ), log(g), radius, parallax, [Fe/H], and C/O for the SPHINX grid, and included a Gaussian process parameterized by a squared-exponential kernel to account for correlated noise between wavelength channels in the SPHERE data (see Section 4.1 in Wang et al.  2), plotted against the same sample of MCMC orbits as in the top left panel.The 1, 2, and 3σ confidence contours on each measurement are shown with decreasing transparency; the covariance of the measurement is related to the u-v plane coverage of the observation.2020).When fitting the GRAVITY data, species accounts for the correlation matrix of the spectrum in the fit.We weight each spectral point across the entire SED according to its wavelength spacing so that the GRAVITY spectrum, with an order of magnitude more data points, does not completely dominate the fit.To account for an apparent offset in the flux of the two spectra, we fit for a scaling parameter on the SPHERE YJH1 spectrum, as is common in the literature (e.g., Mollière et al. 2020).We place the same parallax prior on the atmospheric fits as we do the orbital fits in Section 3.1.
We fit each model grid to our data twice, first with a uniform prior on the mass (related to the sampled parameters ( ) g log and radius), and second with a Gaussian prior on the mass corresponding to the dynamical mass from our orbit fit in Section 3.1, ( )  M 35, 10 J .This experiment allows the sampler to determine the log(g) and radius based solely on the spectral and parallax measurements, and then adds the additional constraint of the dynamical mass.Figure 4 plots the spectrum of the companion and the results of the ATMO model fitting.We find the median sample has T eff , ( ) , and SPHERE scale factor of 2530 K, 3.5, 2.0 R J , −2.8, 0.70, respectively, with a uniform mass prior, and 2600 K, 4.3, 1.9 R J , −2.8, 0.71 with the dynamical mass prior.The residuals to the fit are strongly wavelength dependent, and appear related to regions where the water and iron hydride opacitics are dominant in shaping the spectral slope, especially when the dynamical mass prior is considered.This could indicate model deficiencies (nonsolar abundances, reduced temperature gradients, or disequilibrium chemistry) or data deficiencies (improperly corrected tellurics, phase errors, or coronagraphic effects)., and SPHERE scale factor of 2650 K, 4.6, 1.8 R J , 0.42, 0.47, 2.8, and 0.74, respectively, with the uniform mass prior, and 2640 K, 4.35, 1.9 R J , 0.39, 0.45, 2.8, and 0.74 with the dynamical mass prior.The residuals to the fit are smaller compared to the fit to the model assuming solar abundances (Figure 4), and the impact of the dynamical mass prior is more subtle.With more free parameters, the model can produce an equivalently good fit to the data with or without the dynamical mass prior, and no longer appears biased toward an implausibly low surface gravity without the dynamical mass prior.

Evolutionary Model Comparison
As a by-product of our spectral fits, we use species to integrate under the sampled ATMO and SPHINX spectra and determine the distributions of bolometric luminosity for both the uniform and dynamical mass prior cases.The ATMO model gives 2.8 0.3 for both cases, and the SPHINX model gives 2.8 0.3 for both cases.For substellar objects, there is a frequently observed tension between the evolutionary and spectral effective temperatures (see, e.g., Sanghi et al. 2023, their Figure 23).We leverage our dynamical mass to compare the bolometric luminosity expected for a 16 ± 2 Myr brown dwarf according to the ATMO2020 evolutionary model (Phillips et al. 2020) and the bolometric luminosity we derive from our spectral fits.The model uses the ATMO atmospheric model as a surface boundary condition to calculate the interior structure and evolution of substellar objects.We drew luminosities for a 35 M J brown dwarf at 14, 16, and 18 Myr, finding 2.68 evol 0.05 0.08 .Given the larger uncertainties on our spectrally derived bolometric luminosity (not to mention the potential systematic uncertainty on this estimate, given that we must rely on atmospheric models to extrapolate beyond the 1-2.5 μm range), the evolutionary and spectral bolometric luminosities for HD 136164 Ab agree to within 1σ.Increasing the wavelength coverage on the companion (especially at shorter wavelengths) will decrease the uncertainty on the spectrum-derived luminosity, and could reveal a tension if it exists.These models reveal the unphysical nature of the spectrally derived ( ) g log in the absence of a dynamical mass prior.
Figure 6(a) plots the distributions of mass and radius from our atmospheric fits, compared to the ATMO2020 mass-radius model (Phillips et al. 2020) for an age of 16 ± 2 Myr, showing that these parameters are in agreement within 1-2σ with expectations from the evolutionary model.Figure 6(b) plots the distributions of mass and effective temperature.When the dynamical mass is included as a prior in the spectral fit, the Figure 5.The spectrum of HD 136164 Ab (black data points) fit to the cloudless SPHINX equilibrium chemistry self-consistent model, which varies abundance in addition to effective temperature and surface gravity, with and without a prior on the dynamical mass of the companion (median samples are plotted as blue and red curves, respectively).30 random samples from each posterior distribution are plotted as semitransparent curves.In each fit, the SPHERE YJH1 data are allowed to scale by a factor, so for the fit with the dynamical mass prior, the scaled SPHERE data are shown as black circles, and without the dynamical mass prior, they are shifted to the right for visualization as gray squares.
Figure 6.Left: mass vs. effective temperature for our atmospheric fits (contours) compared to the mass-radius relationship from the ATMO2020 evolutionary model (black curve).The uncertainty on the system's age is represented by gray contours to the black curve, 16 ± 2 Myr.When the dynamical mass is included as a prior in the atmospheric model fit, the mass and effective temperature are consistent with the evolutionary model prediction to within 1σ.Right: mass vs. radius for our atmospheric fits compared to the mass-radius relationship from the evolutionary model.Without a dynamical mass prior, our spectral fits derive a ( ) g log that results in a mass inconsistent with expectations from the evolutionary model.With a dynamical mass prior, the determined ( ) g log is consistent with expectations to within 1σ.
temperatures are consistent within 1σ.For many brown dwarfs, there is significant tension between the evolutionary and spectral effective temperatures.

The Formation of HD 136164 Ab
We have demonstrated that HD 136164 Ab has an eccentric orbit (e = 0.44 ± 0.03), a dynamical mass of 35 ± 10 M J , and a mass ratio with HD 136164 A of q = 0.018 ± 0.005 (see Section 3.1).This makes HD 136164 A the youngest substellar companion with a dynamical mass estimate.The dynamical mass agrees within the uncertainties with the previous evolutionary-model-derived mass based on SPHERE spectrophotometry, which are in perfect agreement when considering the evolutionary models anchored to the SPHERE/IRDIS Kband photometry (Wagner et al. 2020).The dynamical mass and mass ratio is inconsistent with core-accretion models of formation (e.g., Mordasini et al. 2009;Emsenhuber et al. 2021Emsenhuber et al. , 2021c)).The good agreement with the dynamical mass, spectral effective temperature, and the ATMO isochrone with an Sco-Cen age (see Figure 6) can be interpreted as evidence for the similarity in age of the host and companion, the coevality of the two, and therefore additional evidence against core accretion.
The simulations of disk fragmentation formation and fragment-fragment interactions in Forgan & Rice (2013) and Forgan et al. (2015Forgan et al. ( , 2018) ) appear to indicate the result of diskinstability formation is a distribution of single giant planets/ brown dwarfs near the deuterium burning limit at wide separations (10-100 au), with a distribution peaking at low orbital eccentricities, but flattening toward a nonzero relative frequency between e = 0.1 and 0.8.This could tentatively support the interpretation that HD 136164 Ab formed via disk instability, since it is a widely separated companion at 22.5 au.It could be that the companion formed from disk fragmentation with an initially low eccentricity that was then increased via dynamical interactions with the outer M dwarf HD 136164 B at 650 au, but this would take time to evolve dynamically, and the system is young.Demonstrating this would require dynamical modeling and additional orbital monitoring of B and Ab.We could interpret the slightly enhanced metallicity we derive for the companion using the SPHINX models as evidence the companion accreted a significant fraction of metals from the circumstellar environment, if it indeed formed via disk instability.For instance, it could have undergone late accretion of solids (e.g., Nowak et al. 2020).The posterior distribution of abundances from our model fits to the SPHINX grid are wide, and still encompass solar metallicity solutions at the 2-3σ level.In Wagner et al. (2020), it was suggested that a low C/O ratio could imply an earlier period of gas accretion from the disk, because during the Class 0 stage of star/circumstellar evolution, no CO gas would be frozen out, whereas during the Class 1 stage, CO gas would be frozen out beyond 30 au.Our best-fit SPHINX spectrum has a lower C/O ratio of 0.45, although the posterior still spans a wide range of C/O encompassing the solar value at 2-3σ, and similarly this could be loosely interpreted as evidence for formation of the companion within a circumstellar disk.The host star could also have a subsolar C/O abundance, and the lower C/O abundance for HD 136164 Ab could be reflective of a match to the stellar abundance.Kratter et al. (2010) show that disk fragmentation can produce brown dwarf and even planetary mass companions, but provide the additional interpretation that these systems are necessarily "failed stars," i.e., that the distribution of diskinstability-born companions has a low-mass tail, but increases above the deuterium burning limit.They predict an increase in the occurrence of brown dwarf companions to A-type stars; as is the case here, but also appears to be the case more generally (Nielsen et al. 2019).
If it did not form from a fragmenting circumstellar disk, could HD 136164 Ab have formed via the collapse of a molecular cloud, as a product of the IMF?With opacity-limited minimum mass estimates as low as 1-3 M J for core fragmentation (e.g., Boyd & Whitworth 2005;Bate 2009), or a derived minimum mass of around 6-18 M J from radiation hydrodynamical simulations (e.g., Bate 2012), the mass of HD 136164 Ab is not unprecedented as the by-product of a fragmenting core mass function.In Bate (2012), the outcome of the stellar IMF produces a handful of binaries with small mass ratios (0.01-0.02) with a relatively uniform distribution of eccentricities.HD 136164 Ab's high eccentricity and mass ratio of 0.02 could therefore indicate a core fragmentation formation pathway, if eccentricities are expected to be damped during formation in a disk, but not during core fragmentation.In the context of population-level studies of substellar companions, this eccentricity also empirically indicates more of a "brown dwarf/failed star" nature, as directly imaged planets appear to follow a distinct eccentricity distribution with generally lower eccentricities (as their eccentricities are expected to be damped by their parent disks), whereas substellar objects produced by fragmenting cores tend to have a distribution of on average higher eccentricities (Bowler et al. 2020;Doó et al. 2023;Nagpal et al. 2023).Parker & Daffern-Powell (2022) use Nbody simulations of star-forming regions to show that earlytype stars can capture single brown dwarfs or even steal planets from other lower-mass systems within the first 10 Myr of an association's lifetime, so it is a reasonable assumption that, even if the brown dwarf did not form directly into a binary configuration with HD 136164 A, it could have formed nearby and been subsequently captured onto a binary orbit.
It is still unclear whether this brown dwarf companion could have formed from disk fragmentation or core fragmentation, but it appears clear that the companion is a "failed star" coeval with the host, either born or captured into a binary orbit.

Atmospheric Modeling in the Context of a Dynamical Mass Prior
In Balmer et al. (2023), we saw that the atmospheric modeling of a close separation brown dwarf companion can be strongly influenced by the inclusion of a dynamical mass prior.In that work, for a benchmark brown dwarf with a known mass and stellar abundances, the ( ) g log and radius of the atmospheric model appeared to be strongly influenced by regions of increased systematics in the spectrum.It is necessarily the case that spectral and astrometric measurements made at such close separations will suffer, at least initially, from such systematic errors.These could be responsible for the wavelengthdependent residuals between our data and the ATMO models.It is also the case that model deficiencies in solar abundance, chemical equilibrium, cloudless, radiative-convective equilibrium models (e.g., nonsolar abundances, clouds or reduced temperature gradients, vertical mixing and disequilibrium chemistry, etc.) could impact the determined physical parameters and quality of fit.The SPHINX grid, which varies abundances, can provide a better fit to the data, invoking a subsolar C/O ratio and an enhanced metallicity.In Figure 4, the residuals to the ATMO model occur near the edges of the instruments' respective wavelength range, and appear to influence the fitted ( ) g log strongly.In our fit without the dynamical mass prior, the posterior distribution was driven to unphysically low values, ( ) = g log 3.0.When the dynamical mass prior is considered, we find a much more reasonable ( ) = g log 4.4, and a corresponding radius of 1.9 R J .It could be that a circum-secondary disk or dust extinction produces the slope in the GRAVITY data.We briefly considered including an extra blackbody term, or an extinction term, in our spectral fitting to account for these effects, and were unable to constrain its presence with our data convincingly, so we do not consider that here (the disk excess could be investigated with 3-5 μm JWST observations).
Given the good agreement between our atmospherically derived luminosity (based on our fits to the spectrum using a dynamical mass prior) and our evolutionarily derived luminosity (based on the system age and dynamical mass), for relatively cloudless late M-type substellar companions without a dynamical mass prior, it appears reasonable to adopt an evolutionary-model-derived mass prior based on the observed flux of the companion when fitting the spectrum with model grids.Doing so will ensure the fit is not driven to unphysical values of ( ) g log .We note, however, that if a fit is driven to unphysical values, it could be due to a model deficiency since, in the case of HD 136164 Ab, we found a model grid that varied abundances away from solar values derived physical ( ) g log values without the dynamical mass prior.As in Balmer et al. (2023), the dynamical mass prior appears most useful in determining whether there are data or model deficiencies, and what changes (to data treatment or model construction) might alleviate the tension between the two.

Conclusions
In this work we: 1. Observed the brown dwarf companion HD 136164 Ab four times with VLTI/GRAVITY dual-field interferometry.We detected the companion in all four observations despite limited field rotation and exposure time, demonstrating the power of GRAVITY to characterize brown dwarf companions efficiently.2. Extracted precise astrometry of the companion (median σ = 40 μas) and a contrast spectrum at R = 500 resolution from the interferometric observables.We also reextracted the SPHERE IFS YJH1 spectrum of the companion and flux calibrated these contrast spectra using a stellar template.3. Fit the orbit of the companion with a combination of relative astrometry from GRAVITY and SPHERE, and absolute astrometry from the HGCA, determining the eccentricity and dynamical mass of the companion for the first time.We leverage the moderate eccentricity (e = 0.44 ± 0.03) and mass ratio (q = 0.018 ± 0.005) to assess the formation of the object.4. Fit the spectrum of the companion to two grids of selfconsistent, cloudless atmospheric models, ATMO and SPHINX.We discussed the inclusion of the dynamical mass as a prior in atmospheric fits for this object.The SPHINX grid, which varies abundances, indicates a slightly subsolar C/O ratio (0.45) and enriched metallicity (0.4), with uncertainties that are consistent with solar values.Deriving the luminosity from our spectral fits and comparing them to expectations from evolutionary models revealed an apparent agreement, albeit with uncertainties of 0.3 dex in luminosity. 5. Discussed the formation of HD 136164 Ab in the context of the new GRAVITY data and especially the dynamical mass.We rule out formation via core accretion, but present evidence that could be interpreted in favor of either disk fragmentation or cloud fragmentation.
HD 136164 Ab now joins a very select group of young substellar companions with dynamical mass measurements (see Figure 4 in Franson & Bowler 2023).It is currently the youngest substellar companion with a dynamical mass estimate.Like PZ Tel B (Biller et al. 2010;Mugrauer et al. 2010;Schmidt et al. 2014;Maire et al. 2016;Stolker et al. 2020;Franson & Bowler 2023), this object will be crucial for benchmarking evolutionary models in the coming decade.HD 136164 Ab does not need to stand alone, however.We have shown that, even without stellar radial velocities, coupling GRAVITY observations with the HGCA can produce wellconstrained dynamical mass estimates for young companions orbiting stars that are not amenable to radial velocity characterization.Future work should look to estimate the dynamical mass of the young substellar companions to HD 115470 (HIP 64892 B; Cheetham et al. 2018), μ 2 Sco (μ 2 Sco b; Squicciarini et al. 2022), HD 149274 (HIP 81208 B;Chomez et al. 2023;Viswanath et al. 2023), and the planet orbiting HD 116434 (HIP 65426b; Chauvin et al. 2017;Cheetham et al. 2018;Stolker et al. 2020;Petrus et al. 2021;Blunt et al. 2023).While these might not result in dynamical masses of same precision as this result (since the aforementioned companions are on longer-period orbits), a similar analysis might at the very least provide useful mass upper limits for assessing the formation histories of these objects.
For HD 136164 Ab itself, there remains a number of additional investigations that can continue to refine estimates of its bulk and atmospheric properties, to aid in the benchmarking of various models.Future observations of the companion at shorter wavelengths (e.g., z′ or i′ band with MagAO-X; Males et al. 2022), and longer wavelengths (e.g., with JWST/NIRCam at 3-5 μm, via GO program 1902, PI: Kammerer) will produce a precise bolometric luminosity estimate and could probe the existence of a circum-secondary disk with more certainty than our observations.High-resolution spectroscopy with HiRISE (Vigan et al. 2023) could assess whether the companion is truly enhanced in metallicity.Epoch astrometry from Gaia DR4 should improve the precision of the companion's dynamical mass estimate further, and enable comprehensive benchmarking of multiple evolutionary models at very young ages.3. Orbital parameters with respect to HD 136164 Ab are denoted with subscript "1," and parameters with respect to HD 136164 A with subscript "0."The distribution of orbital elements is unimodal, with no obvious degeneracies present between parameters, thanks to the combination of precise relative astrometry from GRAVITY, and long-duration baseline absolute astrometry from the HGCA considered in the fit.
Figure 8.The posterior distributions of atmospheric model parameters from our ATMO spectral fit, presented in Section 3.2.The contours of the posteriors from the fit with a uniform mass prior are plotted in red, and the contours of the posteriors from the fit with a dynamical mass prior are plotted in blue.The major takeaway is that ( ) g log and mass remain relatively unconstrained in the uniform prior fit, but the other parameters remain effectively the same between the two.
Figure 9.The posterior distributions of atmospheric parameters from spectral fit, shown in Section 3.2.The contours of the posteriors from the fit with a uniform mass prior are plotted in red, and contours the posteriors from the with a dynamical mass prior are plotted in blue.The major takeaway is that, with a slightly enhanced metallicity and subsolar C/O, the sampler is no longer driven toward implausibly low ( ) g log .

Figure 1 .
Figure 1.Detections of HD 136164 Ab with VLTI/GRAVITY.Each panel visualizes the χ 2 map calculated after subtracting the stellar contamination.Each epoch in Table1is presented chronologically, left to right.The dashed gray circle indicates the fiber's field of view.The origin is the placement of the science fiber on sky for a given observation, a prediction based on the previous available orbit fit, and the units are in ΔR.A. and Δdecl., e.g., the displacement of the fiber with respect to the host star HD 136164 A. The strongest peak in the χ 2 map (indicated with a blue arrow) reveals the position of the companion, with characteristic interferometric side lobes whose shape and distribution depend on the u-v plane coverage.

Figure 2 .
Figure 2. The spectral energy distribution (SED) of HD 136164 A. Archival photometry is overplotted a BT-NextGen model spectrum with T eff = 8100 K and ( ) = g log 4.2.This spectrum is used to transform the companion's contrast spectra into flux calibrated spectra.

Figure 3 .
Figure 3. HD 136164 orbital analysis.Top left: the visual orbit of HD 136164 Ab around HD 136164 A (black star).1000 randomly chosen orbits from our fit using the MCMC are shown in color, with the color corresponding to the eccentricity of a given orbit.Top middle and right: the proper motion of the host star HD 136164 A over time in R.A. and decl.(middle and right, respectively).The proper motion of the host as predicted by the orbit fit is shown in curves colored corresponding to the mass of the companion.Observations from Hipparcos and Gaia are shown with error bars as recorded in the HGCA.Note that the orbitize!proper-motion loglikelihood function evaluates the time-averaged proper motions, and not the instantaneous proper motions that are overplotted here (see Hinkley et al. 2023); the instantaneous measurements are shown here for visualization purposes only.Bottom: VLTI/GRAVITY measurements of the position of HD 136164 Ab across four epochs (see Table2), plotted against the same sample of MCMC orbits as in the top left panel.The 1, 2, and 3σ confidence contours on each measurement are shown with decreasing transparency; the covariance of the measurement is related to the u-v plane coverage of the observation.

Figure 5
Figure 5 plots the spectrum of the companion and the results of the SPHINX model fitting.We find the median sample has T eff , ( ) g log , R p , [Fe/H], C/O, ( )  L L log, and SPHERE scale factor of 2650 K, 4.6, 1.8 R J , 0.42, 0.47, 2.8, and 0.74, respectively, with the uniform mass prior, and 2640 K, 4.35, 1.9 R J , 0.39, 0.45, 2.8, and 0.74 with the dynamical mass prior.The residuals to the fit are smaller compared to the fit to the model assuming solar abundances (Figure4), and the impact of the dynamical mass prior is more subtle.With more free parameters, the model can produce an equivalently good fit to the data with or without the dynamical mass prior, and no longer appears biased toward an implausibly low surface gravity without the dynamical mass prior.

Figure 4 .
Figure4.The spectrum of HD 136164 Ab (black data points) fit to the cloudless ATMO equilibrium chemistry self-consistent model, with and without a prior on the dynamical mass of the companion (median samples are plotted as blue and red curves, respectively).30 random samples from each posterior distribution are plotted as semitransparent curves.In each fit, the SPHERE YJH1 data are allowed to scale by a factor, so for the fit with a dynamical mass prior, the scaled SPHERE data are shown as black circles, and without a dynamical mass prior, they are shifted to the right for visualization as gray squares.The SPHERE and GRAVITY spectra are available as the data behind the figure.(Thedata used to create this figure are available.)

Figure 7 .
Figure7.Posterior distributions of orbital elements in the orbital analysis presented in Section 3.1, and visualized in Figure3.The median and 1σ CIs from the marginalized 1D histograms show here, for each parameter, are recorded in Table3.Orbital parameters with respect to HD 136164 Ab are denoted with subscript "1," and parameters with respect to HD 136164 A with subscript "0."The distribution of orbital elements is unimodal, with no obvious degeneracies present between parameters, thanks to the combination of precise relative astrometry from GRAVITY, and long-duration baseline absolute astrometry from the HGCA considered in the fit.

Table 2
New Relative Astrometry of HD 136264 Ab around HD 136264 A ΔR.A. × σ Δdecl.on the off-diagonal.
Note.The covariance matrix can be reconstructed using s DR.A. 2 and s Ddecl. 2 on the diagonal, and ρ × σ

Table 3
Orbital Parameters Inferred for HD 136164 Ab in This Work