Migration and Evolution of giant ExoPlanets (MEEP). I. Nine Newly Confirmed Hot Jupiters from the TESS Mission

Hot Jupiters were many of the first exoplanets discovered in the 1990s, but in the decades since their discovery the mysteries surrounding their origins have remained. Here we present nine new hot Jupiters (TOI-1855 b, TOI-2107 b, TOI-2368 b, TOI-3321 b, TOI-3894 b, TOI-3919 b, TOI-4153 b, TOI-5232 b, and TOI-5301 b) discovered by NASA’s TESS mission and confirmed using ground-based imaging and spectroscopy. These discoveries are the first in a series of papers named the Migration and Evolution of giant ExoPlanets survey and are part of an ongoing effort to build a complete sample of hot Jupiters orbiting FGK stars, with a limiting Gaia G-band magnitude of 12.5. This effort aims to use homogeneous detection and analysis techniques to generate a set of precisely measured stellar and planetary properties that is ripe for statistical analysis. The nine planets presented in this work occupy a range of masses (0.55M J < MP < 3.88M J) and sizes (0.967R J < RP < 1.438R J) and orbit stars that have an effective temperature in the range of 5360 K < T eff < 6860 K with Gaia G-band magnitudes ranging from 11.1 to 12.7. Two of the planets in our sample have detectable orbital eccentricity: TOI-3919 b ( e=0.259−0.036+0.033 ) and TOI-5301 b ( e=0.33−0.10+0.11 ). These eccentric planets join a growing sample of eccentric hot Jupiters that are consistent with high-eccentricity tidal migration, one of the three most prominent theories explaining hot Jupiter formation and evolution.


INTRODUCTION
The first exoplanet discovered around a main sequence star, 51 Pegasi b (Mayor & Queloz 1995), immediately HJs.The rapid discovery of HJs has prompted questions about how these giant planets, with orbits far more compact than any in our Solar System, came to be.Several theories have been proposed to explain the preponderance of short-period giant planets, which we group into three categories: in situ formation (Batygin et al. 2016), gas-disk migration (Goldreich & Tremaine 1980;Lin & Papaloizou 1986), and high-eccentricity tidal migration (e.g., Kozai 1962;Lidov 1962;Rasio & Ford 1996;Naoz 2016).
The plausibility of in situ formation, the idea that giant planets can form in their present-day orbital configurations and not require large-scale migration, has been debated throughout the last two decades (Rafikov 2005(Rafikov , 2006; Lee et al. 2014;Boley et al. 2016;Poon et al. 2021).The majority of the currently accepted literature suggests that in situ formation is unlikely to be the dominant mechanism to form HJs, as it would require a rapid build-up of ∼ 10 Earth masses (M E ) of solids in a region of the protoplanetary disk where the feeding zones are tiny (Dawson & Johnson 2018).Additionally, the cores that form would not likely be able to accrete enough gas to form a HJ before the dissipation of the gas disk (Lee et al. 2014;Dawson & Johnson 2018).Batygin et al. (2016) argued that rapid gas accretion onto planets with masses between 10 and 20 M E is possible, but that HJs that form in situ should be accompanied by low-mass planets on short periods, a scenario which has proven to be rare (Huang et al. 2016).
If in situ formation is not the dominant mechanism to produce HJs, then there must be some mechanisms for giant planets to undergo large-scale migration from orbits beyond ∼ 3 au to less than ∼ 0.1 au.The currently favored pathways for this large-scale migration are gas-disk migration, where the planet transfers its orbital energy and angular momentum to the protoplanetary disk (Goldreich & Tremaine 1980;Lin & Papaloizou 1986), and high-eccentricity tidal migration, where the planet exchanges angular momentum with another planet or star, exciting the planet's eccentricity to be later circularized on a much smaller orbit as a consequence of tidal interactions with the host star (Rasio & Ford 1996;Fabrycky & Tremaine 2007).Each of these mechanisms is expected to generate different observable outcomes in planetary systems that have undergone migration.Gas-disk migration does not typically excite the orbital eccentricities or misalign the orbits of the migrating planets, whereas planets that have undergone high-eccentricity tidal migration may still have observable eccentric or misaligned orbits that are remnants of their migration (Dawson & Johnson 2018).Additionally, high-eccentricity tidal migration typically destabi-lizes the orbits of nearby planets, leading to the preponderance of isolated HJs (Huang et al. 2016), while gas-disk migration is more quiescent, enabling nearby companions to survive (Becker et al. 2015;Vanderburg et al. 2017).
This article is the first in connection with a survey intended to construct a magnitude-limited, complete sample of hot and warm Jupiters with precisely measured masses, radii, semi-major axes, and eccentricities.This set of parameters will be useful for future efforts to understand the HJ population and constrain the migration pathways of giant planets.Our survey, which we name the Migration and Evolution of giant ExoPlanets (MEEP) survey, is joining other efforts (Rodriguez et al. 2019(Rodriguez et al. , 2021;;Ikwut-Ukwa et al. 2022;Yee et al. 2022Yee et al. , 2023;;Rodriguez et al. 2023) to build such a catalog of parameters using data from the Transiting Exoplanet Survey Satellite (TESS ; Ricker et al. 2015) and homogeneous analysis techniques, to ensure that statistical studies can be performed on the sample with a good understanding of the selection biases.This sample will be constrained by the host star brightness to ensure that radial velocity (RV) observations yield a precise mass and eccentricity, and so only planets with host star Gaia Gband magnitudes brighter than 12.5 will be included.Yee et al. (2021) finds that there are likely roughly 300-400 discoverable transiting HJs orbiting FGK stars brighter than G = 12.5 mag and, while < 300 of these have been confirmed so far, TESS should be able to detect nearly all of the remaining undiscovered HJs.With a complete sample of HJs orbiting bright FGK stars, we will be able to confidently explore the distribution of planetary and stellar parameters of HJ-hosting systems, nailing down the occurrence rate of HJs around Sun-like stars and probing the likelihood of each formation and evolution mechanism by comparing population synthesis models to the observed distribution of orbital eccentricity and orbital period.Ikwut-Ukwa et al. (2022) found tentative evidence for multiple, separate distributions in mass-period space in the incomplete HJ population.A complete population of HJs could enable follow-up investigations of this and other tentative trends and connect them to the physics of planet formation and evolution.
In this article, we present nine new HJs discovered by TESS , eight of which orbit an FGK star brighter than G = 12.5 mag.These planets were first detected via transit photometry by TESS and classified as TESS Objects of Interest (TOIs) before being confirmed by ground-based photometry, radial velocities, and highresolution imaging.In §2, we discuss these observations in detail.In §3, we discuss the use of Markov Chain Monte Carlo to globally fit each planetary system and obtain relevant stellar and planetary parameters.In §4, we present the results of these global fits and discuss the future of our survey.Finally, in §5, we summarize the key arguments made in this article.

OBSERVATIONS
We used a combination of space-based and groundbased photometric, high-resolution imaging, and spectroscopic observations to characterize each planetary system and rule out possible false-positive scenarios, such as eclipsing binaries, blended nearby eclipsing binaries, and stellar activity.Table 1 shows the relevant parameters of each system from archival observations and ground-based spectroscopy.

TESS Photometry
All nine of the planets confirmed in this work were first classified as candidates after transits were detected by TESS .The TESS spacecraft observes 24 • × 96 • sectors of the sky for 27 days at a time, with a pixel scale of 21 arcseconds per pixel.During TESS 's Prime Mission (Sectors 1 -26), more than 200,000 bright stars were selected for 2-minute cadence, high-precision photometry, while the remainder of the stars in the fields were observed in the 30-minute cadence full frame images.In TESS 's first extended mission (Sectors 27 -55), 2-minute cadence photometry was obtained for ∼15,000 stars per sector, and the cadence of the full frame images was reduced to 10 minutes.In addition, approximately 600 targets per sector received 20-second cadence data.In the ongoing second extended mission (Sectors 56 -97), the full frame images have 200-second sampling, ∼8,000 targets per sector will be selected for 2-minute cadence, and ∼2,000 targets per sector will be selected for 20-second cadence.
The systems studied in this work were observed in at least two sectors each, between Sector 8 and Sector 63.These observations were collected in 2-minute cadence, 10-minute cadence, and 30-minute cadence.The TESS data were reduced using both the TESS Science Processing Operations Center (SPOC) Pipeline (Jenkins et al. 2016;Caldwell et al. 2020) located at NASA Ames Research Center and the MIT Quick-Look Pipeline (QLP; Huang et al. (2020a,b); Kunimoto et al. (2021)).Both pipelines account for contamination from other stars in the pixel mask.Where available, the SPOC lightcurves were chosen over the QLP lightcurves.See Table 2 for a full list of the sectors and cadences in which each system was observed and the pipelines used to process the data.In the table, SPOC lightcurves generated from the fullframe images have their source listed as TESS-SPOC, as opposed to SPOC, which refers to the SPOC lightcurves from the 2-minute data.Both the SPOC and QLP pipelines search for and identify transit-like signals, which are then vetted by the TESS Science Office.The SPOC conducted its transit searches with an adaptive, noise-compensating matched filter (Jenkins et al. 2002(Jenkins et al. , 2010(Jenkins et al. , 2020)), producing a threshold crossing event for which an initial limb-darkened transit model was fitted (Li et al. 2019) and a suite of diagnostic tests were conducted to help make or break the planetary nature of the signal (Twicken et al. 2018).Those that survive this process are identified as candidates, labeled TESS Objects of Interest (TOIs; Guerrero et al. 2021), and are released to the public.The lightcurves were downloaded using the lightkurve1 package (Lightkurve Collaboration et al. 2018), which accesses the Mikulski Archive for Space Telescopes (MAST2 ).The selected lightcurves were then flattened using keplersplinev23 , which fits a spline to any out-of-transit stellar variability in the lightcurve, which can then be divided out to flatten the lightcurve without significantly affecting the transit (Vanderburg & Johnson 2014).In order to optimally choose the spline breakpoints, we adopted the methodology from Shallue & Vanderburg (2018) and used the choosekeplersplinev2 feature, which chooses the breakpoint spacing that minimizes the Bayesian Information Criterion (BIC; Schwarz 1978).We then chop the flattened lightcurves to ensure that only a baseline of one transit duration (T 14 ) is kept on either side of each transit to reduce the computational cost of using 27-day, nearly uninterrupted, lightcurves in our global fits.

Ground-based Photometry
To confirm and better characterize each system, we collected ground-based photometric follow-up of all nine targets through the TESS Follow-up Observing Program (TFOP; Collins et al. 2018) Sub Group 1 (SG1; Collins 2019).This ground-based follow-up serves multiple purposes: to rule out contamination from nearby eclipsing binaries that were blended in TESS 's larger pixels, to collect observations in multiple different bandpasses and rule out that the target star itself is an eclipsing binary, and to extend the observation baselines to provide better constraints on the orbital period and transit epoch.The follow-up photometry from SG1 is included in our global analysis (see §3) and shown in Table 3 and Figures 6-14.Additionally, for more detailed information on the follow-up observations collected for each of these planets, see Table 3.All 40 follow-up lightcurves are available to the public and can be downloaded from ExoFOP4 .
The Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) generated twelve of the follow-up lightcurves of our targets using observations from the following LCOGT sites: Cerro Tololo Inter-American Observatory (CTIO), the Siding Spring Observatory (SSO), the South African Astronomical Observatory (SAAO), and the McDonald Observatory (McD).Additional ground-based observations of our targets were collected using the following facilities: the National The data calibration and reduction for each followup observation (with the exception of the observations from MuSCAT and MuSCAT2; see §2.2.1) were done using AstroImageJ (AIJ) (Collins et al. 2017), a tool designed to enable point-and-click photometry and data reduction.The exposures from the observation are first calibrated using flat field images, dark images, and bias images by the observer or facility that collected the ob-servations.The general procedure5 that is then used to generate a reduced lightcurve for each target is as follows: • Generate calibrated images using AIJ's CCD Data Processor tool, median combining the biases, darks, and flats, before subtracting the master bias and dark and dividing the master flat.
• Locate the target star in the first image and develop an aperture that has a diameter of roughly twice the full width at half maximum (FWHM) of the target.The background annulus should have roughly the same area as the target aperture.
• Place apertures on the target star and at least 5-10 comparison stars of similar brightness.The number of comparison stars varies significantly depending on the field of view of the telescope and the density of stars in the surrounding area of the sky.
• Assess the quality of the data and search for relevant parameters that change significantly throughout the observation and correlate with changes in the relative flux of the host star.Likely culprits usually include the airmass, the target star's X and Y pixel location (due to imperfect guiding), and the width or FWHM of the target source (due to evolving weather conditions).
• Use AIJ's built-in transit-fitting tool to fit the transit with detrending enabled.Choose the detrending parameters from the previous step that minimize the Bayesian Information Criterion.
• Normalize the relative flux and save a lightcurve file containing the Barycentric Julian Date in Barycentric Dynamical Time (BJD TDB ), normalized relative flux, relative flux error, and the selected detrending parameters.
This reduced lightcurve is then ready to be included in our global fit (see §3).

MuSCAT and MuSCAT2 observations of TOI-1855
Two of the follow-up light curves that we obtained for TOI-1855 b were reduced using a reduction strategy different than the one presented above.
We observed TOI-1855 b using MuSCAT, a simultaneous multi-band photometer installed on the 188cm telescope of the National Astronomical Observatory of Japan (NAOJ) in Okayama, Japan (Narita et al. 2015).It has three 1k CCDs each with 6.1 ′ × 6.1 ′ field-of-view (FOV) and 0.358 ′′ pix −1 pixel scale, enabling simultaneous photometry in the g (400-550 nm), r (550-700 nm), and z s (820-920 nm) bands.The data reduction and differential photometry were performed using the pipeline described in Fukui et al. (2011).We optimized both the aperture radii and the set of comparison stars by minimizing the dispersion of the resulting relative light curves.However, there are only two useful comparison stars within the FOV, one of which is a faint star 12.8 ′′ from the target.Using the brighter comparison star and uncontaminated aperture radius of 7.2 ′′ yielded the optimum light curves with a typical rms of 0.5 ppt/10 min across all bands.
We additionally observed TOI-1855 b with the 1.5m Telescopio Carlos Sánchez (TCS) instrument MuS-CAT2 (Narita et al. 2019) from Teide Observatory, Tenerife, Spain.MuSCAT2 is another multi-band imager, equipped with four cameras with a field of view of 7.4×7.4arcmin 2 (pixel scale of 0.44 arcsec per pixel) and is capable of obtaining simultaneous observations in g ′ , r ′ , i ′ , and z s bands.The exposure times were initially set to 8 s, 6 s, 15 s, and 12 s in g ′ , r ′ , i ′ , and z s , respectively.During the observations, the i ′ , and z s band exposure times were changed to 12 s and 9 s respectively to avoid the saturation of the target star.The dedicated MuSCAT2 pipeline (Parviainen et al. 2019) reduced the raw science frames using standard procedures (i.e., dark and flat field corrections) and computed the photometry for the stars in the field using a set of circular apertures.Finally, the pipeline obtained the best light curve by fitting a transit model that includes instrumental effect noise components and choosing the aperture that delivered the lowest scatter in the time series.

Spectroscopy
In order to rule out false-positive scenarios and measure the mass and orbital eccentricity of each planet, we collected spectroscopic observations of each host star to measure the radial velocity (RV) shift imparted on the star by the orbiting planet.The spectra that we used in our final fits of the nine systems in this sample were collected by the Tillinghast Reflector Échelle Spectrograph (TRES) on the 1.5-meter Tillinghast Reflector at the Fred Lawrence Whipple Observatory (Fűrész 2008), the CHIRON spectrograph on the SMARTS 1.5-meter telescope at the Cerro Tololo Inter-American Observatory (Tokovinin et al. 2013), and the NEID spectrograph on the 3.5-meter WIYN Telescope at the Kitt Peak National Observatory (Schwab et al. 2016).The RV of each star as a function of time and orbital phase are shown in Figures 6-14.An example RV measurement is shown for each target and instrument in Table 4. Additionally, the spectroscopic metallicity measurements from TRES and CHIRON were averaged and used as a wide Gaussian prior in our EXOFASTv2 global fits, with a prior width equal to twice the standard deviation of the measured metallicities.

TRES Spectroscopy
Five of the nine host stars presented in this paper (TOI-3894, TOI-3919, TOI-4153, TOI-5232, and TOI-5301) were observed using the Tillinghast Reflector Échelle Spectrograph (TRES) to obtain precise radial velocities.TRES is a fiber-fed, high-resolution échelle spectrograph mounted on a 1.5-meter telescope with a resolving power (R) of 44,000.The reduction of the spectra was done following the work of Buchhave et al. (2010).Multi-order velocities were derived by crosscorrelating the spectrum with the highest signal-to-noise per resolution element against each observation, order by order.To measure the metallicity, effective temperature, surface gravity, and projected rotational velocity of the star, we analyzed the spectra with the Stellar Parameter Classification (SPC) tool (Buchhave et al. 2012).Finally, we calculated the bisector spans of the TRES spectra following the work of Torres et al. (2007) to hunt for correlations between the bisector spans and RVs, an indication that the RVs are being induced by an eclipsing binary.For each of the five targets with TRES spectra, we found no evidence of correlation between the bisector spans and the RVs.

CHIRON Spectroscopy
The other four host stars in this work (TOI-1855, TOI-2107, TOI-2368, and TOI-3321) were observed using CHIRON, a fiber-fed, high-resolution échelle spectrograph.For the majority of the observations in this paper, we used an image slicer that offers a resolving power of ∼ 80,000.The exception to this is TOI-2107, of which we observed 23 spectra in the fiber mode (R ∼ 25,000) before switching to the slicer mode for the remaining 13 observations.In both cases, we bracketed each of the spectra we collected with ThAr calibration spectra.To obtain estimates of the metallicity, effective temperature, surface gravity, and projected rotational velocity, each spectrum is matched to an interpolated catalog of ∼ 10,000 spectra classified by SPC (Buchhave et al. 2012).RVs were derived by the least-squares deconvolution (Donati et al. 1997;Zhou et al. 2021) of the observed spectra against non-rotating synthetic templates generated using ATLAS9 model atmospheres (Kurucz 1992), before fitting the resulting line profile with a rotational broadening kernel prescribed by Gray (2005).

NEID Spectroscopy
We obtained six observations of TOI-3894 with the NEID spectrograph, between UT 2022 April 29 and UT 2022 June 09.NEID is a highly-stabilized, fiber-fed optical spectrograph (Schwab et al. 2016;Halverson et al. 2016) on the WIYN 3.5m telescope at Kitt Peak National Observatory (KPNO).We used NEID in highresolution (HR) mode (R ∼ 110,000), and the data were reduced using versions 1.1.2-1.1.4 of the NEID Data Reduction Pipeline (DRP).The DRP extracts precise RVs by cross-correlating the observed spectra with a stellar line mask (Baranne et al. 1996;Pepe et al. 2002).Each observation was made with exposure times between 300s and 420s, achieving a typical RV precision of ≈ 10 m s −1 .On the nights of 9 and 10 May 2022, the NEID data were affected by a failure of the Fabry-Perot etalon laser, which is used to calibrate nightly offsets.We used observations of standard stars to determine those offsets (35.8 ± 1.2 m s −1 and 31.5 ± 1.3 m s −1 respectively, see Yee et al. 2023) and correct the NEID RVs taken on the affected nights.

Supplementary Spectroscopy
Three additional spectra of TOI-1855 were collected by the Skalnaté Pleso Observatory in the High Tatras (Slovak Republic), using the 1.3 m f/8.36 Astelco Alt-azimuthal Nasmyth-Cassegrain reflecting telescope, equipped with a fiber-fed échelle spectrograph of MU-SICOS design (Baudrand & Bohm 1992).The spectra were recorded using an Andor iKon-L DZ936N-BV CCD camera with a 2048 × 2048 array, 13.5 µm square pixels, 2.9e-readout noise, and a gain close to unity.The spectral range of the instrument is 4250-7375 Å (56 échelle orders) with a maximum resolving power of R = 38,000.The raw data were then reduced using IRAF package tasks, Linux shell scripts, and FORTRAN programs (Pych 2004;Pribulla et al. 2015;Garai et al. 2017).
These RVs were used to increase our confidence that TOI-1855 b is a bonafide exoplanet.Because there were only three observations from this telescope, they were not included in our final EXOFASTv2 global fit.However, these observations were consistent with the 204 +16 −15 ms −1 RV semi-amplitude measured by CHIRON.

High-resolution Imaging
Close stellar companions (bound or line of sight) can confound exoplanet discoveries in a number of ways.The detected transit signal might be a false positive due to a nearby eclipsing binary and even real planet discoveries will yield incorrect stellar and exoplanet parameters if a close companion exists and is unaccounted for (Ciardi et al. 2015;Furlan & Howell 2017).Additionally, the presence of a close companion star leads to the non-detection of small planets residing in the same exoplanetary system (Lester et al. 2021).Given that nearly one-half of solar-like stars are in binary or multi-star systems (Matson et al. 2018), high-resolution imaging (HRI) provides crucial information toward our understanding of exoplanet formation, dynamics, and evolution (Howell et al. 2021).
In order to understand and properly account for the light contribution of nearby field stars and stellar companions, we observed all nine of the host stars in our sample using adaptive optics (AO) and speckle imaging.Both of these techniques are capable of high angular resolution and allow us to place upper limits on the brightness of a nearby star and therefore the potential of any nearby star to contaminate the signal in the photometric and spectroscopic data.
Twenty AO and speckle observations from eight different instruments were used to find or rule out nearby contaminants to the host stars in this paper.The results of these observations are presented in the following sections, organized by the instrument used.An example of a non-detection (TOI-3919) and positive detection (TOI-5232) of a nearby stellar companion is presented in Figure 1.The remainder of the high-resolution images and sensitivity curves are available on ExoFOP and are described in §2.5.1 and §2.5.2 and presented in Table 5.

Adaptive Optics
Three AO instruments were employed to observe the stars studied in this work: the ShARCS camera on the Shane 3-meter telescope at Lick Observatory (Kupke et al. 2012;Gavel et al. 2014;McGurk et al. 2014;Dressing et al. in preparation), the PHARO instrument on the Hale 5-meter telescope at Palomar Observatory (Hayward et al. 2001), and the NIRC2 instrument on the Keck II 10-meter telescope at Keck Observatory (Wizinowich et al. 2000).
The ShARCS observations of TOI-1855 were taken with the Shane adaptive optics (AO) system in natural guide star mode.We collected sequences of observations using a K s filter (λ 0 = 2.150 µm, ∆λ = 0.320 µm) and a J filter (λ 0 = 1.238 µm, ∆λ = 0.271 µm).We reduced the data using the publicly available SImMER pipeline  (Savel et al. 2020(Savel et al. , 2022)). 6This observation ruled out any bright stellar companions other than TIC 81247738, the source with a separation of 12.352 ′′ already identified by Gaia.
Additionally, we observed TOI-1855, TOI-3894, and TOI-3919 with the PHARO instrument (Hayward et al. 2001) on the Palomar Hale (5m) behind the P3K natural guide star AO system (Dekany et al. 2013) in the narrowband Br-γ filter (λ 0 = 2.1686; ∆λ = 0.0326 µm).The PHARO pixel scale is 0.025 ′′ per pixel.A standard 5-point quincunx dither pattern with steps of 5 ′′ † Detection refers to a positive detection of a star within the field of view of the AO or speckle instrument, subject to the maximum contrast possible with the instrument in question.TOI-1855's companion was outside of the field of view of each instrument that observed it.
was repeated twice with each repeat separated by 0.5 ′′ .The final resolutions of the combined dithers were determined from the full-width half-maximum (FWHM) of the point spread functions: 0.09 ′′ .Observations of TOI-5232 were made with the NIRC2 instrument behind the natural guide star AO system (Wizinowich et al. 2000) in the standard 3-point dither pattern that is used with NIRC2 to avoid the left lower quadrant of the detector which is typically noisier than the other three quadrants.The dither pattern step size was 3 ′′ and was repeated twice, with each dither offset from the previous dither by 0.5 ′′ .NIRC2 was used in the narrow-angle mode with a full field of view of ∼ 10 ′′ and a pixel scale of approximately 0.0099442 ′′ per pixel.The NIRC2 observations were made in the Kcont filter (λ o = 2.2706; ∆λ = 0.0296 µm) and the Hcont filter (λ o = 1.5804; ∆λ = 0.0232 µm).The final resolutions of the combined dithers were once again determined from the FWHM of the point spread functions: 0.067 ′′ .
For both the PHARO and the NIRC2 images, the sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 20 • at separations of integer multiples of the central source's FWHM (Furlan et al. 2017).The brightness of each injected source was scaled until standard aperture photometry detected it with 5σ significance.The final 5σ limit at each separation was determined from the average of all of the determined limits at that separation and the uncertainty on the limit was set by the rms dispersion of the azimuthal slices at a given radial distance.Only TOI-5232 shows a close companion -also detected in the optical speckle data and by Gaia.

Speckle Imaging
In addition to our AO observations, five speckle instruments were used to observe the stars in our sample: the speckle polarimeter on the 2.5-m telescope at the Caucasian Mountain Observatory of the Sternberg Astronomical Institute (SAI) at Lomonosov Moscow State University, the HRCam instrument on the 4.1meter Southern Astrophysical Research (SOAR) telescope (Tokovinin 2018) at CTIO, the NN-EXPLORE Exoplanet Stellar Speckle Imager (NESSI; Scott et al. 2018) on the WIYN 3.5-meter telescope at Kitt Peak Observatory, and the Zorro and 'Alopeke instruments on the Gemini-South and Gemini-North 8-meter telescopes, respectively (Scott et al. 2021).
The HRCam observations were made of TOI-1855, TOI-2107, TOI-2368, and TOI-3321 in the Cousins Iband, a similar visible bandpass as TESS.The four observations were typically sensitive to a 5-magnitude fainter star at an angular distance of 1 arcsec from the target.More details of the speckle observations from the SOAR TESS survey are available in Ziegler et al. (2020).A 4.6 magnitude fainter star was detected at a separation of 0.69 ′′ from TOI-3321.No nearby stars were detected within 3 ′′ of TOI-1855, TOI-2107, or TOI-2368 in the SOAR observations.TOI-1855, TOI-3919, TOI-4153, and TOI-5301 were observed with the SAI speckle polarimeter.The speckle polarimeter used an Electron Multiplying CCD (EM-CCD) Andor iXon 897 as its main detector prior to August 2022 (Safonov et al. 2017).Since then, it has used the high-speed, low-noise CMOS detector Hamamatsu ORCA-quest (Strakhov et al. 2023).The atmospheric dispersion compensator was active in all cases.Observations were carried out in the I c band, where the respective angular resolution is 0.083 ′′ .No nearby stars were detected in any of the SAI images.
The NN-EXPLORE Exoplanet Stellar Speckle Imager (Scott et al. 2018) was used to obtain speckle imaging of TOI-1855 and TOI-3894.NESSI is a dual-channel speckle imager at the WIYN 3.5 m telescope on Kitt Peak, Arizona.TOI-1855 was observed with two filters having central wavelengths of λ c = 562 and 832 nm.It was then re-observed using only the 832 nm filter.TOI-3894 was observed using only the 832 nm filter.For each observation, the star was observed by taking 9 sets of 1000 40 ms exposures in a 256×256 pixel section of the EMCCD for a 4.6 × 4.6 arcsecond field of view centered on the target.A nearby single star was observed in close temporal proximity to each science target using a set of 1000 exposures to serve as a measure of the point spread function.The speckle pipeline data reduction followed the description given in Howell et al. (2011).No companion sources were detected for either TOI-1855 or TOI-3894.
Finally, TOI-3321 and TOI-5232 were observed using the Zorro and 'Alopeke speckle instruments mounted on the Gemini South/North 8-m telescopes (Scott et al. 2021).Zorro and 'Alopeke both provide simultaneous speckle imaging in two bands (562 nm and 832 nm) with output data products including Fourier analysis and reconstructed images with robust contrast limits on companion detections.Seven sets of 1000 X 0.06 second images were obtained for TOI-3321 and ten sets of 1000 X 0.06 second images were obtained for TOI-5232.Each of the observations was processed using our standard reduction pipeline (see Howell et al. 2011).Figure 2 shows our final 5σ contrast curve and the reconstructed speckle image from the TOI-3321 observation.We find that TOI-3321 has a close companion with a separation of 0.663 ′′ at a Position Angle of 79.1 degrees and is only detected in the 832 nm filter.The delta magnitude in 832 nm filter, that is the magnitude of the companion relative to the primary star at 832 nm, is 4.89 magnitudes.TOI-5232 also reveals a companion star at a separation of 1.23 ′′ and Position Angle of 5.2 degrees.This companion was detected in both filters and has delta magnitudes of 5.10 ± 0.5 (562 nm) and 4.64 ± 0.48 (832 nm) with both delta magnitude values being fairly uncertain due to the large separation (>1.0 arcsec) for which the speckles are decorrelated.At the distance of TOI-3321 (d = 286 pc) and TOI-5232 (d = 609 pc) the angular limits from the diffraction limit (20 mas) out to 1.2 arcsec correspond to spatial limits of 5.7 to 343 AU and 12 to 727 AU respectively.

Possible Stellar Companions
Of the nine planet-hosting stars in our sample, three (TOI-1855, TOI-3321, and TOI-5232) have nearby stars that are likely associated with our targets.The likely companions of TOI-1855 and TOI-5232 were both observed by Gaia and have distances and proper motions that are consistent with their probable host stars.TOI-3321's nearby star is faint, with ∆mag ∼ 4.6−4.9, and is fairly close to TOI-3321, with a separation of ∼ 0.679 ′′ .This source doesn't appear to be associated with any of the nearby Gaia DR3 (Gaia Collaboration et al. 2023) detections.The Gaia Renormalized Unit Weight Error (RUWE), a metric that is used to determine whether or not the astrometric solution of a star is well-behaved for a single source, of this target is 1.283, which is less than the upper threshold of RUWE = 1.4,indicating that the Gaia solution is a well-behaved solution for a single source.A speckle image of TOI-3321 and the corresponding sensitivity curve from Gemini's Zorro speckle imager are shown in Figure 2, illustrating the detection of TOI-3321's faint companion.
TOI-1855 and its tentative stellar companion have a ∆mag of 1.94 in the Brγ filter, which is small enough to influence TOI-1855's lightcurve; however, the 12.352 ′′ separation is distant enough that the Bergeron, MuS-CAT, MuSCAT2, LCO-CTIO, and SOAR observations were all able to separate the stars and create uncontaminated lightcurves, confirming that the transit signal is not the consequence of a blended eclipsing binary.Additionally, the companion is distant enough that the spectral energy distribution of TOI-1855 (see Figure 6) is unaffected by the light from the companion.TOI-3321's tentative stellar companion is fainter than TOI-3321 by 4.89 magnitudes in Gemini's 832 nm filter, which is too faint to significantly influence our lightcurve.Similarly, TOI-5232's tentative stellar companion is fainter than TOI-5232 by 4.64 magnitudes in the same filter and is also too faint to significantly influence the photometry.

Archival Multi-Band Observations
In order to constrain the spectral energy distribution (SED) of the stars in our sample, we queried VizieR (Ochsenbein et al. 2000) to obtain broadband observations from Gaia EDR3 (Gaia Collaboration et al. 2021), 2MASS (Cutri et al. 2003;Skrutskie et al. 2006), and WISE (Wright et al. 2010;Cutri et al. 2012).All of our targets were observed in the Gaia G, Bp, Rp, 2MASS J, H, Ks, and WISE W1, W2, and W3 bandpasses, which provide wavelength coverage from ∼ 0.3 − 12µm.In addition, astrometric measurements from Gaia EDR3 were obtained in order to place a Bayesian prior on the parallax of each target (see §3.1 for more details).These archival observations were then used in our EXOFASTv2 global fits to determine relevant stellar parameters (see §3).These photometric and astrometric parameters are reported in Table 1.

EXOFASTV2 GLOBAL FITS
We used the publicly available code EXOFASTv27 (Eastman et al. 2019) to determine all of the system parameters for all nine of the targets in this work.EXOFASTv2 is a Differential Evolution Markov Chain Monte Carlo (MCMC) code, written in IDL, designed to globally fit both the star(s) and planet(s) in each system together, ensuring a self-consistent set of parameters describing each system.These fits combine the TESS data (discussed in §2.1), follow-up transit observations from the ground (discussed in Section 2.2) and the archival multi-band photometric observations (discussed in Section 2.6) to generate a SED model using MESA Isochrones and Stellar Tracks (MIST) (Paxton et al. 2011(Paxton et al. , 2013) and a transit model simultaneously at each step of the MCMC.Additional details about the EXOFASTv2 modeling suite can be found in Eastman et al. (2019).
To test for long-term radial velocity trends, we ran an initial fit for each system, allowing for a linear trend in the radial velocities.If the fit found that the slope of this trend was more than 1σ away from zero, we fit for the linear trend in the final fit; otherwise, the slope was fixed at zero.There were only two cases in which a tentative long-term radial velocity trend was found: TOI-2368 and TOI-3919.In all of the fits, we allowed the planet's orbital eccentricity to be a free parameter.In doing this, we provide a measurement and estimated uncertainty for the orbital eccentricity of each of the systems in this study, as other derived parameters depend on eccentricity.However, because of the Lucy-Sweeney bias (Lucy & Sweeney 1971), eccentricities that are smaller than 2.45σ above zero should be considered consistent with a circular orbit (Eastman et al. 2019).
The resulting median and 68% confidence interval of all relevant planetary and stellar parameters for each system are presented in Table 6.Two of the systems that we fit had bimodal posterior distributions in stellar mass and age.These systems are discussed in further detail in §4.2.

NOTES:
The priors listed at the top of the table are labeled as G[mean, standard deviation] if they are Gaussian priors and U [lower limit, upper limit] if they are uniform priors.† TOI-1855's transit is grazing, and as a consequence, the radius measurement is highly uncertain.‡ Initial metallicity represents the metallicity of the star at formation.

Priors and Constraints
In all nine fits, we adopted Gaussian priors for the Gaia parallax measurements from Early Data Release 3 (Gaia Collaboration et al. 2021).We also placed wide Gaussian priors on the stellar metallicity, using averaged metallicity measurements from the TRES and CHIRON spectra, with a prior width equal to twice the standard deviation of the spectroscopic metallicities.Finally, we account for the possible contamination of the TESS target pixel by fitting a dilution term and placing a Gaussian prior centered at 0% with a standard deviation that is 10% of the contamination ratio from the TESS Input Catalog (TIC) v8.2 (Stassun et al. 2018(Stassun et al. , 2019)).This contamination is already corrected by both the QLP and SPOC pipelines (see §2.1), so this prior acts as a conservative assumption that the contamination correction had a precision better than 10%.Finally, we placed an upper limit on the V-band extinction of each target using the Schlegel et al. (1998) and Schlafly & Finkbeiner (2011) dustmaps.Each of these priors is included at the top of Table 6.
In addition to the supplied Gaussian priors, we adopt starting values for the stellar mass, stellar radius, and stellar effective temperature from the TIC.To provide a constraint on the stellar radius, we apply an upper limit on the V-band extinction from Schlafly & Finkbeiner (2011).Finally, the starting values for the transit epoch T C , orbital period P , and the ratio of planetary and stellar radii R P /R * for our initial fit are retrieved from the TESS mission catalog on ExoFOP 8 .One of our targets, TOI-1855, has a grazing transit configuration, in which the planet incompletely transits the limb of the star.As a result, the radius of this planet is poorly constrained.To account for the effectively flat radius posterior and to reduce the computational cost of the fit, we placed a generous physical upper limit on the radius of the planet, equal to 2.5R J .This upper limit is significantly larger than the largest nongrazing, transiting planet (Smalley et al. 2012), ensuring that our interpretation is not affected by the placement of a strict upper limit.The upper limit does, however, likely affect the median of the R P posterior listed in Table 6, and so for that reason, we also highlight the posterior mode, 1.38 R J , as the most likely radius of the planet, illustrated with a vertical dashed line in Figure 3.A more detailed discussion of TOI-1855 and our interpretation of the radius distribution is included in §4.1.In this first installment of the Migration and Evolution of giant ExoPlanets (MEEP) survey, we confirm nine new hot Jupiters (TOI-1855 b, TOI-2107 b, TOI-2368 b, TOI-3321 b, TOI-3894 b, TOI-3919 b, TOI-4153 b, TOI-5232 b, and TOI-5301 b) detected by TESS , orbiting FGK stars.These planets range in mass between 0.554 M J and 3.88 M J and range in size between 0.967 R J and 1.438 R J .The stars they orbit are all brighter than G = 12.7, which makes them suitable for further followup from a variety of ground-based facilities.

RESULTS AND DISCUSSION
The transit photometry, radial velocity, spectral energy distribution, and MIST evolutionary plots for each system investigated within this work are organized by system and presented in Figures 6-14.The relevant stellar and planetary parameters from our EXOFASTv2 global fits are presented in Table 6.

TOI-1855 b's grazing transit
One of our planets, TOI-1855 b, exhibits a V-shaped transit caused by the planet only grazing the limb of the star.This transit shape is similar to what is usually seen for an eclipsing binary.Because of this, we paid special attention to TOI-1855 to ensure that our interpretation of it as a transiting giant planet is correct.TOI-1855 was observed by three different speckle instruments with four different filters and by two different AO instruments in three different filters.These observations detected a stellar companion, previously identified by Gaia (TIC 81247738), at a separation of 12.352 ′′ with a TESS magnitude of 13.524 (∆T = 2.82 mag).No other stellar companions were detected by any of the six HRI instruments.While these sources are blended in the TESS observations, the Bergeron, MuSCAT, MuSCAT2, LCO-CTIO, and SOAR observations used uncontaminated apertures and confirmed the transit signal to be on target.Additionally, the TESS data are corrected for contamination by both the SPOC and QLP pipelines (see §2.1) and we leave the dilution term for the TESS lightcurves as a free parameter in our EXOFASTv2 global fits (see §3) to account for possible errors in the dilution correction.
Another indicator of a blended eclipsing binary that is visible in the TESS data is a centroid offset, where the position of the target star appears to change as the stars enter an eclipse (Batalha et al. 2010).We find from the data validation report9 for the SPOC-processed TESS Sector 50 data that the measured centroid offset is 0.408 ± 2.55 ′′ , well within the expected variance and therefore consistent with no centroid offset.Finally, we searched the TESS data for evidence of a secondary eclipse, another sign of an eclipsing binary, by phasing the TESS data on the optimal ephemeris from our EXOFASTv2 global fit and binning the phased lightcurve.We find no evidence of a secondary eclipse in either sector of data.The combination of an RV orbit, HRI non-detection, secondary eclipse non-detection, and the lack of a centroid offset imply that it is highly unlikely that TOI-1855 b's transit signal is the consequence of an eclipsing binary.
The other difficulty raised by the grazing configuration of TOI-1855 b's transit is that we are unable to precisely constrain its radius from the transit alone.In §3.1, we discussed the R P upper limit placed on TOI-1855's fit and the resulting posterior mode, 1.38 R J .In order to get a different constraint on the radius of TOI-1855 b, we used the code Forecaster10 from Chen & Kipping (2017) to generate a PDF of the radius of TOI-1855 given its mass distribution from our EXOFASTv2 global fit.From this mass-radius relation alone, we find a median planetary radius of 1.21 ± 0.21 R J , which is consistent within 1σ with the posterior mode (1.38 R J ) and the median (1.65 +0.52  −0.37 ) from the fit.Forecaster is built from a sample of planets with a wide range of isolations, but it is possible that TOI-1855 b's size could be underestimated by using only Forecaster if its radius is inflated by insolation.

Bimodal Solutions
Two of our EXOFASTv2 global fits converged on solutions that had bimodal posterior distributions in stel-lar mass and age: TOI-3894 and TOI-5301.In order to characterize each individual solution from both fits, we split the probability density functions of TOI-3894's and TOI-5301's fits, identifying a solution with a higher stellar mass and a separate solution with a lower stellar mass for each (see Figure 4).For TOI-3894, if we separate the PDF using a stellar mass of 1.145 M ⊙ (the center of the valley between the two peaks in Figure 4), the EXOFASTv2 global fit slightly favors the lower mass solution (M ⋆ = 1.089 +0.038 −0.051 M ⊙ , age = 7.2 +1.7 −1.1 Gyr) with a probability of 53.0%.In the case of TOI-5301, we split the PDF using a stellar mass of 1.375 M ⊙ and find that the higher mass solution (M ⋆ = 1.507 +0.068 −0.067 M ⊙ , age = 2.29 +0.42  −0.4 Gyr) is favored with a probability of 78.3%.Both of these uncertain fits are consequences of the stars' observed parameters indicating that they are near an evolutionary transition.To encapsulate both possibilities in both systems, we present the median and standard deviation of each parameter for the higher mass and lower mass solutions of both targets in Table 7.

Early Trends in the Hot Jupiter Population
The purpose of this survey is to generate a large, complete sample in order to statistically assess the origins of HJs.As such, it is too early to make conclusive claims that require robust statistics.Instead, we argue that there are several tentative trends arising in the present sample of HJs.
First, it is notable that seven of our nine planets have orbits that are consistent with circular within 2.45σ (Lucy & Sweeney 1971) and are thus consistent with either migration mechanism: gas-disk migration or higheccentricity tidal migration.Two of the planets in this sample, TOI-3919 b and TOI-5301 b, have significant eccentricities that cannot be explained by quiescent disk migration alone.TOI-3919's radial velocities show signs of a long-term trend that could be evidence of an outer companion and therefore a possible perturber to boost the eccentricity of TOI-3919 b and trigger higheccentricity tidal migration.Both of these planets is shown in Figure 5 and compared to several different avenues and outcomes of planet migration.These two planets join the 39 other HJs with reported eccentricities > 3σ above zero (11.2% of the HJ population) per the NASA Exoplanet Archive 11 (date accessed: 2023 October 23).While the NASA Exoplanet Archive is a heterogeneous source of exoplanet parameters and should 0.9 1.0 1.1 1.2 1.3 1.    and TOI-5301's (right) stellar mass and age posterior distributions from our EXOFASTv2 global fits, illustrating a bimodal distribution in both parameters.The cyan dashed line indicates the median value of each parameter, while the shaded region represents the 68% confidence region.The normalized probability densities were estimated using the Gaussian kernel density estimator in the Python package scipy.stats.TOI-3894's fit slightly favored the lower mass solution with a probability of 53.0% while TOI-5301's fit favored the higher mass solution with a probability of 78.3%.not be used to identify definitive trends, this is a significant minority of the population that must receive further scrutiny.Since many planets with periods less than ten days are expected to have very short tidal circularization timescales (Adams & Laughlin 2006), the fraction of HJs that once had boosted eccentricities could be much larger than 11.2%.This implies, that at the very least, high-eccentricity tidal migration is an important migration mechanism that must be considered for a non-negligible fraction of the HJ population.Rodriguez et al. (2023) and Zink & Howard (2023) both argue that the current eccentricity distribution of known HJs is consistent with high-eccentricity tidal migration being the dominant pathway for the evolution of HJs, an argument that can continue to be tested with a growing sample of HJs.
Additionally, two of the stars that we observed have long-term RV trends: TOI-2368 has a linear slope of 0.21 ± 0.11 m/s/day and TOI-3919 has a slope of −0.79 +0.31  −0.30 m/s/day.If further RV follow-up reveals a continuation of the trend consistent with a giant substellar companion on a wide orbit, then these systems would be further evidence that HJ-hosting systems are not always devoid of other planets, as RV surveys (e.g., Knutson et al. 2014;Bryan et al. 2016;Zink & Howard 2023) have argued.Zink & Howard (2023) also argue that the properties of systems hosting a HJ and a distant giant companion favor a scenario in which the HJ migrated via coplanar high-eccentricity tidal migration.
Several issues with high-eccentricity tidal migration have been raised, however, that must be addressed.For example, Socrates et al. (2012) argued that we should be able to detect a significant number of "super-eccentric" Jupiters that are actively migrating to become HJs, more than have yet been observed.Additionally, they argue that the tidal efficiency of HJs must be at least 10 times larger than Jupiter's.Further studies should address these issues and elucidate the importance of each mechanism of giant planet formation and evolution.

Lithium in TOI-5301
While vetting the spectra of the stars in our sample, we discovered a discernible lithium absorption feature in the spectra of TOI-5301.The strength of a star's lithium absorption feature serves as a valuable tool to determine stellar age (e.g., Vauclair 1973).This feature is also useful in identifying instances where a star may have ingested planets (Montalbán & Rebolo 2002), par- The region labeled "Stellar Collision" corresponds to a giant planet colliding with the star, assuming the star is 1 R⊙ in size.The region labeled "Tidal Disruption" corresponds to a Jupiter-like planet falling within the Roche limit of the Sun.Finally, the red region corresponds to the tidal circularization of a highly eccentric giant planet around a Sun-like star.The blue circles represent the EXOFASTv2 median semi-major axis and eccentricity of the planets in this work, while the gray circles represent substellar bodies with precise eccentricities and masses between 0.25 -13 MJ, obtained from the NASA Exoplanet Archive.Many of the eccentric planets, including those that fall outside of the three shaded regions, can be explained by planet-planet scattering that has not yet or will not be tidally circularized.The equations that describe each of these equations are shown in Dawson & Johnson (2018).
Is the lithium abundance of our target unusual when compared to a control sample of stars?To answer this question, we performed a lithium abundance analysis of the target by co-adding 19 TRES observations taken between 2022-09-01 and 2022-11-09.Before co-adding the spectra to build up the signal-to-noise ratio (SNR), each spectrum was continuum normalized, barycenter corrected, and Doppler corrected.Following the procedures outlined in Jeffries et al. (2023), we measured the equivalent width (EW) of the lithium absorption doublet (Li I) at 6707.8 Å.We measure an SNR=84 for the co-added spectra.Using pymoog12 and pyMOOGi13 , we measured a lithium abundance on the 12-point scale of A(Li) = 2.47 ± 0.04 dex.
To determine if this measured abundance strength is unusual, we compared this value to the lithium abundances of an ensemble of 392 control sample stars with reported lithium abundance measurements in the GALAH DR3 database (Buder et al. 2021).We selected control stars that exhibited stellar properties akin to that of our target.This included stars with a surface gravity of 3.9 ± 0.05 dex, an effective temperature of 6260 ± 50 K, an iron abundance ([Fe/H]) of −0.046 ± 0.05 dex, and a Gaia renormalized unit weight error (RUWE) value < 1.4, which indicates that the target is less likely to be part of a binary system (Belokurov et al. 2020).
The lithium abundances of our control population present a skewed lithium distribution.Therefore, we computed a modified z-score, which incorporates the median and absolute deviation from the median (MAD) for the control sample lithium abundances, as opposed to the population mean and standard deviation.Our target's modified Z-score is measured at 0.003, suggesting that its lithium abundance is nearly identical to the median of our dataset.Therefore, we report no indication of abnormal lithium contamination in the spectra of TOI-5301.

Building a Complete Sample of Hot Jupiters
This paper is the first in a new series of papers aimed at constructing a homogeneously analyzed set of HJ parameters, built off of the works of Rodriguez et al. (2019Rodriguez et al. ( , 2021)) 2021), roughly 300-400 transiting HJs around FGK stars with G < 12.5 mag should be discoverable with TESS , a sample that is now >∼ 60% complete.Once this sample is complete, which is achievable within the next few years with the help of TESS and TFOP, it will be possible to carefully consider the different possible migration mechanisms of HJs using detailed statistics with well-quantified biases.Models of other physical processes, such as the radius inflation of hot Jupiters subject to high instellation, will also be easier to test with this complete sample.
Individual discoveries of HJ systems with confounding features are constantly being made that will weigh on the models that we use to describe the evolution of systems with HJs (e.g., Becker et al. 2015;Wittenmyer et al. 2022;Yoshida et al. 2023).For example, since HJs that underwent high-eccentricity tidal migration are unlikely to have nearby companions due to the high likeli-hood of ejections and engulfment (Mustill et al. 2015), the exceptions to this rule will be key systems to test models of gas-disk migration.Additionally, photometric (e.g., De et al. 2023) and chemical (Montalbán & Rebolo 2002; see §4.4) signatures of planet engulfment may be found that provide insight into the stability of systems containing HJs.Throughout the construction of the sample, we will give greater consideration and care to benchmark systems in order to generate further constraints on models of migration.It is with the combination of a statistical sample of systems hosting HJs and the careful consideration of these benchmark systems that we hope to gain further insight into the origins and evolution of the enigmatic population of HJs.

SUMMARY
In this article, we present nine new hot Jupiters discovered by NASA's TESS mission as part of an ongoing effort to discover and characterize all HJs orbiting FGK stars brighter than G = 12.5 mag.This article is the first in a series of articles titled the Migration and Evolution of giant ExoPlanets (MEEP) survey.This effort aims to use homogenous detection and analysis techniques to generate a population of HJs with well-quantified selection biases.In pursuit of this effort, we are obtaining precise eccentricities, along with a host of other important parameters (such as planetary and stellar mass and radius), to test theories of giant planet formation and evolution.
The HJs we present in this work have masses ranging from 0.554 +0.076 −0.075 M J to 3.88 ± 0.23 M J and radii ranging from 0.967 +0.036 −0.031 R J to 1.438 +0.045 −0.042 R J , excluding the one planet in our sample that exhibits a grazing transit and for which the radius cannot be precisely constrained.Two of the planets in our sample, TOI-3919 b and TOI-5301 b, exhibit significant orbital eccentricities that imply that they may have undergone high-eccentricity tidal migration.We investigated the lithium absorption feature of the spectra of one star, TOI-5301, and found no evidence of an anomalous lithium enrichment that could be associated with youth or planetary engulfment.Each of these planets, however, will be an important member of the rising population of HJs that, once complete, will help answer the longest-standing question of exoplanet science, how do hot Jupiters form and evolve?(Collins et al. 2017), TAPIR (Jensen 2013), astropy (Astropy Collaboration et al. 2013, 2018), numpy (Harris et al. 2020), matplotlib (Hunter 2007), pandas (The pandas development Team 2023), scipy (Virtanen et al. 2020)        m/s/day, which was included in our fit per our prescription to fit for an RV trend if it is favored to greater than 1σ.

Figure 1 .
Figure 1.Adaptive optics images and sensitivity curves for TOI-3919 and TOI-5301.Top: Palomar PHARO AO image and 5σ sensitivity curve of TOI-3919 showing no detection of a companion star, down to the detection limits of the instrument.The magenta-shaded region represents the uncertainty on the contrast curve.Bottom: Keck2 NIRC2 AO image and 5σ sensitivity curve of TOI-5232, showing a clear detection of a stellar companion at a separation of 1.417 ′′ with ∆mag = 2.952 in the Kcont filter.

Figure 2 .
Figure2.Speckle high-resolution image and 5σ sensitivity curve from Gemini's Zorro speckle imager.The observation reveals a source 79.1 • E of N and 0.663 ′′ away from TOI-3321 that is 4.89 magnitudes fainter in Zorro's 832 nm narrowband filter.

8Figure 3 .
Figure 3. Posterior probability density function (PDF) of TOI-1855 b's radius.The black, labeled, vertical dashed line represents the posterior mode, 1.38 RJ.The cyan vertical dashed line and shaded region represent the median and 68% confidence interval, 1.65 +0.52−0.37RJ.The allowed planetary radius was constrained by an upper limit of 2.5 RJ.

Figure 4 .
Figure 4. Gaussian kernel density estimations of and TOI-5301's (right) stellar mass and age posterior distributions from our EXOFASTv2 global fits, illustrating a bimodal distribution in both parameters.The cyan dashed line indicates the median value of each parameter, while the shaded region represents the 68% confidence region.The normalized probability densities were estimated using the Gaussian kernel density estimator in the Python package scipy.stats.TOI-3894's fit slightly favored the lower mass solution with a probability of 53.0% while TOI-5301's fit favored the higher mass solution with a probability of 78.3%.

Figure 5 .
Figure 5. Eccentricity and semi-major axis distribution of the planets discovered in this work, compared to giant planets in the literature and several avenues and outcomes of giant planet migration (similar to Figure 4 from Dawson & Johnson 2018).The region labeled "Stellar Collision" corresponds to a giant planet colliding with the star, assuming the star is 1 R⊙ in size.The region labeled "Tidal Disruption" corresponds to a Jupiter-like planet falling within the Roche limit of the Sun.Finally, the red region corresponds to the tidal circularization of a highly eccentric giant planet around a Sun-like star.The blue circles represent the EXOFASTv2 median semi-major axis and eccentricity of the planets in this work, while the gray circles represent substellar bodies with precise eccentricities and masses between 0.25 -13 MJ, obtained from the NASA Exoplanet Archive.Many of the eccentric planets, including those that fall outside of the three shaded regions, can be explained by planet-planet scattering that has not yet or will not be tidally circularized.The equations that describe each of these equations are shown inDawson & Johnson (2018).

Figure 6 .
Figure 6.TESS , Follow-up and archival observations of TOI-1855 as they compare to the EXOFASTv2 results.Upper left: Unbinned TESS and follow-up ground-based transits, phase-folded and shown in comparison to the best-fit time of conjunction with an arbitrary normalized flux offset.Multiple TESS sectors in the same cadence are stacked on top of each other.Bottom left: The spectral energy distribution of the target star compared to the best-fit EXOFASTv2 model.Residuals are shown on a linear scale, using the same units as the primary y-axis.Upper right: RV observations versus time, including any significant long-term trend.The residuals are shown in the subpanel below in the same units.Middle right: RV observations phase-folded using the best-fit ephemeris from the EXOFASTv2 global fit.The phase is shifted so that the transit occurs at Phase + Offset = 0.The residuals are shown in the subpanel below in the same units.Bottom right: The evolutionary track and current evolutionary stage of the planet according to the best-fit MESA Isochrones and Stellar Tracks (MIST) model.The blue line indicates the best-fit MIST track, while the black contours show the 1σ and 2σ constraints on the star's current T eff and log g from the MIST isochrone alone.The green contours represent the 1σ and 2σ constraints on the star's T eff and log g from the EXOFASTv2 global fit, combining constraints from observations of the star and planet.The red cross indicates the median and 68% confidence interval reported in Table6.

PFigure 8 .
Figure8.Same as Figure6, but for TOI-2368.Notably, the best-fit RVs for TOI-2368 imply a linear RV trend of 0.21 ± 0.11 m/s/day, which was included in our fit per our prescription to fit for an RV trend if it is favored to greater than 1σ.

P
Figure10.Same as Figure6, but for TOI-3894.TOI-3894's EXOFASTv2 global fit was bimodal in stellar mass and age.To represent both modes, we show median values and standard deviations from each stellar mass solution in the bottom right plot.

P
Figure 14.Same as Figure 6, but for TOI-5301.Similarly to Figure 10, TOI-5301 had a bimodal distribution of possible stellar masses, and we therefore represent the median and standard deviation of either solution in the bottom right plot.

Table 2 .
Summary of Observations from TESS

Table 3 .
Follow-up observations

Table 4 .
The first RV measurement of each system, per instrument used.The full table of RVs for each system is available in machine-readable form in the online journal.

Table 5 .
Summary of High-Resolution Imaging Observations NOTE: All images and contrast curves are available on ExoFOP.

Table 6 .
Median Values and 68% Confidence Intervals for Fitted Stellar and Planetary Parameters

Table 7 .
Median Values and 68% Confidence Intervals for Bimodal Fits