Direct-imaging Discovery of a Substellar Companion Orbiting the Accelerating Variable Star HIP 39017

We present the direct-imaging discovery of a substellar companion (a massive planet or low-mass brown dwarf) to the young, γ Doradus (γ Dor)-type variable star HIP 39017 (HD 65526). The companion’s SCExAO/CHARIS JHK (1.1–2.4 μm) spectrum and Keck/NIRC2 L′ photometry indicate that it is an L/T transition object. A comparison of the JHK+L ′ spectrum to several atmospheric model grids finds a significantly better fit to cloudy models than cloudless models. Orbit modeling with relative astrometry and precision stellar astrometry from Hipparcos and Gaia yields a semimajor axis of 23.8−6.1+8.7 au, a dynamical companion mass of 30−12+31 M J, and a mass ratio of ∼1.9%, properties most consistent with low-mass brown dwarfs. However, its mass estimated from luminosity models is a lower ∼13.8 M J due to an estimated young age (≲115 Myr); using a weighted posterior distribution informed by conservative mass constraints from luminosity evolutionary models yields a lower dynamical mass of 23.6−7.4+9.1 M J and a mass ratio of ∼1.4%. Analysis of the host star’s multifrequency γ Dor-type pulsations, astrometric monitoring of HIP 39017 b, and Gaia Data Release 4 astrometry of the star will clarify the system age and better constrain the mass and orbit of the companion. This discovery further reinforces the improved efficiency of targeted direct-imaging campaigns informed by long-baseline, precision stellar astrometry.


INTRODUCTION
Large, ground-based telescopes assisted by facilityand, now, specialized extreme -adaptive optics (AO) systems have allowed the imaging of substellar and self-luminous Jovian companions to nearby stars (e.g.Marois et al. 2008;Lagrange et al. 2010;Macintosh et al. 2015;Chauvin et al. 2018;Currie et al. 2023b).Direct imaging is uniquely poised to probe models of planet formation; not only is it sensitive to wide-separation companions in young systems, it also allows both spectral and dynamical orbital analysis for any detected companion.However, large-scale, "blind" (or unbiased) surveys have found that planetary and brown dwarf companions that can be imaged with these modern, ground-based systems are rare, with only about a few percent of targeted systems yielding companion detections (Nielsen et al. 2019;Vigan et al. 2021;Currie et al. 2023a).
To account for this shortcoming, recent direct imaging surveys have turned toward targeted searches, focusing on stars showing dynamical evidence for a companion.In particular, the 25 yr precision astrometry Hipparcos-Gaia Catalog of Accelerations (HGCA; Brandt 2018Brandt , 2021) ) identifies stars displaying acceleration inconsistent with a single-body solution.Such an acceleration could indicate a massive unseen companion, making young, nearby accelerating systems prime targets for imaging.Surveys of accelerating stars have yielded a significantly higher detection rate (Currie et al. 2021;Bonavita et al. 2022) for low-mass stellar (Steiger et al. 2021;Chilcote et al. 2021), brown dwarf (Currie et al. 2020a;Bonavita et al. 2022;Kuzuhara et al. 2022;Franson et al. 2023a;Li et al. 2023), and even planetary companions of young stars (e.g.HIP 99770 b and AF Lep b; Currie et al. 2023b;De Rosa et al. 2023;Mesa et al. 2023;Franson et al. 2023b).
Given the youth of the substellar and planetary companions accessible with direct imaging, they can serve as laboratories for testing and constraining models of their formation and atmospheres (Spiegel & Burrows 2012).Dynamical mass measurements of imaged companions can provide an independent point of comparison for spectral models.However, due to the large orbits of companions accessible with imaging, the precision of the dynamical masses derived from relative astrometry is often severely limited by small orbital coverage.Selecting imaging targets with existing long-baseline precision astrometric measurements from the HGCA offers the ad-ditional benefit of extending orbital coverage for one of the components of the system (e.g.Brandt et al. 2021a).
Age is another crucial factor in constraining the fidelity of atmospheric models for young substellar and planetary companions.These bodies cool with time, with brown dwarfs manifesting this as an apparent evolution to later spectral types.Uncertainties in system ages can have a significant impact on masses retrieved from atmospheric and evolutionary models (e.g.Moro-Martín et al. 2010).However, if the system is not a member of a coeval stellar association, obtaining precise and reliable independent age constraints can be difficult.While a number of methods to constrain stellar ages from models and empirical relations exist, they may not be applicable to every system (e.g.Dantona & Mazzitelli 1984;Mamajek & Hillenbrand 2008).
In this work, we report the direct imaging discovery of a massive planet/very-low mass brown dwarf companion ∼22 au from the γ Doradus (γ Dor)type variable accelerating star, HIP 39017, with spectrum resembling that of an L/T transition object.An overview of the host star is provided in Section 2.1.Observations and companion detection are described in Sections 2.2 and 2.3.Our spectral and astrometric analyses are presented in Sections 3 and 4, respectively.In Section 4.2.3, we revisit the astrometric analysis with posteriors informed by evolutionary models and the HIP 39017 system's likely age.Finally, further discussion is included in section 5.
While Tetzlaff et al. (2011) estimate an extremely young age of ∼ 11 Myr through comparison to several stellar evolutionary models, they assume an A3 spectral type, which is in conflict with the Henry et al. results.A comparison with the Banyan-Σ tool (Gagné et al. 2018) does not identify HIP 39017 as a member of any moving group or young association.There is currently no evidence that HIP 39017 harbors a luminous debris disk from Wide-field Infrared Survey Explorer, Spitzer, or Herschel results, which are common for early-type stars at such young ages (e.g.Rieke et al. 2005) et al. 2008) than the Pleiades and Hyades, this difference in metallicity is not sufficient to shift an object of Hyades age or older to the position of HIP 39017 on the CMD.Therefore, we adopt ∼ 115 Myr as an approximate age of the system, with a conservative upper limit of ∼ 700 Myr for analysis that will later depend on these values (see Section 4.2.3).1 Despite the star's youth and proximity to the Sun, it was not targeted as a part of major, re-cent high-contrast imaging surveys conducted with the Spectro-Polarimetric High-contrast Exoplanet REsearch or Gemini Planet Imager instruments (Beuzit et al. 2019;Macintosh et al. 2014) and has not been observed with any other AO instrument on any 8-10 m telescope nor with the Hubble Space Telescope.However, the early Data Release 3 (eDR3) version of the HGCA (Brandt 2021) reveals that HIP 39017 has a statistically significant astrometric acceleration (5.5σ), plausibly due to an unseen substellar companion.HIP 39017 is bright enough (R ∼ 6.8 mag; Monet et al. 2002) to achieve high performance of ground-based extreme-AO.Thus, we targeted it as a part of the SCExAO high-contrast imaging survey of nearby, young accelerating stars (Currie et al. 2021(Currie et al. , 2023b)).

High-Contrast Imaging Observations
We conducted AO-assisted, high-contrast imaging observations of HIP 39017 using the Subaru and Keck II telescopes on Maunakea on five nights between February 2022 and March 2023 (Table 1).For the Subaru observations, the facility AO system, AO188 (Hayano et al. 2008), provided a first stage correction of atmospheric turbulence, and SCExAO then delivered a higher order correction of residual wavefront errors (Jovanovic et al. 2015b;Currie et al. 2020b); Keck II observations used either the near-IR Pyramid wavefront sensor (PyWFS Bond et al. 2020) (March 2022) or the facility Shack-Hartmann wavefront sensor (other data sets).The Subaru and Keck observations utilized a Lyot or vector vortex coronagraph, respectively, to suppress residual stellar halo light.For Subaru, sharpened starlight was then sent to the CHARIS integral field spectrograph (IFS; Groff et al. 2016), while the NIRC2 camera recorded Keck images.
All data were obtained in pupil tracking/angular differential imaging (ADI) mode (Marois et al. 2006).CHARIS data were obtained in low resolving power (R ∼ 18) spectroscopic mode to obtain wide wavelength coverage (1.1-2.4 µm), while NIRC2 data were taken in the L ′ -band broadband filter.Satellite spots produced by a modulation of SCExAO's deformable mirror provided spectrophotometric calibration for our CHARIS data and were used to register the star to a common center (Jovanovic et al. 2015a;Sahoo et al. 2020).For NIRC2, we calibrated the photometric measurements by obtaining unsaturated images of the central star before and after the main coronagraphic sequence.For the March 2022 CHARIS data, we offsetted the star from the image plane center by 0.  wide-separation candidate companion seen in the first Keck/NIRC2 epoch that is outside the nominal CHARIS field of view and was shown to be a background star (see Appendix).
We used data cubes extracted from the raw CHARIS data by the Automated Data Extraction, Processing, and Tracking System for CHARIS (ADEPTS; Tobin et al. 2020Tobin et al. , 2022)), which utilizes the CHARIS Data Reduction Pipeline (DRP; Brandt et al. 2017).We performed subsequent processing steps on the extracted CHARIS cubes -sky subtraction, image registration, spectrophotometric calibration, point-spread function subtraction, and spectral extraction -using the CHARIS Data Processing Pipeline (Currie et al. 2020b).For spectrophotometric calibration, we assumed an F0V model atmosphere from the Kurucz library (Castelli &  The integration time of each exposure, the number of exposures used in our analysis, and the total variation of parallactic angle in each sequence are represented by texp, Nexp, and ∆P A, respectively.a These CHARIS data had the star offsetted by 0. ′′ 5 to bring the background star within the field of view b The CFHT seeing monitor did not record values, but the on-sky PSF quality was consistent with median to slightly better-than-median seeing. Kurucz 2003) with the star's Two Micron All Sky Survey (2MASS) photometry (Skrutskie et al. 2006(Skrutskie et al. , 2019)).For the February 2022 and December 2022 data, we subtracted the PSF using A-LOCI (Currie et al. 2012(Currie et al. , 2015) ) in combination with ADI only.For the shallower March 2022 data focused on an outer candidate, we simply subtracted the radial intensity profile of each slice of each cube before rotating all cubes to the north and mediancombining.
We processed the NIRC2 L ′ images with the VIP package from Gomez Gonzalez et al. (2017) and the H-band data using the general-use pipeline from Currie et al. (2011Currie et al. ( , 2014)).For pre-processing, we correct for pixel sensitivity variations through flat-fielding and dark subtraction in our science images.Cosmic ray artifacts are removed using the lacosmic Python package (van Dokkum 2001), and geometric distortion is addressed combining solutions from Service et al. (2016) and the script from Yelda et al. (2010) for the NIRC2 camera's narrow-field mode.The post-alignment data cube has a centering stability of less than 0.5 mas, which we also take into account in our final astrometry.For L ′ -band data post-processing, we applied the annular PCA algorithm in VIP with six principle components to subtract stellar PSF, following the methods outlined in Li et al. (2023).For the H-band data, we applied a radial profile subtraction to each image before combining north-rotated images.

Detection of HIP 39017 b
Figure 2 shows the HIP 39017 images obtained from two CHARIS data sets (February 2022 and December 2022) and one NIRC2 data set (March 2022).Each data set reveals a faint companion about ≈0.′′ 4 southwest from the star.The companion, hereafter HIP 39017 b, is detected with signal-to-noise ratios (SNRs) of 6-242 .The December 2022 CHARIS detection is the highest fidelity one, as HIP 39017 b is visible in each channel.The signal-to-noise and relative astrometry for the companion in each epoch are listed in Table 2.
To correct our HIP 39017 b SCExAO/CHARIS spectrophotometry and astrometry for biasing due to processing, we carried out forward modeling as in Currie et al. (2018).HIP 39017 b's throughput ranged from 27% to 80%, depending on the data set and channel.Our forward-modeling analysis indicates minimal astrometric biasing (∼ 0.05 pixels in both x and y).To estimate astrometric errors, we imputed the planet forwardmodel at the same approximate angular separation but 1000 different azimuthal angles, compared the imputed and recovered astrometry, and computed the scatter in these values (e.g.Kuzuhara et al. 2022).Based on these analyses, we estimate an intrinsic centroiding uncertainty of ∼0.1 pixels in both coordinates, which we add in quadrature with an assumed stellar centroiding uncertainty of 0.25 pixels and converted to values in polar coordinates.
For Keck/NIRC2 L ′ data, we opted for a simpler method of inserting negative copies of HIP 39017 b over a range of separations and with a range of brightnesses until the companion's signal is fully nulled.While the data.Because of its slightly higher SNR and lower levels of residual speckle contamination in J band, we focus our analysis on the December spectrum.
SCExAO/CHARIS detections have a higher SNR, the instrument's centroiding uncertainty coupled with the larger north position angle offset result in an uncertainty comparable to NIRC2's in the companion's position angle.

INFRARED SPECTRUM AND PHOTOMETRY
Figure 3 shows the CHARIS spectrum for HIP 39017 b extracted from both epochs.The February 2022 spectrum is much noisier in J band due to greater residual speckle noise contamination.However, the spectra agree within errors in 18 of the 22 channels, including all but one channel in the H and K passbands: the February data appear to have a systematic offset of ∼ 5-10%, consistent with slight variations in the absolute spectrophotometric calibration.We estimate HIP 39017 b's broadband photometry in the major Mauna Kea Observatories passbands from December 2022 CHARIS data to be J = 18.31 ± 0.20, H = 17.75 ± 0.13, and K s = 17.02 ± 0.14, consistent with a substellar object at the L/T transition (e.g.Dupuy & Liu 2012;Currie et al. 2023b) 3 .
We derive a companion-to-host contrast of ∆L ′ = 9.93 ± 0.07 magnitudes for the NIRC2 L ′ band data.Using the intrinsic colors for main sequence stars from Kenyon & Hartmann (1995) and Pecaut & Mamajek (2013), we adopt a primary star brightness of m L ′ ,⋆ = 6.19 ± 0.044 mag, the same as its W1 magnitude (Cutri et al. 2012), as the two are equivalent in the Rayleigh-3 Photometry derived from the February 2022 data formally agrees within errors: J = 18.91 ± 0.41, H = 17.81 ± 0.13, and Ks = 17.12 ± 0.13.
Jeans tail of A stars.HIP 39017 b's brightness in L ′ is then 16.12 ± 0.18 mag. Figure 4 compares HIP 39017 b's broadband photometry to measurements for directly-imaged planets, and directly-imaged very low-mass brown dwarfs, and field brown dwarfs.The companion's position in both J/J-K and L ′ /H-L ′ CMDs overlaps with objects at the L/T transition (Saumon & Marley 2008).HIP 39017 b's colors also show a strong resemblance to those for directlyimaged companions around accelerating stars HIP 99770 b, HD 33632 Ab, and especially HIP 21152 B.

Comparison to Empirical Spectra
Next, we compared the CHARIS spectrum from the December 2022 data set to the archive of published L and T dwarfs from the SpeX Prism Spectral Libraries (Burgasser 2014).
The template spectra were first resampled to linearlyspaced bins using the FluxConservingResampler in Astropy's specutils package (Astropy-Specutils Development Team 2019; Carnall 2017), before being smoothed and resampled down to the same R as the CHARIS spectrum.Each template spectrum was then fit by minimizing where F ν,k and f ν,k are the flux of the template and target spectra, respectively, in the k th bin, and α is the scaling factor.The covariance matrix, C, includes the uncertainties from both the template and target spectra.The covariance for the target SCExAO/CHARIS spectrum is derived and parameterized following Greco & Brandt (2016), with the December spectrum best fit by A ρ = 0.0675, A λ = 0.0525, A δ = 0.88, σ ρ = 0.8175, and σ λ = 2.8.In each fitting iteration, the uncertainty of the template spectrum, scaled by the amplitude α for that iteration, is added in quadrature to the diagonal of the covariance matrix for the target (Currie et al. 2018).
Figure 5 shows the reduced-χ 2 achieved for each fitted spectrum as a function of its near-infrared (NIR) spectral type, and Figure 6 shows the four best-fit spectra compared to that of HIP 39017b.Within 1σ, the SpeX Prism sources that best fit the HIP 39017b spectrum indicate that it is ∼T1-2, though a closer look at those best-fit sources may be more accurately described as a peculiar T1 − 2 dwarf.
The best-fit source to the HIP 39017 b JHK SCExAO/CHARIS spectrum is 2MASS J13243553+6358281 (Looper et al. 2007).Classified at the time as a peculiar T2 dwarf, spectral matching analyses postulated that its unusually red colors Note-S/N represents the object's signal-to-noise ratio, ρ represents the companion.2021).Indigo error bars with filled circles are the HR 8799 planets (Marois et al. 2008;Metchev et al. 2009;Skemer et al. 2014;Currie et al. 2014;Zurlo et al. 2016;Skrutskie et al. 2006Skrutskie et al. , 2019)).ν values that correspond to the 1, 2, and 3 σ confidence intervals, from bottom to top.(Metchev et al. 2008) could be a result of an unresolved binary (Looper et al. 2007;Burgasser et al. 2010), though a more recent identification of this source as a member of the ∼ 150 Myr old AB Doradus group suggests that its redness may instead be ascribed to its young age (Gagné et al. 2018).HIP 39017's young age is certainly consistent with the appearance of a redder spectrum for HIP 39017b, like that seen in 2MASS J13243553+6358281.
The only other SpeX Prism library source with a χ 2 ν fit within 1σ is SDSS J103931.35+325625.5 (2MASS J10393137+3256263; Burgasser et al. 2010).While its NIR spectrum was closest to that of a T1 dwarf, it was also classified as a possible unresolved binary (Burgasser et al. 2010).
Although L ′ photometry was not included in this fit, we briefly compare the measured or estimated H − L ′ color of the four best-fit SpeX library sources in Figure 6 to that of HIP 39017b, H − L ′ = 1.63 ± 0.22.Of the four sources, only SDSSp J042348.57-041403.5 has measured MKO L ′ photometry (Leggett et al. 2002).For the remaining three, we use the photometric transformation relation from Leggett et al. (2019) to estimate their L′ from Spitzer 3.6 and 4.5µm photometry (Kirkpatrick et al. 2021; Spitzer Science Center (SSC) & Infrared Science Archive (IRSA) 2021; Leggett et al. 2007), combined with measured MKO H-band photometry (Zhang et al. 2021;Chiu et al. 2006).We find that the measured and derived H − L ′ colors of these four sources are somewhat redder than HIP 39017b, with Figure 6.The four smoothed and binned, best fit SpeX Prism spectra (in color) compared to the HIP 39017b Broadband CHARIS spectrum from December 2022 (black), plotted with vertical offset for visualization.Uncertainties on the SpeX Prism library spectra are plotted, but errorbars appear small compared to the marker size, due to CHARIS's larger spectral resolution.NIR spectral types are as given in the SpeX Prism library (Looper et al. 2007;Burgasser et al. 2010Burgasser et al. , 2008Burgasser et al. , 2004)).
values between 2.1 − 2.5 mag.However, a comparison of the absolute magnitudes is complicated by the unknown binarity of most of these library sources.In fact, SDSSp J042348.57-041403.5 has been resolved into a close L6.5+T2 binary (Dupuy & Liu 2017).

Comparison to Atmospheric Model Spectra
We compare the HIP 39017b spectrum from the December 2022 CHARIS observations along with the Keck/NIRC2 L ′ photometry to three model atmosphere grids: the Sonora Bobcat models (Marley et al. 2021), the Sonora Diamondback models (Morley et al., 2024,in preparation), and the Lacy/Burrows models generated for the analysis of HIP 99770b (Currie et al. 2023b).Fitting followed a similar procedure as the empirical spectral fits described in Section 3.1, with a few caveats.Due to the lack of uncertainty arising from the model spectral points, the covariance matrix, C −1 k , used when calculating the fit is just that of the measured CHARIS JHK spectrum.The calculation of χ 2 from Equation 1 also includes the additional photometric term, (f ν,L ′ −αF ν,L ′ ) 2 /σ 2 F ν,L ′ .Each model spectrum is fit to the observed spectrum by varying the scaling parameter for the entire model spectrum, α.Additional parameters for each model spectrum are as defined in the pre-generated model grid.A summary of the best fit model parameters and their χ 2 ν for each model grid is given in Table 3.
The Sonora models are a family of self-consistent atmospheric models for L, T, and Y dwarfs and selfluminous, young planets (Marley et al. 2021;Karalidi et al. 2021).Bobcat, the first of the Sonora models released, invokes equilibrium chemistry and cloudless atmospheres (Marley et al. 2021).The released spectral grid (Marley et al. 2021)  To investigate an alternative cloud parameterization, we also compared the HIP 39017b spectrum to the Lacy/Burrows model grid derived for comparison to HIP 99770b (Currie et al. 2023b).The Lacy/Burrows models expand on those of Burrows et al. (2006) with updated molecular line lists and absorption cross sections, as well as non-equilibrium carbon chemistry (Currie et al. 2023b).While less expansive in their coverage of the T ef f (∼ 1100 − 1600K) and log g (4.0, 4.5, 5.0) parameter space, they do include the (T ef f , log g) that provided the best fit in the cloudy Sonora Diamondback models.Further, cloud parameterization is differentiated into the modal size of the atmospheric dust particles, (in microns, Currie et al. 2023b) and cloud structure:AE, AEE, and E, in order of increasing sharpness of the cutoff at the top of the clouds (Madhusudhan et al. 2011).As these cloud models have the same interior profile, a more gradual cutoff with height means the presence of more cloud material and therefore a higher total optical depth (Madhusudhan et al. 2011).
The Lacy/Burrows model spectra from Currie et al. (2023b) contain grids of AE-and AEE-type clouds with a modal dust size of a 0 = 100µm with default methane abundances, as well as grids of AEE-type, a 0 = 100µm clouds and E-type, 60µm clouds with CH 4 abundance reduced by a factor of 10.
The Sonora Bobcat model that best fits the HIP 39017b spectrum has T ef f = 1700K, log g = 4.0, and C/O= 1.0, albeit with a significantly higher χ 2 ν (2.82) than is achieved by either of the other models.Notably, the best fit Bobcat models all fall short of the high peak flux observed around λ ∼ 2.07 − 2.14µm in the HIP 39017b spectrum (see Appendix B).
The results of the Sonora Diamondback fits are shown in Figures 7 and 8.The addition of clouds with the Sonora Diamondback models provide a smaller χ 2 ν .The  The results of our fit to the Lacy/Burrows models are shown in Figure 9.The best-fit Lacy/Burrows model from the HIP 99770b model grid for HIP 39017b is 1500 K, log g = 4.5 with solar CH 4 and AEE-type clouds with 100µm dust.Compared to the Sonora Diamondback model fits above, this represents a slightly warmer and higher gravity minimum than found for the extensive Diamondback grids (1300 K, log g = 4.0, f sed = 3).While no significant minimum is seen at 1300 K and log g = 4.0 in the Lacy/Burrows grid, this point is at the boundary of the parameter space covered for all but one of the cloud grids.
However, as can be seen from the fine coverage of the Diamondback grids, it is not unusual for the location of the minimum χ 2 ν to vary with cloud strength; in fact, local minima at higher T ef f ∼ 1500 K and log g ∼ 4.0 can be seen in several of the Diamondback grids in Figure 7 (e.g.f sed = 4 and [M/H]= 0).Therefore, the difference in parameters at which the best fit is achieved may be due in part to a difference in parameter space coverage.
The substantial improvement in the best-fit χ 2 ν between the cloudless Sonora Bobcat models and the cloudy Sonora Diamondback and Lacy/Burrows models implies that clouds are not only present on HIP 39017b, but are also having a significant impact on its NIR spectrum.While the best-fit models in all three grids analyzed here have similar log g, [M/H], and C/O, the T ef f of the best-fit model is lower for the cloudy models than the cloudless grid.Among the Sonora grids, the inclusion of clouds causes the best fit T ef f to shift from 1700 K in Bobcat to 1300 K in Diamondback, the latter of which is perhaps more consistent with its apparent position near the L/T transition in Section 3.1 (Kirkpatrick 2005).
In addition, the radius implied for the best-fit spectrum in each suite of model grids is shown in Table 3, as derived from the scaled model flux and the Gaia Data Release 3 (DR3) distance for the system (d = 65.8835± 0.0916 pc; Gaia Collaboration 2022).In calculating R, we include an estimated 7% systematic uncertainty in the amplitude of the measured flux.This combines a ∼ 5% uncertainty from the intrinsic variability of brown dwarf spectra (Cruz et al. 2018) and a ∼ 5% calibration uncertainty in the flux of the CHARIS ν fits for each available model spectrum as a function of T ef f and log g as above.Cloud parameterization is labeled for each subfigure (AE, AEE, or E, as well as modal dust size, in microns), with 'rCH4' indicating grids derived with reduced methane.(Bottom) Six of the Lacy/Burrows model spectra, smoothed to R = 50 at 2.65µm and overplotted on the observed HIP 39017b spectrum used for fitting.The legend gives the T ef f , log g, cloud prescription, whether they have reduced CH4, and its χ 2 ν fit to the observations.The first four spectra shown are the four Lacy/Burrows spectra with the lowest χ 2 ν , while the last two are those with the lowest χ 2 ν of the remaining two cloud prescription/methane abundance grids.
JHK spectral points due to variability in the satellite spot flux ratios (Wang et al. 2022).
Scaled model spectrum flux often yields unphysically small radii for brown dwarfs (e.g.Zalesky et al. 2019;Gonzales et al. 2020;Hood et al. 2023;Zhang et al. 2023), which is generally attributed to uncertainty or absent physical mechanisms in the models (Kitzmann et al. 2020;Lueber et al. 2022).Indeed, for the three model spectrum grids analyzed here, R ranges from 0.51 − 0.87 R J .For comparison, log g = 4.0 (4.5) for a 30M J object would correspond to a radius of ∼ 2.7R J (1.5R J ), while even a 13M J object would correspond to a radius of ∼ 1.8R J (1.0R J ).The two cloudy models do seem to underestimate the radius (or equivalently, overestimate the flux) less than the cloudless Bobcat model, as expected with the improved χ 2 ν of their spectral fits.In comparison, a background object seen at ≈1. ′′ 7 in some data sets follows the expected background object track (see Appendix A).
Following Nielsen et al. (2017) and Currie et al. (2023a), we compute the relative probability of HIP 39017 b being a background star (P (BG)P (ρ|BG)P (ν|ρ)) with a non-zero common proper motion versus a comoving companion P (comp) × P (ρ|comp).We adopt a Besancon model, modeling the background star distribution over a squaredegree area with distances ranging from 0 kpc to 1 kpc (π > 1 mas).Our model considers 3,085 stars as faint as H ∼ 20.6, or ∼2×10 −6 times the brightness of the primary (roughly our contrast limit at HIP 39017 b's separation).
As found in Currie et al. (2023a), the true probability of a background star is lower if HIP 39017 b's spectral type is considered.For a flat detection limit of 2×10 −6 from 0. ′′ 25 to 1. ′′ 1 and absolute H-band magnitudes from Pecaut & Mamajek (2013), a T0 dwarf must be closer than ≈ 250 pc.Adopting the space density of late L/early T dwarfs of ≈2×10 −3 pc −3 , we estimate a contamination rate of ≈10 −6 on the CHARIS field and ≈5×10 −5 for our entire survey (≈60 stars) to date.

Dynamical Modeling
Using the relative astrometry derived in section 2.3 along with proper motion measurements of the stellar host from Gaia DR3 (+30.22 mas yr −1 in R.A., and -42.53 mas yr −1 in decl.), we have verified a common proper motion between the host star HIP 39017 A and HIP 39017 b.Our observations also revealed a wider separation candidate companion near the edge of the field of view in all datasets.Through detailed analysis of our direct imaging observations and relative astrometry, this peripheral object has been conclusively classified as a background star demonstrating proper motions of the same magnitudes but directly contrary to the movement of the HIP 39017 Ab system (see Appendix A for a detailed explanation).
The absolute astrometry from HGCA (eDR3; Brandt (2021)) complements the newly acquired relative astrometry in facilitating the orbital retrieval of HIP 39017 A. HIP 39017 A is a γ Doradus variable known to exhibit large stellar pulsations with flux variations of 0.6-6% or radial velocity (RV) variations in the range of 2-4 km s −1 .The only four epochs of existing RV data (De Cat et al. 2006;Henry et al. 2011) show significant fluctuations within the range of 2-4 km s −1 , rendering these RV data not useful for orbit analysis.In the absence of RV data, we conduct orbital analysis solely based on relative astrometry of HIP 39017 b from high-contrast imaging and absolute astrometry of HIP 39017 A from the HGCA.

Methodology
We use the open-source Bayesian orbit modeling code orvara (Brandt et al. 2021b) to infer the orbital parameters of HIP 39017 b within an Markov Chain Monte Carlo (MCMC) framework.orvara utilizes the emcee (Foreman-Mackey et al. 2013) ensemble sampler for parallel tempering MCMC (PT-MCMC).This technique involves multiple parallel chains with varying temperatures, facilitating efficient sampling near the minimum χ 2 .Chains at different temperatures swap positions periodically to enhance information sharing and sampling efficiency.Typically, the coldest chain in PT-MCMC efficiently explores the parameter space and finds the global optimum likelihood.The parameter space orvara samples is nine-dimensional, including semi-major axis, eccentricity, reference epoch longitude (λ ref ), inclination, argument of periastron (ω), and longitude of ascending node (Ω)), system masses (M pri and M sec ), and an RV jitter term.The code enhances computational efficiency by analytically marginalizing auxiliary model parameters such as RV zero-points, parallax, and barycentric proper motion.
We analyze the orbital properties and mass of HIP 39017 b in several steps as follows.First, we perform a simple joint orvara fit to the full datasets of relative and absolute astrometry to derive posterior distributions for orbital parameters, the primary mass, and companion mass.Over the course of one year, HIP 39017 b moves little between subsequent CHARIS and NIRC2 data sets, in contrast to the behavior for other recently imaged companions like HIP 99770 b and HIP 21152 B (Currie et al. 2023b;Kuzuhara et al. 2022).Thus, we expect that constraints on at least some HIP 39017 b parameters may be poorer.Therefore, we then compare HIP 39017 b's luminosity and estimated age (≲ 115 Myr) to predictions for hot-start luminosity evolution models.These models give another assessment of the plausible mass range for HIP 39017 b and may identify acceptable dynamical solutions that are strongly disfavored.We down-sample the orvara results with informative posteriors, removing solutions whose companion masses are inconsistent with luminosity evolution analyses.

Distribution
We employ a broad informative Gaussian prior for the primary mass, a log-normal prior for mass and semimajor axis, and uninformative priors for other parameters (Table 4).The mass we adopt for HIP 39017 A is chosen to be 1.54 ± 0.25 M ⊙ to appropriately encompass reported literature values and its CMD position as dis-cussed in Section 2. Our PT-MCMC process uses 100 walkers and 30 temperatures across 10 5 steps.The convergence diagnostic of our MCMC simulations is based on the auto-correlation time-based criterion, where we evaluate the ratio of the total number of samples to the integrated auto-correlation time (N/τ int ), with a threshold value of 50 (ac=50; Hogg & Foreman-Mackey 2018).We discard 40% of the coldest chain as burn-in so that the remaining chain converges adequately for posterior orbital parameter inference.
The posterior distributions of selected parameters from our simple joint fit to absolute and relative astrometry are shown by a corner plot in Figure 11.HIP 39017 b's orbit has a semimajor axis of ∼ 23.8 +8.7 −6.1 au and is likely within ∼20 • of being edge on (i ∼ 84 +14 −22 degrees).The dynamical mass for HIP 39017 b retrieved from this raw fit is poorly constrained, with a wide 68% confidence interval of M comp = 30 +31 −12 M Jup , as is the eccentricity (e = 0.65 +0.27  −0.45 ).The implied mass ratio (q ∼ 1.9% +1.9 −0.8 ).The companion's mass, mass ratio, and eccentricity from this modeling are highly uncertain because HIP 39017 b was evidently observed during a period when the orbital motion is slow, resulting in the three epochs of relative astrometry being closely spaced and the orbital phase being sparsely sampled.Consequently, the other parameters we derive from the maximum likelihood orbit given in Figure 11 are also subject to biases.

Dynamical Modeling Results: Informed Posteriors
To further constrain the orbital mass of the companion, we first investigate the use of information from evolutionary models to downsample the original orbital posteriors derived above.Given age and bolometric luminosity, substellar evolutionary models predict the mass based on cooling physics.We use bi-linear interpolation to form a grid of mass predictions as a function of age and luminosity.The primary source of uncertainty in evolutionary model-predicted mass estimates is due to the uncertainty in age.We place a conservative upper age limit of ≤ 700 Myr (see Section 2.1) while we obtain the bolometric luminosity from both NIRC2's singleband L ′ photometry and the entire broadband spectrum, by comparing them to observed L and T dwarf spectra.
We consider two derivations of HIP 39017 b's luminosity.First, we convert the L ′ or W1 band magnitude, 16.12 ± 0.18 mag, to an absolute magnitude in W1 using the relation in Franson et al. (2023a) and subsequently to a bolometric luminosity following the same a Derived using uniform, informative mass posteriors over a range of Mp = [0,46] MJ anchored by hot-start evolutionary model predictions (see Figure 12).
approach as Li et al. ( 2023) using a comparative sample of L and T dwarfs from Filippazzo et al. (2015).
The bolometric luminosity obtained using this method is log(L bol /L ⊙ ) = −4.5 ± 0.1 dex for HIP 39017 b.Second, we consider an atmospheric modeling-derived estimate the bolometric luminosity from the best fit Sonora Diamondback (Morley et al., 2024,in preparation) and Lacy/Burrows models for HIP 99770b (Currie et al. 2023b) described above in Section 3.2.The flux of each best fit model is scaled by the fitted amplitude and integrated over the full extent of wavelengths available for each model (0.30 − 250 µm for Sonora Diamondback, 0.43−300 µm for Lacy/Burrows 99770).The luminosity is then calculated using the flux radius, R F ν , in Table 3.This results in a log(L bol /L ⊙ ) ∼ −4.7 dex for the best fit Sonora Diamondback model, and ∼ −4.6 dex for the best fit Lacy/Burrows 99770 model.These estimates are subject to model uncertainties, and could be slightly lower than the true bolometric luminosity due to unmeasured contribution from wavelengths outside of the spectral range of the models.However, they are consistent with the empirically-determined bolometric luminosity of HIP 39017 b determined above.
Considering both empirical and model-derived estimates, we derive a luminosity of log(L bol /L ⊙ ) = −4.6 ± 0.2 for HIP 39017 b.We compare the measured luminos-ity to several evolutionary tracks: ATMO2020 (Phillips et al. 2020), AMES-COND (Allard et al. 2001), AMES-DUSTY (Allard et al. 2001), BHAC15 (Baraffe et al. 2015), and BT-Settl (Allard et al. 2012(Allard et al. , 2013)).Using both the provisional age of 115 Myr and the upper limit age of 700 Myr, we determine the mass inferred from each of these models given HIP 39017 b's luminosity and its associated uncertainty (see Figure 12).The median of these comparisons yields a modelimplied mass of 13.8 +1.3  −1.6 M J at 115 Myr or 38.2 +7.9 −5.2 M J at 700 Myr.The corresponding model-inferred parameters are T ef f = 1200 +100 −120 K, log(g) = 4.41 +0.08 −0.05 dex, and R = 0.85 +0.53  −0.0047 M J if assuming an age of 115 Myr.Alternatively, an age of 700 Myr yields parameters of T ef f = 1330 +170 −150 K, log(g) = 5.02 +0.08 −0.07 dex, and R = 0.72 +0.0092 −0.014 M J .The evolutionary-model inferred masses are then used to truncate the posterior orbital parameter distributions for HIP 39017 b by down-sampling orbits within a range between zero and the conservative upper mass limit of 46 M J . Figure 13 shows the joint fit to the relative astrometry and absolute astrometry data points after resampling the original fit to the full datasets.The period of the companion is comparable to the 25 yr temporal baseline between Hipparcos and Gaia, which limits our ability to characterize the orbital shape (or eccentric- Figure 11.Posterior distributions of selected parameters from a joint fit with unweighted companion mass posteriors and using the full datasets of relative astrometry from direct imaging and complementary absolute astrometry from the HGCA.The contours show 1, 2, and 3σ contour levels at 68%, 95% and 99.7% confidence levels. ity) and dynamical mass.Posterior distributions of relevant model parameters after down-sampling are shown in Figure 14.A summary of the modeled and derived parameters before and after sampling are given in Table 4.
Compared to the previous fits, the mode of the mass distribution remains unchanged but the long tail of higher masses in conflict with luminosity evolution estimates is truncated.Thus, while most parameter distributions are consistent with results from the full posterior distribution fit, the dynamical mass is reduced from 30 +31 −12 M Jup to 23.6 +9.1 −7.4 M Jup after downsampling over the range of masses allowed by evolutionary track predictions.The derived mass ratio drops to 0.0148 +0.0065 −0.0049 .In contrast to the full posterior fit, the mass and mass ratio median and 68% confidence intervals lie below empirically-motivated thresholds -<25 M J and q <0.025 -to identify massive planets and distinguish them from low-mass brown dwarfs (Currie et al. 2023b,a).

DISCUSSION
We report the discovery of a substellar companion orbiting the young accelerating γ Doradus-type variable star, HIP 39017.Empirical comparisons to HIP 39017 b's photometry and spectrum indicate that it lies at the L/T transition with a best-estimated spectral type of T0.5-T2.Atmospheric modeling with the Sonara Diamondback grid achieves a good fit and indicates that HIP 39017 b is a cloudy object with a bestestimated temperature, gravity, and radius of T eff = 1300 K, log g = 4.0, and R= 0.87R J .Best fit values from the Lacy/Burrows cloudy atmosphere models are similar.
Modeling astrometric data for HIP 39017 b and the primary constrain the companion's semimajor axis to be ≈ 23 au and its orbit to be near edge-on, but otherwise yield poor constraints on system properties.The full posterior distribution for HIP 39017 b's dynamical mass and mass ratio have a wide range of 30 +31 −12 M J and 1.9% Figure 12.Mass evolution isochrones for 0.12-Gyr and 0.7-Gyr objects given by several different hot-start evolutionary models.The positions of HIP 39017 b are highlighted with yellow and cyan points along the substellar iso-age evolutionary tracks.With an estimated bolometric luminosity of -4.6±0.2 dex, the predicted masses at ages of 0.12 Gyr and 0.7 Gyr have median values of 13.8 +1.3 −1.6 MJ and 38.2 +7.9 −5.2 MJ , respectively.The mass uncertainties were estimated by adding the dispersion across the considered evolutionary tracks and the luminosity-related mass uncertainty in quadrature.We adopt 0.7 Gyr/38.2+7.9 −5.2 MJ as the upper age/mass limit.The final conservative upper mass limit from evolutionary models is then 46MJ .We adopt 0.12 Gyr/13.8+1.3  −1.6 MJ -the closest grid value to 115 Myr -as a more realistic characterization for HIP 39017 b.We note that a younger age than 115 Myr may indicate a even lower mass planetary companion as suggested by evolutionary model predictions.
tation.However, given HIP 39017 b's best-estimated age (≲ 115 Myr), luminosity evolution models favor a lower mass of ≲ 13.8 M J .Adopting the highlyconservative assumption that HIP 39017's age is simply "younger than the Hyades," we estimate an upper mass limit of ≈ 38 M J .Informed by the latter results -i.e.removing solutions from the posterior distribution with masses incompatible with luminosity evolution estimates -the revised posterior distributions for HIP 39017 b yield a lower mass of 23.6 +9.1 −7.4 M Jup and a lower mass ratio of 0.0148 +0.0065 −0.0049 .These combined values place this companion in a region of the parameter space populated by massive planets detected by RV and direct imaging, favoring a planetary interpretation as described in Currie et al. (2023b,a).
Current CHARIS and NIRC2 data for HIP 39017 b suggest that our images capture a part of its orbit where it undergoes little astrometric motion.Continued astrometric monitoring of HIP 39017 b over a long time baseline will substantially improve constraints on its orbit and, potentially, its dynamical mass.The Gaia Data Release 4 will yield more precise astrometric data for the star, which will further improve mass estimates and thus clarify its nature (i.e.whether it shares more of its properties with low-mass brown dwarfs or planets).
A more precise age for the system is needed for comparison to evolutionary models.Our estimate of ≲ 115 Myr implies a lower mass of ≲ 13.8M J from comparison with several hot-start evolutionary models (Figure 12), placing HIP 39017b closer to the deuterium burning limit.Its host star's multi-frequency, γ Dortype variability opens the door to deriving independent age constraints for the system through asteroseismology.Further observations with higher spectral resolution would allow the derivation of stronger constraints on its spectral properties (e.g.T ef f , log g, metallicity, C/O ratio) via atmospheric retrieval.
The discovery of HIP 39017b reinforces the improved efficiency of targeted direct imaging campaigns informed by long-baseline, precision stellar astrometry instead of so-called blind or unbiased surveys.Between the system's youth and potential for independent dynamical mass constraints, HIP 39017b provides a new testing ground for substellar evolutionary and atmospheric models.HIP 39017 b joins recently-discovered planetary companions HIP 99770 b and AF Lep b and brown dwarfs HIP 21152 B and HD 33632 Ab as companions with spectra that can constrain their atmospheres and direct dynamical mass measurements.With improved age constraints for all companions, new targeted surveys will establish a sample of companions that will help map how the atmospheres of planets and brown dwarfs of a given mass evolve with time.The 2022 NIRC2 data identify a second point source in the northwest at a wider separation (ρ ∼ 1. ′′ 75), which is recovered in the March 2022 CHARIS epoch and 2023 NIRC2 epoch (Figure 15).If a bound companion at a projected physical separation of 115 au, this object would contribute to the HIP 39017 primary's astrometric acceleration and thus substantially the orbital solutions for HIP 39017 b.However, our analysis below shows that this object is an unrelated background object whose colors are consistent with an unreddened K star.
Table 5 lists astrometric measurements for this object.While basic processing for each data set listed follows methods described in Section 2.2, we did not apply PSF subtraction techniques to detect the wider-separation object in the March 2022 CHARIS and March 2023 NIRC2 data.For these data sets, we extract an object spectrum and/or measure its astrometry directly.The detection signifances are sufficiently high (SNR ∼ 47-222) that systematics in the plate scale and north position angle calibration for CHARIS and NIRC2 dominate our astrometric error budget.
The wider-separation object (Figure 16) shows a spectrum that is smooth and featureless except for channels with substantial telluric contamination (1.3-1.4 µm, 1.8-1.9µm).We measure its broadband photometry to be J = 14.38 ± 0.08, H = 13.89 ± 0.04, and K s = 13.81 ± 0.04.At L ′ , we measure a brightness of 13.83 ± 0.04 magnitudes.The resulting J-K s color is 0.57 ± 0.09, roughly consistent with a background K2-K3 star (Pecaut & Mamajek 2013) at a distance of ≈850 pc.Aside from channels with significant telluric contamination, the 5000 K, log g = 4.5 Kurucz atmosphere model overplotted on Figure 16 matches this object's spectrum.Similarly, a proper motion analysis shows that this object's motion relative to HIP 39017 A is consistent with a background star  For comparison, the fits to the SpeX Prism L/T dwarf library and the Sonora Diamondback model grid discussed in Sections 3.1 and 3.2, respectively, were also performed using the February 2022 SCExAO/CHARIS spectrum.As seen in Figure 3, the shape of the ∼ 2.1µm peak in the February 2022 spectrum more closely resembles that of, for example, the empirical SpeX spectra in Figure 6 than the December 2022 spectrum did.However, the weaker J-band peak in the December spectrum is obscured in the February spectrum.In general, we find that the parameters derived using the February spectrum are more poorly constrained and the data are somewhat overfit.Therefore, we defer to  Note-S/N represents the object's signal-to-noise ratio, ρ represents the companion.
the parameters derived from the December spectrum.However, for completeness, we briefly discuss the analysis of the February spectrum here.
Consistent with the increased spectral noise of this data set (Section 3), the NIR spectral type implied by the SpeX Prism library spectra is more poorly constrained for the February 2022 data than for the December 2022 data.As can be seen in Figure 18, it is consistent with spectral types from ∼L5-T1 within 1σ, or from ∼L2-T2 within 2σ.Further, the higher uncertainty on this spectrum results in the observations being overfit by the single parameter fit to the library spectra, with many having χ 2 ν < 1.The February CHARIS spectrum, in combination with the Keck/NIRC2 L ′ photometry, is also slightly overfit by the Sonora Diamondback grid, but to a much lesser extent, with the minimum χ 2 ν = 0.94 (Figure 19).The Diamondback model that best fits the February 2022 spectrum and L ′ photometry is T ef f = 1400 K, log g = 4.5, [M/H]= 0.0, with f sed = 4 clouds.This is adjacent to the best fit Diamondback model using the December spectrum, but slightly higher temperature and log g.

Figure 5
Figure 5. χ 2 ν of each SpeX Prism template spectrum fit to the HIP 39017b CHARIS spectrum versus its NIR spectral type as given in the SpeX Prism library.Template spectra with 'peculiar' NIR spectral type are shown as diamonds.Dotted lines indicate the χ 2ν values that correspond to the 1, 2, and 3 σ confidence intervals, from bottom to top.

a
Whether or not reduced CH4 abundances were reduced by a factor of 10.Grids with reduced CH4 only available for the Lacy/Burrows 99770 models.b Only value of this parameter available in the model grid.Note-Companion radius, R, is derived from the fitted model flux using the Gaia DR3 system distance (d = 65.8835± 0.0916 pc; Gaia Collaboration 2022).

Figure 7 .
Figure 7. Maps of the χ 2 ν derived from fitting the HIP 39017b spectrum to the Sonora Diamondback (Morley et al., 2024,in preparation) model spectrum grids as a function of model log g and T ef f .Columns of subfigures show the three available metallicities: [M/H] = −0.5 (left column), 0.0 (middle column), and +0.5 (right column).Rows of subfigures show four of the six different cloud prescriptions: f sed = 1 (top row), 3 (second row), 4 (third row), and no clouds (fourth row).

Figure 8 .
Figure 8. Select Sonora Diamondback (Morley et al., 2024,in preparation) model spectra, smoothed to R = 50 at 2.65µm, overplotted with the HIP 39017b observations.Left: The best fit Sonora Diamondback spectra with each metallicity provided in the spectral grid.The T ef f , [M/H], and χ 2 ν for each model are given in the legend.All three best-fit spectra have log g = 4.0.Right: Model spectra with the same T ef f (1300 K), log(g) (4.0), and [M/H] (0.0) as the best fit spectrum, but showing variation with the cloud prescription.Each spectrum's f sed value (or 'nc' for no clouds) and its corresponding χ 2 ν are given in the legend.HIP 39017b spectrum is best fit by the 1300 K, log g = 4.0, [M/H]= 0.0 model with f sed = 3 (χ 2 ν = 1.099).The only other model spectrum that falls within 1σ of this minimum has the same T ef f and log g, but [M/H]= +0.5 and f sed = 4.A comparison of these two model spectra (as well as the best fit model spectrum with [M/H]= −0.5) can be seen in the first panel of Figure8.These two best fit models both show a much more accurate fit to the peak flux observed at ∼ 2.07 − 2.14µm than the cloudless models, as well as a much higher peak ∼ 3.8µm.The second panel of Figure8shows how different model spectra with the same T ef f (1300 K), log g (4.0), and [M/H] (0.0) as the best fit model compare to the observations when the cloud f sed varies.As can be seen in the figure, clouds with f sed ∼ 3 − 8 are necessary to achieve the heightened K-band peak and slightly damped J-band peak detected in the HIP 39017b CHARIS spectrum from 2022 December under the Diamondback model.Clouds with f sed ∼ 1 − 2, on the other hand, induce too much reddening of the spectrum compared to our observations of HIP 39017b, flattening the JHK peaks significantly.The results of our fit to the Lacy/Burrows models are shown in Figure9.The best-fit Lacy/Burrows model from the HIP 99770b model grid for HIP 39017b is 1500 K, log g = 4.5 with solar CH 4 and AEE-type clouds with 100µm dust.Compared to the Sonora Diamondback model fits above, this represents a slightly warmer and higher gravity minimum than found for the extensive Diamondback grids (1300 K, log g = 4.0, f sed = 3).While no significant minimum is seen at 1300 K and log g = 4.0 in the Lacy/Burrows grid, this point is at

Figure 9 .
Figure 9.The fits results of the Lacy/Burrows models (Currie et al. 2023b) to the December 2022 SCExAO/CHARIS spectrum and the NIRC2 L ′ photometric point of HIP 39017b.(Top) χ 2ν fits for each available model spectrum as a function of T ef f and log g as above.Cloud parameterization is labeled for each subfigure (AE, AEE, or E, as well as modal dust size, in microns), with 'rCH4' indicating grids derived with reduced methane.(Bottom) Six of the Lacy/Burrows model spectra, smoothed to R = 50 at 2.65µm and overplotted on the observed HIP 39017b spectrum used for fitting.The legend gives the T ef f , log g, cloud prescription, whether they have reduced CH4, and its χ 2 ν fit to the observations.The first four spectra shown are the four Lacy/Burrows spectra with the lowest χ 2 ν , while the last two are those with the lowest χ 2 ν of the remaining two cloud prescription/methane abundance grids.
Figure10(top panel) compares HIP 39017 b's astrometry to the motion expected for a background star.We rule out a stationary background star at a >12-σ level.In comparison, a background object seen at ≈1. ′′ 7 in some data sets follows the expected background object track (see Appendix A).FollowingNielsen et al. (2017) andCurrie et al. (2023a), we compute the relative probability of HIP 39017 b being a background star (P (BG)P (ρ|BG)P (ν|ρ)) with a non-zero common proper motion versus a comoving companion P (comp) × P (ρ|comp).We adopt a Besancon model, modeling the background star distribution over a squaredegree area with distances ranging from 0 kpc to 1 kpc (π > 1 mas).Our model considers 3,085 stars as faint as H ∼ 20.6, or ∼2×10 −6 times the brightness of the primary (roughly our contrast limit at HIP 39017 b's separation).About 55 stars match HIP 39017 b's kinematics to within the 5-σ level (Figure10, bottom panel) for P (ν|ρ) = 55/3085.For a 2-σ range over our entire set of observations of ρ ∼ 0. ′′ 385-0.′′ 417 and assuming we would detect HIP 39017 b between 0. ′′ 25 and 1. ′′ 1 from the star, we find P (BG)P (ρ|BG) ∼ 1.9×10 −5 for a combined background probability of 3.4×10 −7 .For a frequency of 0.8% for 13 M J -80 M J companions and mass and semimajor axis distributions of dN comp /dM ∝ M −0.47 and dN comp /da ∝ a −0.65 , we derive a planet/brown dwarf probability of ∼8.6×10 −5 .Thus, the relative probability of a background star versus a planet/brown dwarf is ∼ 0.004: a bound companion is more probable at the ≈ 3-σ level.As found inCurrie et al. (2023a), the true probability of a background star is lower if HIP 39017 b's spectral type is considered.For a flat detection limit of 2×10 −6 from 0. ′′ 25 to 1. ′′ 1 and absolute H-band magnitudes fromPecaut & Mamajek (2013), a T0 dwarf must be closer than ≈ 250 pc.Adopting the space density of late L/early T dwarfs of ≈2×10 −3 pc −3 , we estimate a contamination rate of ≈10 −6 on the CHARIS field and ≈5×10 −5 for our entire survey (≈60 stars) to date.

Figure 10 .
Figure 10.Top: Proper motion analysis for HIP 39017 b compared to the motion of a background star with zero proper motion.Bottom: Proper motion analysis for a synthesized Besancon population compared to the 1, 3, and 5-σ contours for background stars needed to match HIP 39017 b's relative astrometry.

Figure 13 .Figure 14 .
Figure13.Weighted MCMC orbital posteriors using an uniform, informative range of [0,46] MJ on the companion mass.This is a result from down-sampling the original fit to the relative astrometry from direct imaging presented in Table2(left panel) and stellar absolute astrometry in R.A. and decl.from the HGCA eDR3 (right two panels).The red data points in the left panel are relative astrometry data driven by direct imaging observations while the two data points in cyan and orange in the right panels denote proper motion measurements from the cross-calibrated Hipparcos and Gaia eDR3 catalogs.

Figure 15 .
Figure 15.Detections of a wider-separation point source around HIP 39017 with Keck/NIRC2 (left and right panels) and SCExAO/CHARIS (middle panel) B. SPECTRAL ANALYSIS: COMPARISON TO CLOUDLESS SONORA BOBCAT GRID As discussed in Section 3.2, the results of fitting the December 2022 SCExAO/CHARIS data set and Keck/NIRC2 L′ photometry to the Sonora Bobcat model grid are shown in Figure 17.
C. SPECTRAL ANALYSIS OF THE FEBRUARY 2022 SCEXAO/CHARIS SPECTRUM Figure 16.Left: CHARIS spectrum of the wider-separation point source shown to be a background star (green) with a Kurucz model atmosphere corresponding to a K2-K3 star overplotted (magenta).Right: Proper motion analysis showing that this point source's astrometry is consistent with the motion expected for a background star.The dashed line connects the object's position in the first detected epoch to its predicted position in the 2023 Keck/NIRC2 data if it is a background star, a prediction that matches the matches the observations within errors (upper-right green dot).

Figure 17 .
Figure 17.The Sonora Bobcat (Marley et al. 2021) fit results to the HIP 39017b spectrum using the December 2022 CHARIS dataset and the Keck/NIRC2 L ′ photometric point.(Left:) Map of the χ 2 ν as a function of model log g and T ef f for Sonora Bobcat models with solar metallicity and C/O.The best fit Sonora Bobcat model is marked with the black diamond, with 1, 2, and 3σ contours around the minimum.(Right:) The best fit model spectrum (smoothed to R = 50 at 2.65µm) for each C/O value provided in the Sonora Bobcat model grid overplotted with the HIP 39017b observations.The legend denotes the T ef f , log g, C/O, and χ 2 ν for each model spectrum shown.Note: Sonora Bobcat spectra with C/O= {0.5, 1.5} are only available with log g = 5.0.

Figure 18 .
Figure 18.Same as Figure 5 but using the χ 2 ν derived from fitting the empirical spectra to the February 2022 CHARIS spectrum rather than the December 2022 data.

Figure 19 .
Figure 19.Similar to Figure 7, but using the February 2022 CHARIS spectrum rather than the December 2022 data to fit the Sonora Diamondback model grid.
, although available data for HIP 39017 are scarce.
(Jones 2016r hand, the star's Gaia CMD pMacintosh et al. 2015) P − R P = 0.42 -is nearly identical to that of the planet-hosting stars HR 8799 (40 Myr;Marois  et al. 2008) and 51 Eri (23 Myr;Macintosh et al. 2015)and consistent with or slightly below the Pleiades sequence, providing evidence of youth (Figure1).HIP 39017's position is slightly blueward of one Ursa Majoris member (blue star) with an interferometricallyestimated age of 440 Myr(Jones 2016) and thus, by implication, younger.It is significantly bluer than the entire Hyades (750 Myr) locus.Given its CMD position with respect to the Pleiades sequence, an age ∼ 115 Myr or lower for HIP 39017 is favored.While HIP 39017 is slightly lower metallicity ([Fe/H] = −0.26±0.13;Bruntt

Table 2 .
HIP 39017 b Astrometry a The wavelength range for CHARIS is 1.16-2.37µm, while the L ′ -filter's central wavelength for NIRC2 is 3.78 µm.

Table 3 .
Best Fit Spectral Models

Table 4 .
Orbital parameters for HIP 39017 b derived with orvara.

Table 5 .
HIP 39017 Background Star Detections The wavelength range for CHARIS is 1.16-2.37µm, while the L ′ -filter's central wavelength for NIRC2 is 3.78 µm.