The Epoch of Giant Planet Migration Planet Search Program. II. A Young Hot Jupiter Candidate around the AB Dor Member HS Psc

We report the discovery of a hot Jupiter candidate orbiting HS Psc, a K7 (≈0.7 M ⊙) member of the ≈130 Myr AB Doradus moving group. Using radial velocities over 4 yr from the Habitable-zone Planet Finder spectrograph at the Hobby–Eberly Telescope, we find a periodic signal of Pb=3.986−0.003+0.044 days. A joint Keplerian and Gaussian process stellar activity model fit to the radial velocities yields a minimum mass of mpsini=1.5−0.4+0.6 M Jup. The stellar rotation period is well constrained by the Transiting Exoplanet Survey Satellite light curve (P rot = 1.086 ± 0.003 days) and is not an integer harmonic nor alias of the orbital period, supporting the planetary nature of the observed periodicity. HS Psc b joins a small population of young, close-in giant planet candidates with robust age and mass constraints and demonstrates that giant planets can either migrate to their close-in orbital separations by 130 Myr or form in situ. Given its membership in a young moving group, HS Psc represents an excellent target for follow-up observations to characterize this young hot Jupiter further, refine its orbital properties, and search for additional planets in the system.


INTRODUCTION
The processes by which the closest-in (P orb ≲ 10 d) giant planets arrive at their present-day locations remain challenging to observationally constrain.Several formation and migration mechanisms have been introduced to account for this population of hot Jupiters (HJs; Daw-son & Johnson 2018).These include in-situ formation (e.g., Boley et al. 2016;Batygin et al. 2016), disk migration (e.g., Goldreich & Tremaine 1980;Ward 1997;Kley & Nelson 2012), or more "violent" scenarios of three-body dynamical interactions such as planet-planet scattering or von Zeipel-Lidov-Kozai interactions coupled with high-eccentricity tidal migration (e.g., Wu & Murray 2003;Fabrycky & Tremaine 2007;Triaud et al. 2010;Naoz et al. 2011).However, as the timescales of these processes (≈1 Myr-1 Gyr) are shorter than the typical ages of observed systems, our knowledge of HJ evolutionary pathways is limited.
Young giant planets are valuable because they provide a means to distinguish among these mechanisms, which operate on different characteristic timescales.However, the detection and characterization of young HJs has Tran et al. been hindered by the presence of stellar activity, making these systems challenging to find and validate.As a result, many proposed young hot Jupiters have been disputed, rejected, or lack independent confirmation (e.g., Hernán-Obispo et al. 2010;Figueira et al. 2010;van Eyken et al. 2012;Hernán-Obispo et al. 2015;Tuomi et al. 2018;Carleo et al. 2018Carleo et al. , 2020;;Bouma et al. 2020a;Carmona et al. 2023).Notable examples of young HJs include V830 Tau b, TAP 26 b, and CI Tau b (Donati et al. 2016;Johns-Krull et al. 2016;Yu et al. 2017).
More recently, a growing number of large transiting planets have been found orbiting young stars, such as K2-33 b (David et al. 2016;Mann et al. 2016), HIP 67522 b (Rizzuto et al. 2020), and TOI-837 b (Bouma et al. 2020b), with radii consistent with Jupiter-sized planets (R p ≥ 0.5 R Jup ).However, these systems lack precise mass measurements.As such, it is difficult to determine whether these systems are truly young HJs or instead are inflated Neptunes that have yet to shrink, either from gravitational contraction as they cool or through atmospheric mass loss (e.g., Fortney et al. 2007;Lopez et al. 2012;Owen & Wu 2013;Gupta & Schlichting 2019, 2020), as has been suggested for systems discovered in young clusters (e.g., Mann et al. 2017;Rizzuto et al. 2018).
Given these challenges and the intrinsically low occurrence rate of HJs at field ages (1.2 ± 0.38%; Wright et al. 2012), each young system provides valuable information about the origin and dynamical evolution of this class of giant planets.Building a statistically robust sample of young HJs provides clues about the migration history and physical mechanism producing close-in giant planets.
The Epoch of Giant Planet Migration planet search program is an ongoing long-baseline, near-infrared (NIR) precision radial velocity (RV) survey targeting young, nearby Sun-like stars (Tran et al. 2021).Our targets comprise bona fide and high-probability candidate members of ten young moving groups (YMGs)-AB Dor, β Pic, Carina, Carina-Near, Octans, Tuc-Hor, Pleiades, 32 Ori, Argus, and Pisces-with ages between 20 Myr and 200 Myr.Known binaries and fast rotators (v sin i * > 35 km s −1 ) are removed to produce a sample consistent with giant planet RV search programs of older field stars (e.g., Johnson et al. 2010).By observing in the NIR, as opposed to the optical, our program leverages the wavelength dependence of spot-driven activity signals to minimize RV contributions from astrophysical jitter.This opens the possibility of detecting massive, close-in young planets with RV semi-amplitudes comparable to the reduced activity signals.The goal of this survey is to measure the occurrence rate of close-in gi-ant planets at young ages, compare this frequency with similar measurements at older ages, and constrain the dominant timescale and mechanism of giant planet migration.
Here, we present the discovery of a HJ candidate orbiting the young star HS Psc with the Habitable-zone Planet Finder spectrograph (HPF) at McDonald Observatory's Hobby-Eberly Telescope (HET).In Section 2, we summarize the physical parameters and previous observations of HS Psc.We describe the Transiting Exoplanet Survey Satellite (TESS) photometry and the HPF RV observations of HS Psc in Section 3. Characterization of the system, including the host star, and modeling of the HPF RVs is presented in Section 4. In Section 5, we discuss the implications of this discovery and follow-up observations that would be helpful to further validate and characterize HS Psc b.Finally, we summarize our results in Section 6.
As a young, bright, nearby star, HS Psc has been observed several times by direct-imaging exoplanet search programs.Brandt et al. (2014) targeted HS Psc with with HiCIAO differential imaging instrument (Suzuki et al. 2010) at the Subaru Telescope as part of the Strategic Exploration of Exoplanets and Disks with Subaru (SEEDS) program.They detected a faint point source at a projected distance of approximately 200 AU (≈5.′′ 3) but did not recover it in follow-up observations.Bowler et al. (2015) also acquired deep observations of HS Psc with Keck/NIRC2 in K s band, reaching contrasts of ∆K s = 9.0 mag at 0. ′′ 5 and ∆K s = 11.6 mag at 1 ′′ .Bowler et al. (2019) obtained shallow optical high-resolution imaging of HS Psc with Robo-AO at the Palomar 1.5 m telescope.No nearby substellar or stellar companions were present within the detection limits of both datasets.
Gaia DR3 reports a re-normalised unit weight error (RUWE) of 1.055 for HS Psc (Gaia Collaboration et al. 2021, 2022).RUWE values greater than 1.4 indicate that the single-star model is a poor fit to the astromet- ric solution (Lindegren 2018;Stassun & Torres 2021), so the RUWE value for HS Psc does not indicate the presence of a binary companion.Furthermore, the RMS of our precise radial velocities over 4 years in Section 3.2 implies that we can rule out a close-in massive companion for most orbital orientations.Altogether HS Psc appears to be a single young star with no nearby binary companions in the stellar or brown dwarf (BD) regimes, within the detection limits of current imaging and RV surveys.

TESS Photometry
The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) observed HS Psc at 2-minute cadence in Sectors 17 (UT 2019 October 8 to 2019 November 2) and 57 (UT 2022 September 30 to 2022 October 29).The light curve was retrieved and processed following the procedure described in Bowler et al. (2023).We downloaded the Science Processing Operations Center (SPOC; Jenkins et al. 2016) Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) light curve (Smith et al. 2012;Stumpe et al. 2012Stumpe et al. , 2014) ) from the MAST data archive1 using the lightkurve (Lightkurve Collaboration et al. 2018) software package.All photometric measurements flagged as poor (DQUALITY > 0) or are listed as NaN are removed.Each TESS sector is first median-normalized and then stitched together.Outlier points from flares and artifacts are removed by flattening the light curve using a high-pass Savitzky-Golay filter (Savitzky & Golay 1964) and excluding all data outside of three standard deviations of the flattened light curve.The TESS light curve for Sectors 17 and 57 are displayed in the left panels of Figure 1.

HPF Spectroscopy
As part of our initial results from the Epoch of Giant Planet Migration program, Tran et al. (2021) measured an RV rms of over 200 m s −1 for HS Psc with 8 RV epochs.This system was an outlier among our sample of young stars, which had a median RV rms of only 34.3 m s −1 .As this anomalously high RV scatter could potentially be explained by a close-in giant planet, we increased our monitoring cadence of this system over the following 3 years.
Observations are obtained and RVs are measured following the procedures detailed in Tran et al. (2021).HET's flexible scheduling system allows observations to be taken in queue mode (Shetrone et al. 2007).During each epoch, 3 contiguous exposures are obtained to sample high frequency variations.All exposures have a 300 s integration time and the average S/N at 1.07 µm is 115±30 per pixel.1D spectra are optimally extracted using the custom HPF data-reduction pipeline following the procedures described in Ninan et al. (2018), Kaplan et al. (2019), andMetcalf et al. (2019).
From UT 2022 May to June, HPF underwent a vacuum warm-up and cool-down cycle so maintenance work could be carried out on the detector.This has introduced a velocity offset in the RV time series.To account for this unknown offset, RVs obtained prior to and after this downtime are treated as if they are measurements from two different instruments when searching for planets and modeling Keplerian signals.In this work, we label the RV measurements prior to this maintenance "pre-warmup" HPF 1 and measurements after "post-warmup" HPF 2 to distinguish between the two observing periods.Altogether there are N RV, HPF1 = 62 and N RV, HPF2 = 21 observations associated with HPF 1 and HPF 2 , respectively.The pre-warmup RVs have an rms of RV rms, HPF1 = 262 m s −1 and an average RV measurement uncertainty of σRV, HPF1 = 124 m s −1 .The post-warmup RVs have an rms of RV rms, HPF2 = 264 m s −1 and an average RV measurement uncertainty of σRV, HPF2 = 150 m s −1 .These large measurement uncertainties are primarily driven by the large rotational velocity (see Section 4.1) of HS Psc.

Spectroscopic Parameters
We determine the spectroscopic properties of HS Psc following the empirical spectral matching procedure as described in Stefansson et al. (2020b), which is based on the SpecMatch-Emp algorithm of Yee et al. (2017).This technique compares a science spectrum with a li-brary of high-S/N spectra of slowly rotating stars and identifies the best match between the stellar and library spectra.During each fit, the reference star spectrum is convolved with a broadening kernel (Gray 1992) in order to estimate the rotational velocity.The final values are estimated from a composite spectra of the five best-fitting reference spectra.
Using the highest-S/N HPF (S/N = 150) spectrum as the science target, we measure the effective temperature (T eff ), metallicity ([Fe/H]), surface gravity (log g), and projected rotational velocity (v sin i * ) of HS Psc.The procedure is applied to the same 8 HPF échelle orders used for RV extraction.Uncertainties on T eff , [Fe/H], and log g are calculated for each échelle order using a cross-validation method following Stefansson et al. (2020b), which iteratively removes each library spectrum from the sample, fits for T eff , [Fe/H], and log g using the interpolated grid, and compares the best-fit measurement to the known library value.The standard deviation of the residuals between the recovered and library values is adopted as the uncertainty of the parameters for each échelle order.The weighted mean and weighted standard deviation of the spectral order measurements are used as the final values and uncertainties for each parameter.This cross-validation method cannot determine uncertainties on v sin i * as all library spectra are of slowly rotating stars.Thus, for v sin i * , we adopt the median value and standard deviation of the 8 échelle order values.We find an effective temperature of T eff = 4203 ± 116 K, a metallicity of [Fe/H] = −0.05± 0.09 dex, a surface gravity of log g = 4.66 ± 0.03 dex, and a stellar rotational velocity of v sin i * = 29.7 ± 3.1 km s −1 for HS Psc.
These values are consistent with measurements reported in the literature.McCarthy & White (2012) found T eff = 4400 ± 105 K using low-resolution (R ∼ 3575) optical spectra.Stassun et al. (2019) report an effective temperature of T eff = 4215 ± 128 K and a surface gravity of log g = 4.64 ± 0.11 dex.Barenfeld et al. (2013) found a near-Solar representative metallicity for the AB Dor moving group of [Fe/H] = 0.02 ± 0.02 dex.Schlieder et al. (2010) report a much lower rotational velocity of v sin i * = 10 ± 2 km s −1 using spectra from CSHELL (Greene et al. 1993).This discrepancy between v sin i * measured using our HPF and their CSHELL spectra may be due to the lower resolution of CSHELL, which can be as low as R = 5,000, depending on the slit choice.In Tran et al. (2021), we previously measured a relatively uncertain projected rotational velocity of v sin i * = 22±8 km s −1 .Our new measurement of v sin i * = 29.7 ± 3.1 km s −1 is more precise and based on empirical templates instead of synthetic models.

Stellar Mass and Radius
We estimate the stellar radius by fitting the spectral energy distribution (SED) of HS Psc using the open software package ARIADNE (Vines & Jenkins 2022).ARIADNE employs a Bayesian model averaging approach to convolve four stellar atmosphere model grids-PHOENIX v2 (Husser et al. 2013), BT-Settl (Allard et al. 2011), Kurucz (1993), andCastelli &Kurucz (2003)-with response functions of common broadband photometric filters.The model leaves distance, stellar radius, line-ofsight extinction (A V ), and excess photometric uncertainty terms as free parameters.Synthetic SEDs are generated by interpolating the model grids in T eff -log g-[Fe/H] space and fit to available photometry.For our fitting process, we use the B, V , G BP , G, G RP , J, H, K s , W 1 , and W 2 bandpasses.
For T eff , log g, and [Fe/H], we adopt Gaussian priors with means and standard deviations set to the values and uncertainties reported in Section 4.1.The distance prior uses a Gaussian distribution with the mean and standard deviation set to the Gaia DR3 distance estimate and the larger of the upper and lower uncertainties, N (37.667,0.044)2 pc, the stellar radius prior set to a Gaussian distribution with the mean and standard deviation set to the Gaia DR3 Final Luminosity Age Mass Estimator (FLAME; Gaia Collaboration et al. 2018) derived radius and uncertainty, P(R * ) = N (0.6561, 0.019) R ⊙ , and the line-of-sight extinction prior set to a uniform distribution with an upper limit set to the maximum line-of-sight (A V, max = 0.235 mag) extinction from the Schlegel et al. (1998, SFD) galactic dust map (Schlafly & Finkbeiner 2011) as determined by the software package dustmaps (Green 2018).The excess photometric noise terms all have Gaussian priors centered at 0 with a standard deviation equal to 10 times the photometric uncertainty of each band.
The posterior distributions of all fitted parameters in the model are sampled using the dynamic nested sampler dynesty (Skilling 2004(Skilling , 2006;;Speagle 2020;Koposov et al. 2023).The sampling is initialized with 5000 live points and terminates when the evidence tolerance reaches a threshold of dlogz < 0.1.From this fit to the SED, we infer a stellar radius of R * = 0.65 ± 0.01 R ⊙ .The parameters derived from the grid-interpolation are T eff = 4232 +45 −38 K, log g = 4.66 ± 0.03 dex, and [Fe/H] = −0.08 +0.08 −0.07 dex, which are consistent with the values measured in Section 4.1.The bolometric luminosity computed from the SED fitting procedure is L * = 0.121±0.006L ⊙ .Using the derived radius and surface gravity, we estimate a mass of M * = 0.70±0.05M ⊙ for HS Psc.However, this mass is sensitive to the inferred surface gravity.
We infer the mass of HS Psc using the open software package isochrones (Morton 2015).isochrones determines fundamental stellar properties by fitting combinations of photometric bandpasses and physical values to synthetic values generated using interpolated grids of evolutionary models.For the fitting routine, we utilize the MESA Isochrones and Stellar Tracks (MIST; Dotter 2016; Choi et al. 2016) evolutionary model grids, all broadband photometry listed in Table 1, all parameters measured in Section 4.1, and the Gaia DR3 parallax.
As isochrones samples the age parameter as log 10 (age), we set a flat prior in log-space based on the AB Dor moving group age, P(τ * ) = log U[113, 148]3 Myr.We also adopt a Gaussian prior based on the representative metallicity of the AB Dor moving group from Barenfeld et al. (2013), P([Fe/H]) = N (0.02, 0.02) dex.For the distance and extinction priors, we adopt the same priors as in the SED fitting, P(distance) = N (37.667,0.044) pc and P(A V ) = U[0, 0.235] mag.Following the default distribution adopted by isochrones, the mass prior follows the lognormal initial mass function from Chabrier (2003, Equation 17).
The posterior distributions of all fit and derived parameters are sampled with the multimodal nested sampling algorithm MultiNest (Feroz et al. 2009(Feroz et al. , 2019) ) using the open software package PyMultiNest (Buchner et al. 2014).We initialize the sampling with 5000 live points.The best-fit stellar mass and radius are M * = 0.686 ± 0.003 M ⊙ and R * = 0.628 ± 0.002 R ⊙ .The best-fit inferred parameters are T eff = 4338 ± 9 K, log g = 4.679 ± 0.001 dex, and [Fe/H] = 0.10 ± 0.02 dex, which are consistent with the measurements adopted in Section 4.1 to within 2σ.Altering the age, metallicity, and distance priors do not substantially change these results.
As a comparison to these measurements, we also report masses and radii estimated using relationships with other parameters such as the effective temperature, metallicity, and luminosity.Combining the parameters derived in Section 4.1 with the empirical functions calibrated using 190 stars from Torres et al. (2010) relating T eff , log g, and [Fe/H] to M * and R * yields a stellar mass and radius of M * = 0.60 ± 0.04 M ⊙ and R * = 0.59 ± 0.02 R ⊙ .Similarly, using the Stefan-Boltzmann law, the SED-computed luminosity (L * = 0.121 ± 0.006 L ⊙ ), and our adopted effective temperature (T eff = 4203±116 K) returns a radius of R * = 0.66 ± 0.04 R ⊙ .Stassun et al. (2019) estimates a mass and radius of M * = 0.66±0.08M ⊙ and R * = 0.647±0.062R ⊙ , respectively, for HS Psc.
We adopt the stellar radius and mass inferred with ARIADNE and isochrones, respectively.Given the range of inferred masses and radii in our analysis and compared with similar dispersion in estimates for other K7V members of the AB Dor moving group, we conservatively adopt a flat 10% estimate for the uncertainty in both parameters.This is larger than the characteristic scatter of inferred masses and radii for comparatively old main-sequence stars (∼5%; Tayar et al. 2022).The final adopted stellar mass and radius for HS Psc is M * = 0.69 ± 0.07 M ⊙ and R * = 0.65 ± 0.07 R ⊙ .

Rotation Period, Stellar Inclination, and Transit Search
The TESS light curve of HS Psc exhibits clear and consistent modulation at the levels of 5% and 2.5% during Sectors 17 and 57, respectively.While the amplitude (and to a lesser degree the phase) of the signal changes between the two observational windows, its periodicity does not.This suggests that the observed variability is driven by long-lived spots on the stellar surface and therefore can be used to accurately estimate the stellar rotation period.
The rotation period is derived by independently fitting a quasi-periodic (QP) Gaussian process (GP) to each TESS Sector light curve.QP GPs have been shown to accurately infer stellar rotation periods from photometric time series (e.g., Angus et al. 2018).The TESS Sectors are treated separately as significant evolution has occurred in the 2.91 years between Sectors 17 and 57 (Figure 1).We bin the TESS photometry to 30-minute cadence to improve computational efficiency while retaining large-scale, activity-driven structure.We use a QP kernel of the form: where τ is the time interval between any two points in time t i and t j , |t i − t j |, A is the amplitude, l is the local correlation timescale, P rot is stellar rotation period, and θ is the smoothness of the periodic correlation.
We use the emcee open Python package (Foreman-Mackey et al. 2013, 2019) to sample the posteriors of the kernel hyperparameters and an additional whitenoise jitter term, σ jit .We impose non-informative uniform priors for each hyperparameter; for A, l, θ, σ jit , and P rot , we adopt P(A) = U[0.01,5000] ppt, P(l) = U[0.001,10], P(θ) = U[0.001,10], P(σ jit ) = U[log(0.001),log(100)] ppt, and P(P rot ) = U[0.1,10] d, respectively.Sampling is initialized with 50 walkers and 10,000 steps, for a total of 5 × 10 6 samples, the first 50% are removed as burn-in, and convergence is confirmed using the autocorrelation time.We adopt the P rot MAP value and standard deviation of the posteriors as the rotation period measurement and uncertainty.We find separate rotation periods of P rot, 17 = 1.086 ± 0.003 d and P rot, 57 = 1.086 ± 0.004 d for TESS Sectors 17 and 57, respectively.
These rotation periods are consistent with the values derived from the SuperWASP (P rot = 1.0852 d; Norton et al. 2007) and KELT (P rot = 1.0859 d; Oelkers et al. 2018) photometric surveys.Thus, we safely adopt the weighted mean and mean error of the two TESS Sectors, P rot = 1.086 ± 0.003 d, as the stellar rotation period of HS Psc.The processed TESS light curves and marginalized P rot posteriors for Sectors 17 and 57 are displayed in Figure 1.
Using the rotation period, stellar radius, and projected rotational velocity, we infer the line-of-sight stellar inclination, i * , following the Bayesian framework from Masuda & Winn (2020).We utilize the analytical expressions described in Bowler et al. (2023, Equation 9) to obtain the posterior distribution of i * while accounting for the correlation between the projected and equatorial rotational velocities.Here, we adopt the TESS rotation period of P rot = 1.086 ± 0.003 d, a stellar radius of R * = 0.65 ± 0.07 R ⊙ , and a projected rotational velocity of v sin i * = 29.7 ± 3.1 km s −1 .It is evident, given the projected rotational velocity, v sin i * = 29.7 ± 3.1 km s −1 , and the equatorial velocity, 2πR * /P rot = v eq = 29.7 ± 2.8 km s −1 , that HS Psc is seen approximately equator-on.We find a maximum a posteriori (MAP) value of 83 • , a 68% credible interval spanning 67-90 • , and a 95% credible interval spanning 52-90 • for the stellar inclination.
As a close-in giant planet, HS Psc b's geometric transit probability is relatively high at ≈6%.We use the Notch filter and Locally Optimized Combination of Rotations (LOCoR) detrending algorithms (Rizzuto et al. 2017) to search for potential transit events in the TESS light curve.The Notch and LOCoR filters use a combination of a transit-shaped notch and quadratic continuum to search for transit-like events in small windows of the light curve.These algorithms have been used to discover planets smaller than Jupiter transiting young (τ * < 100 Myr) stars (e.g., Rizzuto et al. 2018;Newton et al. 2019;Rizzuto et al. 2020).
We run the Notch filter to detrend the original TESS light curve4 and, as HS Psc is a fast rotator, we apply the LOCoR algorithm to account for aliases of its rotation period (P rot = 1.086 ± 0.003 d).Using various window sizes (0.05-0.5 d), we do not detect a signal at the periodicity of the planet (P b = 3.986 d) in the detrended light curves.Several strong signals do emerge at other periodicities, but inspection of the phase-folded, detrended light curves at those observed periods and times of inferior conjunction indicate that these signals are artifacts of the detrending process or remnants of stellar activity that were not detrended.

Periodogram and Stellar Activity Analysis
To search for periodic signals from planets and investigate the stellar activity of HS Psc, we compute the GLS periodogram of the RV observations and associated activity indicators over the frequency range 0.0005-1.0526d −1 (0.95-2000 d).This upper frequency limit near 1 day corresponds to the first peak of the spectral window function, which is caused by a strong sampling alias resulting from nightly observations frequently seen in ground-based observations.This frequency represents the lowest limit at which the data can reliably identify periodicities.The periodograms are calculated on the concatenated sets of measurements (HPF 1 and HPF 2 ).
To test if the choice to combine HPF 1 and HPF 2 measurements impacts our results, we calculate the GLS periodogram with both the LombScargle and the LombScargleMultiband classes in the open-source Python package astropy.timeseries(Astropy Collaboration et al. 2013Collaboration et al. , 2018Collaboration et al. , 2022)), which are designed to treat single time series and multiband time series, respectively.We also calculate GLS periodograms on the concatenated RVs after subtracting the average of each dataset and after subtracting the derived RV offsets.Furthermore, we also calculate periodograms on the RV data separated by even and odd observations and in two halves.We find that all periodograms are similar to each other and periodogram features remain robust against these changes.We thus safely choose to calculate the single-band periodogram on the concatenated measurements for further analysis.tion 4.5), and the spectral window function over the period range of 0.95-50 d.The FAPs at the 1% and 0.01% levels are calculated using the bootstrap method and are shown in dotted and dashed orange lines, respectively.Five peaks rise above the 0.01% threshold in the RV periodogram.A peak at 1.003 d corresponds to an alias caused by nightly sampling, which is also seen in the window function periodogram.Significant peaks appear near the rotation period (1.09 d) and a likely alias sequence of the rotation period (∼1.2 d).Finally, two significant peaks are observed at 3.99 d and a possible harmonic alias of this signal, 1.91 d.We interpret the ≈4-day signal as a planet candidate, as no significant peaks emerge in the activity indicators at or near this frequency.No additional periodic signals that may have been overwhelmed by the initial planetary signal emerge in the periodogram of the residuals; only the ∼1-day synoptic nightly sampling alias remains significant at the FAP = 0.01% level.
To further test for correlations between the RVs and associated spectral indicators, we compute the Pearson's correlation coefficient (Pearson's r) and pvalues between the RVs and activity indicators (see Section 3.2).This test is conducted as spectral indicators that trace distortions in the spectral line profile induced by starspots may correlate with RV measurements (e.g., Queloz et al. 2001;Boisse et al. 2011;Meunier & Lagrange 2013).Here, r values measure correlation strength, ranging from −1 for perfect anti-correlation and 1 for perfect correlation.p-values report the probability of observing a corresponding r value if the data are not correlated, so a lower p-value represents a higher probability that the data are correlated.For the HPF RVs and associated activity indicators, all Pearson's r values have an absolute value less than 0.5 and all pvalues are greater than 0.05, indicating that none of the activity indicators are significantly correlated with the HPF RVs.The maximum correlation coefficient is r = 0.4, with a p-value of 0.07, for the HPF 2 RVs and the corresponding Ca II IRT 2 line index.The other Ca II IRT indices have similar r values.

Radial Velocities and Gaussian Process Modeling
The youth and persistent photometric variability of HS Psc signifies that strong, correlated starspot-driven modulations should be present in the RV observations of HS Psc b.Red-noise GP models are now commonly employed to fit RV variability arising from stellar activity (e.g., Affer et al. 2016;Damasso & Del Sordo 2017;Faria et al. 2020;Stefànsson et al. 2022).To take into account the expected correlation of RV signals caused by stellar activity, we perform three separate model fits to the HPF RVs that vary in their treatment of stellar activity: Model 1. Keplerian-only, Model 2. Keplerian and Quasi-periodic GP, Model 3. Keplerian and Matérn-5/2 GP.
All three models fit for a Keplerian orbit with 5 parameters: orbital period (P b ), time of inferior conjunction (T 0,b ), RV semi-amplitude (K b ), and the parameterized forms of eccentricity and argument of periastron ( √ e b sin ω * and √ e b cos ω * ).Furthermore, for each set of HPF RV observations, we add a "jitter" term (σ HPF1 and σ HPF2 ) to account for any additional instrumental noise not represented in the measurement uncertainties, as well as zero-point velocity terms (γ HPF1 and γ HPF2 ) to account for any systematic offsets.Models 2 and 3 additionally include GP components with different kernels, as described in Sections 4.5.2 and 4.5.3.For all models, the parameter posterior distributions are sampled with a Markov chain Monte Carlo Metropolis-Hasting algorithm as described by Sharma (2017).Each distribution is compiled from 50 independent chains of 200,000 iterations and a thinning factor of 10.Convergence is determined using the Gelman-Rubin statistic, with all chains having values under 1.02 (Gelman & Rubin 1992).

Model 1: Keplerian-only
For Model 1, we perform a Keplerian orbit fit to the RVs for the 9 parameters described above.No stellar activity mitigation scheme is adopted in this fiducial fit.Moreover, to search for longer-term accelerations, we also separately carry out two additional fits that model linear and and linear-plus-quadradic terms ( γ and γ) for Model 1 only.
We impose uniform priors on all Keplerian model parameters.Based on the strong 4-day signal from the periodogram, we adopt P(P b ) = U[0.5, 10.0] d and P(T 0,b ) = U[2458483.0,2458487.5]d as priors for the orbital period and time of transit, respectively.We set the prior on K b to P(K b ) = U[1.0,1000.0]m s −1 .For the eccentricity parameters, we allow the full range of which results in a uniform sampling of P(e b ) = U[0, 1) and P(ω * ) = U[0, 2π], for e b and ω * respectively.For the jitter terms, we adopt modified Jeffreys priors, J (1, 100) m s −1 , as defined by Equation 16of Gregory (2005), The modified Jeffreys prior behaves like a uniform prior for σ < σ a and a Jeffreys prior for σ > σ a .The Jeffreys prior is scale invariant in that it treats each decade bin as having equal probability.For the offset terms, we use uniform priors bounded by the minimum and maximum of the RV measurements, P(γ HPFi ) = U[min(HPF i ), max(HPF i )].For the model fits that incorporate long-term trends, we adopt uniform priors of We report the results of the Model 1 fit in Appendix B. Table B2 summarizes the parameter constraints and Figure B1 displays the best-fit Keplerian signal, 1σ and 3σ confidence intervals, and associated joint posterior distributions in the form of a corner plot.We note that the model fits with linear and quadratic terms are consistent with zero acceleration and return nearly identical solutions.The resulting Bayesian Information Criterion (BIC; Schwarz 1978;Raftery 1986;Tran et al. 2022) for the Keplerian orbit-only, Keplerian orbit and acceleration, and Keplerian orbit, acceleration, and curvature models are −12.0,−8.2, and −3.8, respectively, which indicates that the data prefer the Keplerian orbit-only model.Thus, we report only the single-planet model without long-term acceleration or curvature terms.
The median and 1σ posterior values for the Keplerian semi-amplitude and orbital period are K b,Model 1 = 301 ± 27 m s −1 and P b,Model 1 = 3.986 +0.001 −0.001 d, respectively.Using our adopted stellar mass, M * = 0.69±0.07M ⊙ (Section 4.2), we find a minimum mass of m b sin i Model 1 = 1.78 +0.22  −0.21 M Jup for HS Psc b.The bestfit eccentricity points to a modest value of e b,Model 1 = 0.18, although the posteriors are broad and consistent with circular orbits.The rms of the RV residuals after subtracting the best-fit Keplerian curve is 178 m s −1 .
The inferred posteriors of the RV jitter terms are σ HPF1, Model 1 = 161.3+21.0 −24.6 m s −1 and σ HPF2, Model 1 = 17.4 +27.2 −17.4 m s −1 .The large difference between σ HPF1, Model 1 and σ HPF2, Model 1 is likely driven by stronger stellar activity during the HPF 1 observations.The evolution of stellar activity over time is evident from the TESS light curves.The amplitude of photometric variability is approximately a factor of 2 greater in Sector 17 compared to Sector 57.As the HPF 1 and HPF 2 RVs are approximately coincident with TESS Sectors 17 and 57, respectively, this evolution could be driving differences between the jitter terms.Furthermore, the lower number and shorter time baseline of the HPF 2 RVs results in a diminished exploration of the stellar activity.This effect is amplified as a nightly cadence over this short observational window results in a limited sampling of the stellar rotation period with the HPF 2 observations.When phased to P rot = 1.086 d, the HPF 2 RVs cover only ∼1/3 of the stellar rotation phase space.

Tran et al.
This poor phase coverage can artificially dampen activity signals in the HPF 2 measurements, leading to a lower jitter estimate.Altogether these factors suggest that stellar activity is prominent in the RVs and necessitate a mitigation strategy.

Model 2: Keplerian and Quasi-Periodic GP
In Model 2, we simultaneously fit for a Keplerian signal and a GP model defined by a QP kernel.Under similar assumptions as for photometry (see Section 4.3), the QP GP has been widely leveraged to model stellar activity in RVs (e.g., Haywood et al. 2014;Grunblatt et al. 2015;Faria et al. 2016;Cloutier et al. 2017;Dai et al. 2017).We adopt a GP model with a QP kernel as implemented in pyaneti: (3) where t and t ′ represent the pairs of observations in time, and A, λ e , P GP , and λ p are the same or analogous to A, l, P rot , and θ from Equation 1.
We adopt uniform priors on all Keplerian parameters.Based on the results from Model 1, we impose P(P b ) = U[3.5, 4.5] d and P(T 0,b = U[2458485.0,2458488.0]d as priors for the orbital period and time of transit, respectively.All other priors for planetary parameters are the same as in Model 1. For the QP kernel hyperparameters, we adopt P(A) = U[0.0,1000.0]m s −1 , P(λ e ) = U[0.01,5.0], P(λ p ) = U[0.01,5.0], and P(P GP ) = N (1.09, 0.02) d for A, λ e , λ p , and P GP , respectively.This last prior uses the MAP value and five times the standard deviation of P rot measured from the TESS photometry.This is based on the assumption that the stellar activity signal is starspotdriven and will modulate at the stellar rotation period in both photometric and RV time series (e.g., Aigrain et al. 2012;Tran et al. 2023).The larger range allow the model to explore potential deviations in the rotation period across the large observational window, for instance from mid-latitude spots that could produce a spread of periodic signals as a result of differential rotation.
Table B3 and Figure  For Model 3, we perform a similar joint Keplerianplus-GP model fit as Model 2, but with a GP defined by a Matérn-5/2 (M 5 /2) kernel instead of a QP kernel.Our goal is to assess how changes to the GP stellar activity model can impact the inferred planetary parameters.Furthermore, this change allows us to account for the possibility that the stellar activity signals manifest differently in the RVs as compared to the TESS photometry.The Matérn family of kernels are also widely adopted in GP regression to model more stochastic behavior.GP models incorporating the M 5 /2 kernel and its derivative have been found to reasonably match stellar activity in solar RV data (Gilbertson et al. 2020).
We apply a GP model with a M 5 /2 kernel as defined in pyaneti: , where λ is the length of local variations.
We adopt uniform priors for all parameters in the joint Model 3 fit.For the planetary parameters, we impose the same priors as in Model 2. For the M 5 /2 kernel hyperparameters, we choose P(A) = U[0.0,1000.0]m s −1 for A and P(λ) = U[0.001,10.0] for λ.
Table 2 details the results of the Model 3 fit.a U[a, b] refers to the uniform distribution bounded by a and b.J (a, b) refers to the modified Jeffreys prior as defined in Equation 6of Gregory (2005), b MAP refers to the maximum a posteriori value.
c Planetary mass is derived assuming a stellar mass of M * = 0.69 ± 0.07 M⊙.All inferred posteriors for planetary and instrumental parameters are consistent among the three models within 1σ except for the HPF 1 RV jitter term.The posterior values of σ HPF1 are similar for both Models 2 and 3, which are reduced by a factor of approximately 20 compared to the inferred Model 1 values.This suggests that the higher Model 1 σ HPF1 values are driven by correlated stellar activity signals as both GP models are able to remove this additional scatter.Furthermore, Models 2 and 3, which employ GPs with different kernel choices, yield nearly identical planetary parameters.This provides evidence that the GP models are identifying similar correlated signals and are robust against kernel selection, which can affect the inferred planetary properties (e.g., Benatti et al. 2021).The GP models also decrease the rms of the RV residuals.The residual rms values are 178, 109, and 99 m s −1 for Models 1, 2, and 3, respectively.Altogether the lower overall residuals, the reduced σ HPF1 parameter, and the consistent planetary values suggest that the GP models are effectively mitigating stellar activity contributions to the RVs of HS Psc.
Finally, the uncertainties in the RV semi-amplitudes and eccentricities for Models 2 and 3 are larger than for Model 1.The 1σ uncertainties of Models 2 and 3 increase by a factor of approximately 3 for K b and 2 for e b as compared to Model 1.The lower Model 1 uncertainties can be attributed to the large planetary amplitude, long time baseline, and high number of RV observations.However, the posterior width may not fully capture the effects of stellar activity, which are not explicitly built into the model.These larger uncertainties present in Models 2 and 3 better reflect the difficulty of detecting and appropriately characterizing a young planet candidate such as HS Psc b.
Altogether we find that a joint GP model better predicts the HPF RVs compared to a Keplerian-only model.However, GP regression can be susceptible to overfitting (e.g., Aigrain & Foreman-Mackey 2022;Blunt et al. 2023), leading to systematic biases that affect the interpretation of a planetary signal.This is particularly true for young systems where stellar activity signals are large.
To assess the possibility of Models 2 and 3 overfitting the data, we apply a cross-validation, or "train-testsplit", test as described in Blunt et al. (2023).Crossvalidation is a procedure designed to test the predictive performance of a model (Gelfand et al. 1992;Aigrain & Foreman-Mackey 2022).As implemented in Blunt et al. (2023), a validation sample is first constructed by splitting the data into training and test sets, which are comprised of 80% and 20% of the data, respectively.The model is then conditioned on the training set and the best-fit model predictions at the times of the test set are compared to the test set values.If the predictions produce similar residuals for the withheld data as with the conditioned data, then the model is predictive and can be considered robust against overfitting.
We randomly withhold 17 RVs out of the 83 observations (amounting to 20% of the data) and conduct this cross-validation test for both Models 2 and 3. Figure 4 displays the training and test data sets, best-fit model predictions, and fit residuals from this cross-validation test for Models 2 and 3.For both models, the residuals of the test set have consistent structure and spread with the training set residuals.This test is repeated twice more with different data split variations.The results in all three instances are consistent.We interpret the similarity between the residuals of the training and test sets as evidence supporting the predictive nature of the best-fit models, indicating that Models 2 and 3 do not overfit the RV data.
To determine which model is favored by the data, we further evaluate comparison metrics between Models 1, 2, and 3. Equation 3 reports different model selection criteria for each of the three models, including the BIC, both the uncorrected and corrected Akaike Information Criterion (AIC and AIC c ; Akaike 1998; Tran et al. 2022) and the Akaike weights (Akaike 1981;Burnham & Anderson 2004;Liddle 2007).The BIC comparison suggests that there is strong evidence in favor of Model 3 (∆BIC > 5 against Models 1 and 2).The Akaike weights, or relative likelihoods, also prefer Model 3 over Models 1 and 2 (92% favorability).
The consistent inferred parameters, reduced residual rms, cross-validation tests, and model selection criteria all indicate that Model 3 performs the best at robustly detecting the Keplerian signal and modeling correlated stellar activity signals.As a result, we adopt the results from Model 3, the joint Keplerian and M 5 /2 GP model, as the most appropriate values for HS Psc b.

Exclusion of a Low-Inclination Brown Dwarf or Stellar Companion
RV observations alone cannot measure the true mass of a Keplerian signal.RV surveys are thus susceptible to false-positive scenarios in which a BD or stellar companion on a face-on orbit masquerades as a planet (e.g., Wright et al. 2013).In Section 4.5, we report a minimum mass of m sin i * = 1.46 M Jup .For HS Psc b to have a true mass larger than the BD limit (≥ 13 M Jup ), its orbital inclination must be i ≤ 6.4 • .Assuming an isotropic distribution for inclinations, the a priori prob-ability of this scenario is P ≈ 0.63%.For a stellar mass companion (≥ 75 M Jup ), this probability decreases to P ≈ 0.019%.However, these probabilities only reflect the chances of a false positive planet signal for this particular star and do not take into account our broader survey size; the larger the survey, the more false positive planets should be found.Binomial statistics can be used to correct for this sample size by determining the probability of a success (planet) given an underlying rate (the false positive probability) and a sample size (104 stars, in this case).This result will assume that each star in the sample has a short-period stellar or substellar companion, so the resulting probability needs to be multiplied by the actual occurrence rate of these companions.
Combining this probability of an inclined BD, binomial statistics based on our survey, and a conservative estimate of 1% for the occurrence rate of BDs, we calculate the probability of observing at least 1 false-positive event in our sample of 104 systems from this system as P FP = 0.48%.This false-positive rate is an upper limit as the occurrence rate of close-in (a ≲ 5 au) BDs is ≲ 1% (e.g., Grether & Lineweaver 2006;Sahlmann et al. 2011;Santerne et al. 2016) and even lower at closer separations (P ≲ 100 d; Ma & Ge 2014;Ranc et al. 2015;Csizmadia et al. 2015;Kiefer et al. 2019).The true probability of a false alarm scenario is therefore likely to be < 0.48%.Moreover, we note that the stellar inclination is likely high (Section 4.3) and hot Jupiters around cool stars usually have low obliquities (e.g., Winn et al. 2010;Albrecht et al. 2022), so it would be unusual for the planet inclination to differ substantially from 90 • .Based on these estimates, we confidently exclude the possibility that HS Psc b is a BD or star on an inclined orbit.

Timescales for Tidal Circularization and High-eccentricity Migration
As a young hot Jupiter, HS Psc b could offer rare constraints on the timescales of different migration processes such as tidal circularization and high-eccentricity migration.Tidal interactions between close-in planets and their host stars lead to an exchange of angular momentum and, as a result, evolution in orbital separation and eccentricity.The characteristic timescale for this mechanism is typically short compared to the characteristic age of several Gyr for field stars.A general form for the tidal circularization timescale is given by Equation 3 of Adams & Laughlin (2006), where Q P is the tidal quality factor, generally 105 -10 6 , m p and r p are the planetary mass and radius, respectively, M * is the stellar host mass, and a is the orbital separation.This relation holds for systems where the planetary orbital period is greater than the stellar rotational period (Goldreich & Soter 1966;Jackson et al. 2008).However, HS Psc b operates in the supersynchronous rotation regime, where P orb > P rot (Ferraz-Mello et al. 2008).In this state, the assumptions in Equation 5 break down and tidal forces operate less efficiently.As a result, circularization timescales estimated by Equation 5 should be treated as lower limits (Jackson et al. 2008).
Coupled with the age of the system, this lower limit can provide us with a sense of the minimum time it should take for a planet like HS Psc b to tidally circularize if it had a high eccentricity after migrating to its current orbital distance.Adopting the median values reported in this work and assuming radii of {1.0, 1.25, 1.5} R Jup , 5 the ranges of tidal circularization timescales of HS Psc b are τ cir = {0.18-1.79,0.09-0.89,0.05-0.48}Gyr.
These estimates suggest the earliest timescale for tidal circularization is on the order of several tens of Myr.This means that if the current eccentricity of HS Psc b is zero, then if it migrated recently through higheccentricity migration, the earliest that process could have occurred was a few tens of Myr ago.It could also have migrated this way soon after formation, setting an upper limit to this process of ≈130 Myr.A low eccentricity is also consistent with disk migration early on.On the other hand, if HS Psc b currently has a modest or high eccentricity, then it should still be in the process of circularizing and high-eccentricity migration could have occurred at any point in the planet's history.
We stress that while the eccentricity posteriors of the adopted model (Model 3) prefer lower values, there is power across all eccentricities and this parameter is not well constrained.Follow-up RV observations are needed to more confidently rule out higher eccentricities and better constrain the evolutionary history and migration timescale of HS Psc b.Irrespective of the specifc timing or channel, the discovery of HS Psc b implies that HJs can form or migrate to their current locations by 130 Myr.

Future Observations
Optical RV observations can further validate and characterize HS Psc b.Modulations arising from the presence of starspots are wavelength dependent; RV amplitudes of starspot-driven variability have been shown to be higher in the optical than in the NIR by a factor of ≈2 for systems at this age (e.g., Prato et al. 2008;Mahmud et al. 2011;Crockett et al. 2012;Bailey et al. 2012;Gagné et al. 2016;Tran et al. 2021).Confirming that the RV amplitude, orbital period, and orbital phase of optical RVs are consistent with the NIR HPF observations would further support the planetary nature of HS Psc b.Additional RVs can also be used to search for other companions, such as smaller closer-in or more distant giant planets, to further inform the formation and migration of the system.Moreover, more precise optical RVs can refine the eccentricity of HS Psc b and establish whether it is actively undergoing tidal circularization.• We derive an effective temperature of T eff = 4203± 116 K, a metallicity of [Fe/H] = −0.05± 0.09 dex, a surface gravity of log g = 4.66 ± 0.03 dex, and a stellar projected rotational velocity of v sin i * = 29.7 ± 3.1 km s −1 for HS Psc.We infer a stellar mass of M * = 0.69 ± 0.07 M ⊙ and radius of R * = 0.65±0.07R ⊙ for HS Psc using the ARIADNE and isochrones packages.These values are consistent with its spectral type of K7V, stellar age, and previous literature values.
• Our HPF RVs over 4 years reveal a periodicity at P = 3.99 d.This period is not commensurate with an integer harmonic of the stellar rotation period measured from TESS photometry (P rot = 1.086 ± 0.003 d).Furthermore, these RV measurements do not correlate significantly with associated activity indicators, supporting a planetary origin for the observed signal.
• A joint Keplerian and M 5 /2 GP stellar activity model fit to the HPF RVs yields a minimum mass of m b sin i = 1.46 +0.57−0.44 M Jup , an orbital period of P b = 3.986 +0.044 −0.003 d, and a broad eccentricity constraint with a slight preference for low values.No evidence of a longer-term acceleration is evident.HS Psc b is unlikely to be a BD or star on a face-on orbit.
• As a young, close-in giant planet, HS Psc b may have undergone high-eccentricity tidal migration.
If so, we estimate a lower limit of several tens of Myr for the tidal circularization timescale of HS Psc b.The age of HS Psc places an upper limit on the migration timescale of ≈130 Myr.Disk migration is also possible if HS Psc b has a low eccentricity.A modest or high eccentricity would imply that it is still undergoing circularization.Additional high-precision RV observations will help confirm HS Psc b, refine the orbit and minimum mass, and constrain its orbital evolutionary history.
HS Psc b joins only a small handful of other hot Jupiter candidates that have both robust age constraints and planetary (minimum) mass measurements.HS Psc b is an excellent target for future observations with precision optical and IR spectrographs.If confirmed with follow-up radial velocities, HS Psc b will be one of the youngest hot Jupiters discovered to date.

ACKNOWLEDGEMENTS
The authors would like to thank Marvin Morgan, Kyle Franson, and Adam Kraus for insightful discussions on the high-eccentricity tidal migration, orbital circularization, and thermal evolution of close-in giant planets.The authors would also like to thank Chad Bender, Steven Janowiecki, Greg Zeimann, the HPF team, and all the resident astronomers and telescope operators at the HET for supporting these observations and data processing.The authors are grateful to the referee for their helpful comments, which improved the quality of this manuscript.
Q.H.T. and B.P.B. acknowledge the support from a NASA FINESST grant (80NSSC20K1554).B.P.B. acknowledges support from the National Science Foundation grant AST-1909209, NASA Exoplanet Research Program grant 20-XRP20_2-0119, and the Alfred P. Sloan Foundation.GS acknowledges support provided by NASA through the NASA Hubble Fellowship grant HST-HF2-51519.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.
These results are based on observations obtained with the Habitable-zone Planet Finder Spectrograph on the HET.The HPF team acknowledges support from NSF grants AST-1006676, AST-1126413, AST-1310885, AST-1517592, AST-1310875, ATI 2009889, ATI-2009982, AST-2108512, and the NASA Astrobiology Institute (NNA09DA76A) in the pursuit of precision radial velocities in the NIR.The HPF team also acknowledges support from the Heising-Simons Foundation via grant 2017-0494.The Hobby-Eberly Telescope is a joint project of the University of Texas at Austin, the Pennsylvania State University, Ludwig-Maximilians-Universität München, and Georg-August Universität Gottingen.The HET is named in honor of its principal benefactors, William P. Hobby and Robert E. Eberly.We acknowledge the Texas Advanced Computing Center (TACC) at The University of Texas at Austin for providing high performance computing, visualization, and storage resources that have contributed to the results reported within this paper.Computations for this research were also performed on the Pennsylvania State University's Institute for Computational and Data Sciences Advanced CyberInfrastructure (ICDS-ACI, now known as Roar), including the CyberLAMP cluster supported by NSF grant MRI1626251.
We would like to acknowledge that the HET is built on Indigenous land.Moreover, we would like to acknowledge and pay our respects to the Carrizo & Comecrudo, Coahuiltecan, Caddo, Tonkawa, Comanche, Lipan Apache, Alabama-Coushatta, Kickapoo, Tigua Pueblo, and all the American Indian and Indigenous Peoples and communities who have been or have become a part of these lands and territories in Texas, here on Turtle Island.
This paper includes data collected by the TESS mission.Funding for the TESS mission is provided by the NASA's Science Mission Directorate.This work presents results from the European Space Agency (ESA) space mission Gaia.Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC).Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA).
This research has made use of the VizieR catalogue access tool, CDS, Strasbourg, France (DOI:  Tables B2 and B3 summarize the prior choices and results of the Models 1 and 2 fit, respectively.Figures B1 and B2 display the best-fit phased RV curve and posterior distributions of the planetary parameters Models 1 and 2, respectively.See Sections 4.5.1 and 4.5.2 for details.

Figure 1 .
Figure 1.TESS Sector 17 (top) and 57 (bottom) light curves for HS Psc.The best-fit mean GP model and 1σ variance are displayed as solid blue lines and shaded regions, respectively.The posterior distributions and MAP values (dashed vertical line) for the QP GP Prot hyperparameter for TESS Sectors 17 and 57 are shown in the right panels.The GP model separately recovers rotation periods of Prot, 17 = 1.086 ± 0.003 d and Prot, 57 = 1.086 ± 0.004 d for Sectors 17 and 57, respectively.

Figure 2 .
Figure2.Generalized Lomb-Scargle periodograms of the spectral time series.Beginning with the top panel, the periodograms are calculated for the HPF RVs, the differential line width (dLW) activity indicator, the chromatic index activity indicator (CRX), the Ca II IRT indices for the Ca II 1, 2, and 3 IRT emission lines, the RV residuals after subtracting the best-fit planetary model, and the spectral window function (WF).The planet orbital period (P b = 3.99 d) and the stellar rotation period (Prot = 1.086 d) are highlighted in the blue and red vertical shaded regions, respectively.False alarm probabilities (FAPs) of 1% and 0.01%, calculated using the bootstrap method, are shown as orange dotted and dashed lines, respectively.There are five peaks above the FAP = 0.01% threshold in the periodogram of the RVs; these represent the nightly sampling alias (seen at 1.0 day in the bottom panel with the WF periodogram), the stellar rotation period, the orbital period, and likely aliases of the stellar rotation and orbital periods.

Figure 3 .Figure 4 .
Figure 3. Top: RV curve of HS Psc b phase folded to the best-fit orbital period from Model 3. The blue and orange points denote the two different observing seasons, HPF1 and HPF2, respectively, and the best-fit RV model is shown as the solid black line.The fit residuals are plotted in the lower panels.The colored error bars are nominal RV errors and the grey error bars include the systematic jitter terms.1σ and 3σ confidence intervals are plotted as grey shaded regions.The right panel displays the same best-fit orbit with the phased RV points representing median values in bins of ≈0.1 and ≈0.2 of the phase for the HPF1 and HPF2 RVs, respectively.The bottom panels display the joint posterior distributions of T 0,b , P b , K b , and e b from the Model 3 fit of HS Psc b.The diagonal panels show the marginalized distribution for each parameter.The upper right panels show the marginalized distributions for the M 5 /2 kernel parameters.
the discovery and characterization of a young hot Jupiter candidate orbiting a member of the 133 +15 −20 Myr AB Dor moving group as part of the Epoch of Giant Planet Migration planet search program.Below, we summarize our main conclusions: • We obtained 83 NIR RVs of HS Psc with HPF at the HET.Using SERVAL-based leastsquares matching and SpecMatch-Emp-based empirical spectral matching algorithms, we extract the relative RVs, dLWs, CRXs, and the Ca II infrared triplet line indices.

Table 1 .
Properties of HS Psc

Table 3 .
Metrics for best-fit model comparison.

Table B3 .
Parameter Priors and Posteriors from Keplerian and Quasi-Periodic GP Fit (Model 2) to HPF RVs of HS Psc.