Observations and statistical studies have shown that giant planets are rare around M dwarfs compared with Sunlike stars. The formation mechanism of these extreme systems has remained under debate for decades. With the help of the TESS mission and groundbased followup observations we report the discovery of TOI4201b the most massive and densest hot Jupiter around an M dwarf known so far with a radius of 1.22  0.04 R SUBJSUB and a mass of 2.48  0.09 M SUBJSUB about 5 times heavier than most other giant planets around M dwarfs. It also has the highest planettostar mass ratio q  4  10SUP3SUP among such systems. The host star is an early M dwarf with a mass of 0.61  0.02 M SUBSUB and a radius of 0.63  0.02 R SUBSUB. It has significant supersolar iron abundance FeH  0.52  0.08 dex. However interior structure modeling suggests that its planet TOI4201b is metalpoor which challenges the classical coreaccretion correlation of stellarplanet metallicity unless the planet is inflated by additional energy sources. Building on the detection of this planet we compare the stellar metallicity distribution of four planetary groups hotwarm Jupiters around GM dwarfs. We find that hotwarm Jupiters show a similar metallicity dependence around Gtype stars. For Mdwarf host stars the occurrence of hot Jupiters shows a much stronger correlation with iron abundance while warm Jupiters display a weaker preference indicating possible different formation histories.
['M dwarf stars', 'Radial velocity', 'Photometry', 'M stars', 'Extrasolar gaseous planets', '982', '1332', '1234', '985', '2172', 'Astrophysics - Earth and Planetary Astrophysics', 'Astrophysics - Solar and Stellar Astrophysics']
A Massive Hot Jupiter Orbiting a Metalrich Early M Star Discovered in the TESS Fullframe Images
{'A massive hot Jupiter orbiting a metal-rich early-M star discovered in the TESS full frame images': "Tianjun Gan , 1 Charles Cadieux , 2 Farbod Jahandar , 2 Allona Vazan , 3 Sharon X. Wang , 1 Shude Mao , 1, 4 Jaime A. Alvarado-Montes , 5, 6 D. N. C. Lin , 7, 8 ' Etienne Artigau , 2, 9 Neil J. Cook , 2 Ren'e Doyon , 2, 9 Andrew W. Mann , 10 Keivan G. Stassun , 11, 12 Adam J. Burgasser , 13 Benjamin V. Rackham , 14, 15, ∗ Steve B. Howell , 16 Karen A. Collins , 17 Khalid Barkaoui , 18, 14, 19 Avi Shporer , 20 Jerome de Leon , 21 Luc Arnold, 22 George R. Ricker, 20 Roland Vanderspek, 20 David W. Latham, 17 Sara Seager , 20, 14, 23 Joshua N. Winn, 24 Jon M. Jenkins, 16 Artem Burdanov , 14 David Charbonneau , 17 Georgina Dransfield, 25 Akihiko Fukui , 26, 27 Elise Furlan , 28 Michael Gillon , 18 Matthew J. Hooton , 29 Hannah M. Lewis , 30 Colin Littlefield , 31, 16 Ismael Mireles , 32 Norio Narita , 26, 33, 27 Chris W. Ormel , 1 Samuel N. Quinn , 17 Ramotholo Sefako , 34 18 20 14  \nMathilde Timmermans, Michael Vezie, and Julien de Wit  \n- 1 Department of Astronomy, Tsinghua University, Beijing 100084, People's Republic of China  \n2  \nUniversit'e de Montr'eal, D'epartement de Physique, IREX, Montr'eal, QC H3C 3J7, Canada  \nAstrophysics Research Center (ARCO), Dept. of Natural Sciences, Open University of Israel, Raanana, 4353701 Israel  \n3  \n- 4 National Astronomical Observatories, Chinese Academy of Sciences, 20A Datun Road, Chaoyang District, Beijing 100012, People's Republic of China  \n5 School of Mathematical and Physical Sciences, Macquarie University, Balaclava Road, North Ryde, NSW 2109, Australia 6 The Macquarie University Astrophysics and Space Technologies Research Centre, Macquarie University, Balaclava Road, North Ryde, NSW 2109, Australia  \n- 7 Department of Astronomy and Astrophysics University of California, Santa Cruz, CA 95064, USA  \n8 Institute for Advanced Studies Tsinghua University, Beijing 100086, People's Republic of China 9 Observatoire du Mont-M'egantic, Universit'e de Montr'eal, Montr'eal, QC H3C 3J7, Canada 10 Department of Physics and Astronomy, The University of North Carolina at Chapel Hill, Chapel Hill, NC 27599-3255, USA 11 Department of Physics and Astronomy, Vanderbilt University, 6301 Stevenson Center Ln., Nashville, TN 37235, USA 12 Department of Physics, Fisk University, 1000 17th Avenue North, Nashville, TN 37208, USA 13 Department of Astronomy & Astrophysics, UC San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA Department of Earth, Atmospheric and Planetary Science, Massachusetts Institute of Technology, 77 Massachusetts Avenue,  \n- 14 Cambridge, MA 02139, USA  \n15 Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA  \n16 NASA Ames Research Center, Moffett Field, CA 94035, USA  \n17 Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  \n18 Astrobiology Research Unit, Universit'e de Li'ege, 19C All'ee du 6 Aoˆut, 4000 Li'ege, Belgium  \n19 Instituto de Astrof'ısica de Canarias (IAC), Calle V'ıa L'actea s/n, 38200, La Laguna, Tenerife, Spain  \n20 Department of Physics and Kavli Institute for Astrophysics and Space Research, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  \n- 21 Department of Multi-Disciplinary Sciences, Graduate School of Arts and Sciences, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan  \n22 Canada-France-Hawaii Telescope, 65-1238 Mamalahoa Hwy, Kamuela, HI 96743, USA  \n- 23 Department of Aeronautics and Astronautics, MIT, 77 Massachusetts Avenue, Cambridge, MA 02139, USA\n- 24 Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA\n- 25 School of Physics & Astronomy, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK  \n26  \nKomaba Institute for Science, The University of Tokyo, 3-8-1 Komaba, Meguro, Tokyo 153-8902, Japan  \n27 Instituto de Astrof'ısica de Canarias (IAC), 38205 La Laguna, Tenerife, Spain  \n28 NASA Exoplanet Science Institute, Caltech/IPAC, Mail Code 100-22, 1200 E. California Blvd., Pasadena, CA 91125, USA  \n29 Cavendish Laboratory, JJ Thomson Avenue, Cambridge CB3 0HE, UK  \n30  \nSpace Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD, 21218, USA  \n31 Bay Area Environmental Research Institute, Moffett Field, CA 94035, USA  \nCorresponding author: Tianjun Gan  \ntianjungan@gmail.com  \n32 Department of Physics and Astronomy, University of New Mexico, 210 Yale Blvd NE, Albuquerque, NM 87106, USA 33 Astrobiology Center, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  \n34 South African Astronomical Observatory, P.O. Box 9, Observatory, Cape Town 7935, South Africa", 'ABSTRACT': 'Observations and statistical studies have shown that giant planets are rare around M dwarfs compared with Sun-like stars. The formation mechanism of these extreme systems remains under debate for decades. With the help of the TESS mission and ground based follow-up observations, we report the discovery of TOI-4201b, the most massive and densest hot Jupiter around an M dwarf known so far with a radius of 1 . 22 ± 0 . 04 R J and a mass of 2 . 48 ± 0 . 09 M J , about 5 times heavier than most other giant planets around M dwarfs. It also has the highest planet-to-star mass ratio ( q ∼ 4 × 10 -3 ) among such systems. The host star is an early-M dwarf with a mass of 0 . 61 ± 0 . 02 M ⊙ and a radius of 0 . 63 ± 0 . 02 R ⊙ . It has significant super-solar iron abundance ([Fe/H]=0 . 52 ± 0 . 08 dex). However, interior structure modeling suggests that its planet TOI-4201b is metal-poor, which challenges the classical core-accretion correlation of stellar-planet metallicity, unless the planet is inflated by additional energy sources. Building on the detection of this planet, we compare the stellar metallicity distribution of four planetary groups: hot/warm Jupiters around G/M dwarfs. We find that hot/warm Jupiters show a similar metallicity dependence around G-type stars. For M dwarf host stars, the occurrence of hot Jupiters shows a much stronger correlation with iron abundance, while warm Jupiters display a weaker preference, indicating possible different formation histories.  \nKeywords: planetary systems, planets and satellites, stars: individual (TIC 95057860, TOI-4201)', '1. INTRODUCTION': "As the most abundant stellar population in our Milky Way (Henry et al. 2006), M dwarfs are popular targets for exoplanet research. The Kepler mission (Borucki et al. 2010) has revealed that early-M stars host close-in small planets about 3 times more frequent than Sunlike stars (Petigura et al. 2013; Dressing & Charbonneau 2013, 2015; Gaidos et al. 2016). At the same time, however, recent work from Gan et al. (2023) found that hot Jupiters (defined as orbital period P ≤ 10 days, radius R p ≥ 7 R ⊕ ) are depleted around early-M dwarfs with an occurrence rate of 0 . 27 ± 0 . 09% (see also Bryant et al. 2023), in contrast to the frequency ∼ 0 . 6% around FGK stars (e.g., Fressin et al. 2013; Petigura et al. 2018; Zhou et al. 2019; Beleznay & Kunimoto 2022). This difference is even more significant if comparing with the result ∼ 1 . 0% around FGK stars from radial velocity (RV) surveys (e.g., Mayor et al. 2011; Cumming et al. 2008; Wright et al. 2012; Wittenmyer et al. 2020).  \nThe rarity of giant planets around M dwarfs is probably a natural outcome of the core accretion planet formation mechanism (Pollack et al. 1996). Under this hypothesis, the low mass of the protoplanetary disk around a low-mass star makes it difficult to form a massive enough solid core to start runaway gas accretion be-  \nfore disk dissipation (Laughlin et al. 2004; Ida & Lin 2005; Kennedy & Kenyon 2008). Indeed, several simulation works based on this paradigm reported a similar increasing trend of hot Jupiter occurrence rate as the stellar mass going up from 0.5 towards 1 M ⊙ (e.g., Liu et al. 2019; Burn et al. 2021). In addition, some initial studies found that most confirmed hot Jupiters are orbiting metal-rich M dwarfs with a median metallicity around 0.3 dex, which is higher than that of Sun-like stars harboring gas giants (Gan et al. 2022a; Kanodia et al. 2022; Kagetani et al. 2023). This is also supposed to be evidence in line with the core accretion scenario as the stellar metallicity is generally correlated with the mass of solids in the protoplanetary disk that are available for forming planets (Santos et al. 2004; Fischer & Valenti 2005; Sousa et al. 2011; Wang et al. 2018). In turn, the core accretion framework probably simultaneously explains the fact that small planets are common around M dwarfs because a forming outer giant planet around an FGK star may suppress the formation of close-in super-Earths due to the cut-off of the inward flow of solids (Mulders et al. 2021).  \nOn the other hand, there are still two relevant types of planetary systems that remain challenging for core accretion. The first population is gas giants around midto-late M dwarfs which theoretical works do not expect to form directly (e.g., Liu et al. 2019; Burn et al. 2021), unless a massive core was first formed through, for exam-  \nple, planet-planet collisions (Frelikh et al. 2019). Nevertheless, a few such systems were recently discovered (Morales et al. 2019; Hobson et al. 2023; Kanodia et al. 2023; Kagetani et al. 2023). Additionally, Bryant et al. (2023) also reported an occurrence rate of about 0.1% for hot Jupiters orbiting stars with M ∗ ≤ 0 . 4 M ⊙ , indicating that they are not vanishingly rare. Another group is cold Jupiters around M dwarfs. In particular, gas giants beyond the ice line are frequently discovered through the microlensing method (Mao & Paczynski 1991; Suzuki et al. 2016). By comparing the observations and predictions from planet population synthesis models, Schlecker et al. (2022) suggested that core accretion theories cannot reproduce the RV-detected giant planets around M dwarfs. Moreover, cold gas giants around M dwarfs detected by RV surveys show a much weaker metallicity dependence compared with hot Jupiters (Gan et al. 2022a), which hints that they follow a different formation path. In these cases, gravitational instability probably plays a role (Boss 2000).  \nThe TESS (Transiting Exoplanet Survey Satellite; Ricker et al. 2015) mission, which performs a full-sky photometric survey, has been rapidly increasing the number of gas giants around M dwarfs. Plenty of hot Jupiters have been detected over the last three years (e.g., HATS-71b, Bakos et al. 2020; HATS-74Ab and HATS-75b, Jord'an et al. 2022; TOI-3629b and TOI3714b, Ca˜nas et al. 2022). They have enlarged the sample size of short-period giant planets around M stars by a factor of three. A few warm Jupiters were also reported (e.g., Ca˜nas et al. 2020). All of these findings gradually enable glimpses of their formation and evolution history.  \nIn this manuscript, we report the discovery and confirmation of TOI-4201b, a massive and dense hot Jupiter around an early-M star, which was previously classified as a verified planet candidate by Gan et al. (2023). The rest of the paper is organized as follows. In Section 2, we detail all observations we use to confirm the planetary nature of TOI-4201b. We summarize the stellar properties in Section 3. Section 4 describes the joint-fit analysis to derive the planet parameters. In Section 5, we discuss the interior structure of TOI-4201b, tidal evolution analysis, the prospects for future characterizations as well as the comparison of stellar metallicity distribution of hot/warm Jupiters around G/M stars. We conclude with our findings in Section 6.", '2.1. TESS': "TOI-4201 (TIC 95057860) was first monitored by TESS in Sector 6 between 2018 December 15 and 2019 January 6 with the 30-minute cadence mode, and it was  \nrevisited every 10 minutes in Sector 33 from 2020 December 18 to 2021 January 13 during the first extended mission. The full frame images (FFIs) of TOI-4201 from Sector 6 were processed by the Quick Look Pipeline (QLP; Huang et al. 2020a,b; Fausnaugh et al. 2020). Due to the faintness of TOI-4201 ( T mag = 13 . 5), the light curve was not examined initially because the TESS vetting team only inspects targets with T mag ≤ 10 . 5. Its transiting signal was alerted by the QLP faint-star search program (Kunimoto et al. 2022), which searches for planet candidates around dimmer stars with 10 . 5 ≤ T mag ≤ 13 . 5 that have QLP light curves.  \nHowever, we note that the light curve of TOI-4201 from Sector 33 was not produced by QLP. Therefore, we carried out an independent and uniform aperture photometry for both Sectors using the lightkurve package (Lightkurve Collaboration et al. 2018; Barentsen et al. 2019). We first downloaded the 15 × 15 pixels FFI cutouts of TOI-4201 from each Sector (see the top panels of Figure 1), and extracted a raw target light curve using a 3 × 3 pixel custom aperture, centering at the location of TOI-4201. Next, we excluded a 5 × 5 pixel region around the target and constructed a 0 . 001 σ aperture mask to estimate the background flux. In this way, we picked up all pixels with flux smaller than 0.001 times the standard deviation of flux within the overall 15 × 15pixel region. Finally, we corrected the sky background variation by subtracting the background flux from the raw target flux. We show the light curves we obtained in Figure 1.  \nCombining all TESS data, we performed a transit search using the Transit Least Squares (TLS; Hippke & Heller 2019) algorithm after smoothing the light curve using a median filter with a window size of 0.3 days. We did not find additional significant periodic signals other than the 3.58-day one from TOI-4201b. After masking out all in-transit data, we fit a Gaussian Process (GP) model with a Mat'ern-3/2 kernel using celerite (Foreman-Mackey et al. 2017) for the light curve from each Sector. We divided the light curve by the best-fit GP model for detrending and normalization.", '2.2. Ground-based photometry': 'In order to rule out false positive scenarios such as nearby and blended eclipsing binary, and refine the transit ephemeris as well as transit depth measurement, we collected a total of eight ground-based follow-up light curves for TOI-4201, as part of the TESS Follow-up Observing Program (TFOP). We scheduled all observations based on the transit information from the TESS Transit Finder (TTF), which is a customized version of the Tapir software package (Jensen 2013). We summarize  \nFigure 1. Top panels: The 15x15 pixel TESS full frame image from Sectors 6 and 33. The black-shaded region represents the 3x3 pixel custom aperture we used to conduct photometry. TOI-4201 is marked as a white star at the center. Different sizes of red circles represent different magnitudes in contrast with TOI-4201. Bottom panels: The normalized TESS light curves of TOI-4201 from Sectors 6 (left) and 33 (right) that we extracted. The red solid lines are the GP models used for detrending. The detrended light curves are shown below. The transits of TOI-4201b are marked in blue ticks.  \n<!-- image -->  \n(c) TESS light curve from Sector 6  \n(d) TESS light curve from Sector 33  \nthe details of all observations in Table 1, and describe each of them below.', '2.2.1. LCOGT': "We acquired two SDSS g ' and i ' band alternating time-series observations using 1-meter telescopes from the Las Cumbres Observatory Global Telescope network (LCOGT; Brown et al. 2013). The first partial transit (ingress-only) observation was done at South African Astronomical Observatory (SAAO) on 2021 September 26, while the second one was obtained at Cerro Tololo Interamerican Observatory (CTIO) on 2021 October 4 covering the full event. Both observations were carried out with the Sinistro cameras, which have a field of view of 26 ' × 26 ' and a plate scale of 0 . 389 '' /pixel. The exposure times of g ' and i ' band observations are 300s and 180s, respectively. After the raw images were  \ncalibrated by the automatic BANZAI pipeline (McCully et al. 2018), we conducted a photometric analysis using the AstroImageJ software (Collins et al. 2017) with a 22-pixel (8 . 5 '' ) and 16-pixel (6 . 2 '' ) aperture for two observations. We confirmed the transit signal on target at the ephemeris provided by TESS.", '2.2.2. MuSCAT': "We observed a full transit of TOI-4201b on 2022 Jan 30 with MuSCAT using exposure times of 120, 50, 50 s for g , r , and z s bands, respectively. MuSCAT is a multi-band simultaneous camera installed on the 188-cm 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 ' FOV, enabling simultaneous photometry in the g (400-550 nm), r (550-700 nm), and z s (820-920 nm) bands. The data  \nreduction 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. Five comparison stars and aperture radii of 12 pixels (4.2 '' ) yield the optimum light curves for all bands.", '2.2.3. SPECULOOS-North': "We obtained a full transit of TOI-4201 b by the SPECULOOS-North telescope on 2023 February 4 in the Sloanz ' filter with an exposure time of 20s. SPECULOOS-North is a 1.0-m Ritchey-Chretien telescope equipped with a thermoelectrically cooled 2K × 2K Andor iKon-L BEX2-DD CCD camera with a pixel scale of 0 . 35 '' and a field-of-view of 12 ' × 12 ' (Burdanov et al. 2022). It is a twin of the SPECULOOS-South located in ESO-Paranal in Chile (Jehin et al. 2018; Delrez et al. 2018; Sebastian et al. 2021) and SAINT-EX located at the Sierra de San Pedro M'artir in Baja California, M'exico (Demory et al. 2020). Data reduction and photometric measurements were performed using the PROSE 1 pipeline (Garcia et al. 2021). The resulting light curve is presented in Figure 2.", '2.3. Spectroscopic follow-up': "2.3.1. MagE  \nWe obtained medium-resolution optical spectra of TOI-4201 on 2022 October 06 (UT) with the Magellan Echellette (MagE) spectrograph (Marshall et al. 2008) on the 6.5-m Magellan Baade Telescope. Conditions were clear with seeing of 0 . '' 7, and we observed the target during nautical twilight. We used the 0 . '' 7 × 10 '' slit to obtain resolution λ/ ∆ λ ≈ 6000 over 3200-10 000 ˚ A. We collected two, 290-s exposures at an airmass of 1.04. We also observed the spectrophotometric calibrator Feige 110 during the night for flux calibration (Hamuy et al. 1992, 1994), and obtained bias exposures as well as ThAr arc lamp and Xe flash and incandescent flatfield lamp exposures at the start of the night for wavelength and flux calibration, respectively. We did not observe a telluric absorption calibrator for these observations, hence telluric features remain in the spectra. Data were reduced using PypeIt (Prochaska et al. 2020; Prochaska et al. 2020) using standard settings. The final calibrated spectrum has a median signal-to-noise ≈ 100 in the 8000 ˚ A region.  \nFigure 3 compares the MagE spectrum of TOI-4201 to the best-match spectral template constructed from  \nFigure 2. The TESS and ground-based light curves folded in phase with the orbital period of TOI-4201b. The red solid lines are the best-fit transit models. Instrument name and filter information are given at the bottom left of each panel.  \n<!-- image -->  \nTable 1. Ground-based photometric follow-up observations for TOI-4201  \nFigure 3. MagE optical spectrum of TOI-4201 in normalized F λ units (black line, uncertainties in grey), compared to the best-match M1 dwarf spectral template (magenta line; data from Bochanski et al. 2007). We also label typical spectral absorption features for M dwarfs and regions of uncorrected telluric absorption ( ⊕ ). The inset box shows a closeup of the 6530-6780 ˚ A region indicating H α in absorption and no Li I detection.  \n<!-- image -->  \nSloan Digital Sky Survey data (Bochanski et al. 2007). Note that the template has a lower mean resolution ( λ/ ∆ λ ≈ 2000), but is nonetheless an excellent match, indicating a dwarf classification. This template match is confirmed by index-based classifications from Reid et al. (1995); Gizis (1997); and L'epine et al. (2003), which span M0.5-M1.5. We therefore adopt a spectral classification of M1.0 ± 0.5 for this source. Inspection of the data show H α in absorption with an equivalent width of 0.426 ± 0.017 ˚ A, and no detection of Li I, ruling out a young ( ≲ 10 Myr) source. Using the calibration of L'epine et al. (2013), we measure a ζ metallicity index of 1.166 ± 0.008, indicating a slightly supersolar metallicity. The metallicity relation of Mann et al. (2013) yields [Fe/H] = 0.22 ± 0.20 for this ζ value, where the uncertainty is dominated by systematics among the calibration sample. We conduct a more detailed abundance analysis using the SPIRou data in the following.", '2.3.2. SPIRou': "We obtained a total of 28 spectroscopic observations for TOI-4201 using SPIRou (SpectroPolarim'etre InfraROUge; Donati et al. 2020) between 2022 November 6 and 2023 January 1. SPIRou is a near-infrared  \nhigh resolution ( R ≈ 75 , 000) fiber-fed echelle spectrograph installed on the 3.6m Canada-France-Hawaii Telescope (CFHT), covering a wavelength range from 0.98 to 2.5 µ m. Since TOI-4201 is a faint target ( H mag > 11), instead of taking simultaneous drift calibration with the thermalized Farby-P'erot (FP) etalon, we opted to use the Dark mode to avoid contaminating science spectra. We conducted all of our observations under an environment of airmass around 1.3 and seeing about 0 . 6 '' . Every night, we took two continuous sequences of TOI-4201 at a random phase, with a 1200s exposure time for each one. For our 28 spectra, the median signal-to-noise ratio (SNR) per pixel in the middle of H band is about 40.  \nThe SPIRou data were reduced using APERO version 0.7.275 (Cook et al. 2022). The main steps of data reduction are described in the following but readers are referred to Cook et al. (2022) for a complete description of the APERO pipeline and modules. Raw frames (science and calibration) are first pre-processed to remove known structures on the 4096 × 4096 detector. Nightly calibration images are taken to locate the echelle orders and to perform a series of corrections (e.g., flat, thermal background, blaze, etc.). The flux is extracted from science and calibration frames to produce per-order spectra. A nightly pixel-to-wavelength solution is derived following the method of Hobson et al. (2021). Finally, a three steps telluric correction is applied to the science data. A first cleaning is done with a TAPAS (Bertaux et al. 2014) atmospheric transmission model. The remaining residuals (percent-level in deep H 2 O bands) are then fitted in this second step using a grid of telluric models generated from observations of hot stars (fast rotators) with SPIRou covering ranges of airmass and water column. The last step mitigates finite resolution effects, i.e., the fact that deconvolution between stellar+telluric spectrum with instrumental profile is always imperfect (Wang et al. 2022). This telluric correction method yields final residuals at the level of the PCA-based approach of Artigau et al. (2014) and will be described in more details in Artigau et al. in prep.  \nThe radial velocity (RV) measurements were obtained from the telluric corrected spectra using the line-byline (LBL) method of Artigau et al. (2022). In this  \nframework, Doppler shifts are measured on small chunks of the spectrum called 'lines' ( ∼ 16,000 lines for an M dwarfs observed with SPIRou) using the Bouchy et al. (2001) formalism. With thousands of independent velocity measurements over the full domain, outlying spectral features are readily identified and removed. This is particularly important in the near-infrared where spurious features introduced by complex detector architecture and imperfect telluric correction can significantly bias RV observations (Artigau et al. 2022). The final RV is obtained from the error-weighted average of all valid per-line velocities. The LBL method enabled the mass detection of several TESS planets with SPIRou in recent years (e.g., TOI-1759 b, Martioli et al. 2022; TOI2136 b, Gan et al. 2022b; TOI-1452 b, Cadieux et al. 2022; TOI-1695 b, Kiefer et al. 2023). After subtracting the systemic velocity of about 42.1 km s -1 , we list the LBL RVs of TOI-4201 in Table 2, which have a median uncertainty of 24 m s -1 . As seen in Figure 5, the Doppler signal of the transiting companion is clearly detected.  \nWe then performed a generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kurster 2009) analysis on the SPIRou RVs (see Figure 4). A significant peak around 3.58 days is detected at the orbital period of the transiting planet TOI-4201b. In addition, another two signals at 0.78 and 1.38 days can be clearly seen in the periodogram. These two signals come from the alias of the orbital period with the window function, which correspond to the frequency of f = 1 + 1 / 3 . 58 day -1 and f = 1 -1 / 3 . 58 day -1 . After subtracting the best Keplerian model from the joint-fit (see Section 4), we find that they disappears and no other significant peaks, which may point to the existence of a secondary planet, show up with FAP below 0.1%.", '2.4. High Angular 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 background 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, 2020). Additionally, the presence of a close companion star leads to the non-detection of small planets residing with the same exoplanetary system (Lester et al. 2021). Given that nearly one-half of solar-like stars are in binary or multiple star systems, (Matson et al. 2018) high-resolution imaging provides crucial information toward our understanding of exoplanet formation, dynamics and evolution (Howell et al. 2021).  \nFigure 4. Top panel: The GLS periodogram of SPIRou radial velocities of TOI-4201. Bottom panel: The GLS periodogram of SPIRou RVs but after subtracting the bestfitting Keplerian model of TOI-4201b. In both panels, the theoretical FAP levels of 10, 1, and 0.1 percent are marked as horizontal dashed, dot-dashed and dotted lines. The orbital period of TOI-4201b is marked as a vertical solid line. Another two significant peaks around 0.78 and 1.38 days (vertical dashed lines) are the signals due to frequency beating.  \n<!-- image -->  \nTOI-4201 was observed on 2023 April 25 UT using the Zorro speckle instrument on the Gemini South 8-m telescope (Scott et al. 2021; Howell & Furlan 2022). Zorro provides simultaneous speckle imaging in two bands (562nm and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections. Nine sets of 1000 × 0.06 second images were obtained and processed in our standard speckle imaging reduction pipeline (see Howell et al. 2011). Figure 6 shows our final 5 σ contrast curves and the 832 nm reconstructed speckle image. The observations of TOI-4201 revealed it to have no close companions within the angular and contrast levels achieved. We find no companion brighter than 4-5 magnitudes below that of the star itself from the 8-m telescope diffraction limit (20 mas) out to 1.2 '' . At the distance of TOI-4201 (d=189 pc), these angular limits correspond to spatial limits of 3.8 to 227 AU.", '3.1. Stellar parameters': 'We apply three different methods to determine the stellar properties of TOI-4201 including SED fitting, spectroscopic analysis and empirical relations. We take the weighted stellar radius, mass, effective temperature, and list them in Table 3.  \n3.1.1. SED fitting  \nFigure 5. Left panel: The time-series SPIRou radial velocities of TOI-4201. Right panel: The RVs phased-folded to the orbit period of TOI-4201b. In both panels, the black solid lines are the best Keplerian model from joint-fit. The orange shaded region denotes the 1 σ credible intervals of the RV model. Residuals are plotted below.  \n<!-- image -->  \nFigure 6. Gemini high-resolution speckle imaging 5 σ contrast limits as a function of angular separation. We show a separate curve for each filter and the 832 nm reconstructed image. TOI-4201 has no close companions within the angular and contrast levels achieved.  \n<!-- image -->  \nWe carry out a spectral energy distribution (SED) analysis of the star together with the Gaia DR3 parallax (with no systematic offset applied; see, e.g., Stassun & Torres 2021), in order to determine an empirical measurement of the stellar radius, following the procedures described in Stassun & Torres (2016); Stassun et al.  \n(2017); Stassun & Torres (2018). We pulled the JHK S magnitudes from 2MASS (Cutri et al. 2003), the W1W3 magnitudes from WISE (Wright et al. 2010), the zy magnitudes from Pan-STARRS (Chambers et al. 2016), and the G , G BP , G RP magnitudes from Gaia DR3 (Gaia Collaboration et al. 2023). Together, the available photometry spans the full stellar SED over the wavelength range 0.4-10 µ m (see Figure 7).  \nWe perform a fit using PHOENIX stellar atmosphere models (Husser et al. 2013), with the free parameters being the effective temperature ( T eff ) and metallicity ([Fe/H]), as well as the extinction A V , which we limited to maximum line-of-sight value from the Galactic dust maps of Schlegel et al. (1998). The resulting fit (Figure 7) has a best-fit A V = 0 . 03 ± 0 . 02, T eff = 3800 ± 75 K, [Fe/H] = 0 . 2 ± 0 . 1, with a reduced χ 2 of 1.4. The metallicity in particular is further constrained by the Gaia spectrum (see grey swathe in Figure 7 and inset). Integrating the (unreddened) model SED gives the bolometric flux at Earth, F bol = 6 . 80 ± 0 . 16 × 10 -11 erg s -1 cm -2 . Taking the F bol and T eff together with the Gaia parallax, gives the stellar radius, R ∗ = 0 . 63 ± 0 . 03 R ⊙ .  \nWe also independently derive the stellar effective temperature ( T eff ) and radius R ∗ following the method outlined in (Mann et al. 2016). To briefly summarize, we fit the SED by comparing the available photometry to a grid of late-K and early-M dwarf templates from (Gaidos et al. 2014). To fill in gaps in the templates, we  \nTable 2. SPIRou radial velocities for TOI-4201  \nuse PHOENIX BT-SETTL models (Allard et al. 2013). The free parameters are the template choice, the scale factor, reddening, model selection ( T eff and [M/H]), and three parameters meant to account for systematics if the flux calibration of the templates. We locked the reddening to E ( B -V ) < 0 . 1 based on 3D extinction maps (Green et al. 2019), although this did not change the final fit. We estimate the luminosity from the integral of the absolutely-calibrated spectrum and Gaia DR3 parallax, then estimate R ∗ from the Stefan-Boltzmann relation as well as the infrared-flux method (Blackwell & Shallis 1977). The final resulting fit is T eff = 3784 ± 83 K and R ∗ = 0 . 65 ± 0 . 03 R ⊙ , consistent with other estimates.  \nThe planet host sits high on the color-magnitude diagram for a main-sequence M dwarf (Mann et al. 2015), suggesting it is metal-rich ([Fe/H] > +0 . 2 dex). However, deriving a metallicity from this method assumes a single star with no reddening. Instead, we estimate  \nFigure 7. Spectral energy distribution of TOI-4201. Red symbols represent the observed photometric measurements, where the horizontal bars represent the effective width of the passband. Blue symbols are the model fluxes from the bestfit PHOENIX atmosphere model (black). The Gaia spectrum is shown as a grey swathe; the inset focuses on the Gaia spectrum.  \n<!-- image -->  \n[Fe/H] using the JHK S photometry as outlined in Newton et al. (2014). This yielded [Fe/H]= 0 . 3 ± 0 . 1, although it may be an underestimate as the host star sits above most of the calibration stars (see Figure 21 in Newton et al. 2014).', '3.1.2. Spectroscopic analysis': 'The SPIRou high-resolution combined spectrum of TOI-4201 offers an independent determination of its chemical composition. Such characterization on SPIRou data was also performed on the M4 dwarf TOI-1452 in Cadieux et al. (2022) and additional details on the methodology presented below are given in Jahandar et al. in prep. Briefly, the combined spectrum of TOI-4201 is compared to a 2D grid ( T eff and [M/H]) of ACES stellar models (Allard et al. 2012a; Husser et al. 2013) deconvolved to match SPIRou resolution. We select strong absorption lines matching the models to avoid continuum mismatch and with a known origin, i.e., atomic or molecular and found in the PHOENIX/BT-Settl (Allard et al. 2012b; 2013) and NIST (Kramida et al. 2022) line lists. Fixing the T eff to the value derived above, we performed fits on individual spectral lines to determine the abundances of several chemical species by varying the overall metallicity [M/H]. Since the combined spectrum of TOI-4201 has a relatively low SNR of ∼ 200 in H band, we also ran Monte Carlo simulations (varying flux values) to assess the robustness of the fits. We find a high stellar iron abundance of 0 . 52 ± 0 . 08 dex for TOI-4201 with 12 Fe I absorption lines, which is consistent with  \nFigure 8. SPIRou observation of OH and Fe I lines of TOI4201 (black dashed line). The solid lines represent the ACES models for different overall metallicities [M/H] and the fixed T eff of 3800 K. These plots illustrate the good sensitivity of near-IR high-resolution spectroscopy for constraining the chemical abundances of OH and Fe lines.  \n<!-- image -->  \n1911.70 1911.75 1911.80 1911.85 1911.90 1911.95 1912.00 1912.05 1912.10 Wavelength [nm]  \n<!-- image -->  \nthe estimations above from the photometry within 2 σ . The final [Fe/H] from the spectroscopic analysis is based on the average of best-fit overall metallicity [M/H] of each Fe line. The uncertainty of [Fe/H] provided in the work is the standard variation of these [M/H] divided by √ N -1, where N is the total number of lines (12). We note that the error bar might be underestimated since this methodology has an intrinsic RMS of about 0.2 dex (see Figure 11 in Cristofari et al. 2022, which used the same synthetic model and observed data, although different choice of lines.). We use the [Fe/H] from SPIRou spectrum analysis for further investigation in the rest of this work and list all metallicity measurements from different sources in Table 3. Meanwhile, we see very shallow OH lines, which indicates a low oxygen abundance (see Figure 8). High-resolution spectroscopic observations in optical bands are required to confirm these abundance measurements.', '3.1.3. Empirical relation': 'Based on the m K from 2MASS and parallax from Gaia DR3, we obtain the absolute magnitude M K of TOI4201 to be 4 . 98 ± 0 . 03 mag. We then estimate the stellar radius by employing the empirical relation between R ∗ and M K derived by Mann et al. (2015), and we find an R ∗ of 0 . 61 ± 0 . 02 R ⊙ . This is consistent with the value 0 . 65 ± 0 . 03 R ⊙ within 1 σ measured using the relation between stellar angular diameter and color V -K s described in Boyajian et al. (2014).  \nNext, we estimate the stellar effective temperature T eff using two methods. First, we calculate the bolometric correction BC K using its relation with V -J from Mann et al. (2015), and we find BC K = 2 . 55 ± 0 . 04 mag. Therefore, we measure a bolometric magnitude of M bol = 7 . 54 ± 0 . 05 mag, which results in a bolometric luminosity of L ∗ = 0 . 076 ± 0 . 004 L ⊙ . The effective temperature 3800 ± 77 K is then derived using the Stefan-Boltzmann law. Furthermore, we also obtain an empirical estimation of T eff = 3803 ± 82 K using a polynomial relation with stellar colors V -J and J -H (Mann et al. 2015).  \nFinally, we obtain the stellar mass M ∗ using Equation 2 in Mann et al. (2019) according to the M ∗ -M K relation 2 . We find M ∗ = 0 . 60 ± 0 . 01 M ⊙ , which agrees with 0 . 63 ± 0 . 02 M ⊙ based on the result from Benedict et al. (2016). We in turn use the mass-radius relation from Boyajian et al. (2012) (see their Eq. 10), and we find a stellar radius of 0 . 61 ± 0 . 02 R ⊙ , consistent with other results above.', '3.2. Stellar rotation and age estimation': 'The TESS light curves of TOI-4201 that we extracted using a custom aperture from both Sectors 6 and 33 have an obvious baseline variation, which implies a high stellar rotation speed. A similar modulation is also shown in the QLP light curve of Sector 6. We conduct a frequency analysis using the generalized Lomb-Scargle (GLS) periodogram (Zechmeister & Kurster 2009) after removing all in-transit data to measure the stellar rotation period. The result is presented in Figure 9. The light curve from Sector 6 shows a signal around 16.4 days as well as a significant aliasing signal peaking at 8.2 days. While the one from Sector 33 has a broad signal around 20.4 days, similar aliasing signals could also be seen around 8 days. Both the 16.4-day and 20.4-day signals have a false alarm probability (FAP) below 0.1%.  \nDue to the limited length of the TESS observation baseline, the constraint on the rotation period is poor using TESS-only data. Therefore, we further investigate the long-term ground photometry from Zwicky Transient Facility (ZTF; Bellm et al. 2019; Masci et al. 2019), and search for periodic modulations. ZTF is mounted on the 48-inch aperture Schmidt-type Telescope at the Palomar Observatory. Using a wide-field camera consisting of 16 CCDs, which results in a field of view of 47 squared deg, ZTF scans the entire Northern sky every two days, providing a large amount of data for timedomain science. TOI-4201 is observed by ZTF using  \nTable 3. Summary of stellar parameters for TOI-4201  \n[1]Stassun et al. (2019); [2]Gaia Collaboration et al. (2023), [3]Henden et al. (2016); [4]Cutri et al. (2003); [5]Wright et al. (2010); [6] From photometric analysis; [7] From spectroscopic analysis of MagE data; [8] From spectroscopic analysis of SPIRou data.  \ntwo different CCDs (id = 2 and 16) in both g and r bands. All of these observations were done with an exposure time of 30s and data are publicly available under DR16. After removing all flux measurements flagged as bad quality, we find that ZTF has (i) 170 measurements from CCD-2 in g band between 2018 Marth 27 and 2022 December 27; (ii) 184 measurements from CCD-16 in g band between 2018 Marth 27 and 2022 December 27; (iii) 206 measurements from CCD-2 in r band between 2018 September 4 and 2022 December 27; and (iv) 183 measurements from CCD-16 in r band between 2018 September 6 and 2022 December 27. We compute the GLS periodogram for each of these datasets, and we  \nFigure 9. The generalized Lomb-Scargle periodograms of the TESS and ZTF photometry. The black solid and dashed lines are the results of TESS while the colored lines are from ZTF, with filter and CCD id shown in the legend. The vertical gray line represents the ∼ 17 d rotational signal of TOI4201. A zoomed plot is shown on the top right.  \n<!-- image -->  \nshow the results in Figure 9. Except for the first one that shows no significant signals with FAP below 10%, all the other three light curves have a peak around 17.34 days with a FAP of ∼ 0.1%, ∼ 10% and ∼ 1%, which is consistent with the aforementioned findings from TESS data. A forest of other peaks can also be seen between 15 and 20 days, probably due to the poor sampling.  \nBased on Kepler observations, McQuillan et al. (2014) analyzed the rotation periods of a series of main sequence stars with effective temperatures below 6500 K. They suggested that a typical M star like TOI-4201 with T eff around 3800 K has a rotation period between 10 and 40 days. In addition, Newton et al. (2016, 2018) made use of ground long-term photometric monitoring from the MEarth survey (Irwin et al. 2009), which leads to a conclusion that the rotation periods of early-M stars are generally within 40 days, although their sample has a limited number of early-M stars. Therefore, we attribute this 17.34-day signal to the stellar rotation.  \nBuilding on the stellar rotation period, we next evaluate the stellar age. We first adopt the empirical relation from Engle & Guinan (2018), which yields 1 . 6 ± 0 . 9 Gyr, in agreement with the estimation 785 ± 85 Myr from (Mamajek & Hillenbrand 2008). Additionally, we use the gyrochronology relation reported by Barnes (2007). Without taking reddening into consideration, we also find an age below 2 Gyr given the stellar color B -V from APASS (Henden et al. 2016). We thus provide a conservative age estimation of TOI-4201 between 0.7 and 2.0 Gyr.', '3.3. Stellar kinematic properties': "Following the method described in Johnson & Soderblom (1987), we calculate the 3-dimensional ve-  \nlocity of TOI-4201 with respect to the local standard of rest (LSR). We utilize the stellar proper motion ( µ α and µ δ ) and parallax ( ϖ ) from Gaia DR3 as well as the systemic velocity RV from SPIRou measurement. Taking the solar velocity components relative to the LSR ( U ⊙ , V ⊙ , W ⊙ )=(9.58, 10.52, 7.01) km s -1 from the Large Sky Area Multi-Object Fiber Spectroscopic Telescope (LAMOST; Tian et al. 2015) into consideration, we obtain a Galactic space motion of ( U LSR , V LSR , W LSR )=( -23 . 9 ± 0 . 2, -17 . 6 ± 0 . 2, 6 . 0 ± 0 . 1) km s -1 . We further make use of the BANYAN Σ algorithm (Gagn'e et al. 2018) to determine the membership probability of TOI-4201 within young associations. Using the stellar kinematic parameters above, we find that TOI-4201 is likely a field star and it shows no evidence of cluster membership. We finally follow the procedure in Bensby et al. (2003, 2014), and estimate a low probability ratio about P thick /P thin =1.3% for TOI-4201 belonging to the Galactic thick and thin disk, indicating that TOI-4201 has a thin-disk origin.", '4. JOINT-FIT': 'We carry out a joint-fit of all photometry and SPIRou radial velocities to derive the planetary physical parameters using the juliet package (Espinoza et al. 2019). Basically, we fit light curves using batman (Kreidberg 2015) while the RV model is generated with radvel (Fulton et al. 2018). We obtain the posteriors of all parameters through a dynamic nested sampling with dynesty (Speagle 2020).  \nWe set Gaussian normal priors on both the orbital period and mid-transit epoch, centering at the values we found in our transit search, with a 1 σ value of 0.1 days. Following the parametrization for the planet-to-star radius ratio p and impact parameter b in Espinoza (2018), we efficiently sample their physical parameter spaces by fitting r 1 and r 2 instead and allow them uniformly vary between 0 and 1. Since we use a relatively large box aperture to extract the TESS photometry, light contamination from other stars especially the nearby star Gaia 2997312063605005952 located within the aperture with G = 16 . 8 mag (∆ G = 2 . 3 mag) should be considered (see Figure 1). Therefore, we fit a dilution factor 3 for the TESS photometry while we fix it at 1 for ground-based observations as nearby stars around TOI-4201 in ground images are deblended. We opt to adopt a quadratic limb darkening law for the TESS photometry and a linear law for ground data (Kipping 2013). Finally, we include a  \njitter term for every photometric dataset to account for additional white noise.  \nIn terms of the RV modeling, we fit a standard Keplerian orbit, leaving e sin ω and e cos ω as free parameters. We place wide informative uniform priors on the RV semi-amplitude K , the systematic velocity µ SPIRou as well as the RV jitter σ SPIRou . Due to the short time span, we do not consider linear and quadratic RV trends here, both of which are fixed at 0. The joint fit reveals that the planet has a radius of 1 . 22 ± 0 . 04 R J with a mass of 2 . 48 ± 0 . 09 M J on a nearly circular orbit. We list the prior settings and the posteriors of all parameters in Table 4. We show the phase-folded data along with the best-fit transit and RV models in Figures 2 and 5.', '5.1. TOI-4201b: A massive and dense giant planet': "With a mass of around 2 . 5 M J , TOI-4201b is the most massive hot Jupiter around an M dwarf known so far (see the top panel of Figure 10). It stands out from other similar systems (M dwarf + hot Jupiter) with a much lower companion mass around 0 . 5 M J . Since most host stars of those systems are also early-M dwarfs, the reason why TOI-4201b is more massive than others is unclear. In addition, TOI-4201b also has the highest planet-to-star mass ratio of about 4 × 10 -3 among M dwarf hot Jupiter systems (see the bottom panel of Figure 10). Compared with warm and cold Jupiters detected by RV and microlensing surveys, it turns out that such massive close-in Jupiters are lacking around M dwarfs. Together with HATS-74Ab (Jord'an et al. 2022) and TOI-5205b (Kanodia et al. 2023), TOI-4201b is the third planet situated in the mass ratio paucity region ( q ≥ 2 × 10 -3 and a ≤ 0 . 1 AU), an unpopulated area proposed by Gan et al. (2022a), making it an interesting object for theories of planet formation.  \nFigure 11 shows the Mass-Radius diagram of all confirmed transiting gas giants around FGKM dwarfs. Unlike those giant planets around FGK stars that have both diverse radius and mass, gas giants around M dwarfs have masses spanning orders of magnitude but with relatively concentrated radius, which is mainly due to the low incident flux they received. Because these giant planets around M dwarfs are not strongly irradiated by the host star, they are good targets to study interior structures.", '5.2. Interior structure of TOI-4201b': 'The mean density of TOI-4201b is about 1 . 82 ± 0 . 19 g cm -3 , significantly higher than other observed gas giants around M-dwarfs, as is shown in Figure 11. The  \nTable 4. Parameter priors and the best-fit values along with the 68% credibility intervals in the final joint fit for TOI-4201. N ( µ , σ 2 ) means a normal prior with mean µ and standard deviation σ . U (a , b) stands for a uniform prior between a and b . LU (a , b) is a log-uniform prior between a and b . T N ( µ , σ 2 , a, b) represents a truncated normal prior ranging from a to b .  \n- [1] We do not consider heat distribution between the dayside and nightside here and assume albedo A B = 0.  \n<!-- image -->  \nFigure 10. Top panel: The mass and semi-major axis diagram of all confirmed gas giants around M dwarfs from transit (red) and RV (blue) surveys. Bottom panel: The planet-to-star mass ratio and semi-major axis distribution of gas giants around M dwarfs colored by the detection techniques. The blue shaded region is the mass ratio paucity region proposed by Gan et al. (2022a). In both panels, TOI4201b is marked as a brown star.  \n<!-- image -->  \nhigh density of the planet can be entirely attributed to hydrogen compression in the massive planet, and not of high metal content, as in the case of lower mass planets, since high Z would make the planet radius too small. Interior-evolution modeling suggests metal content between 0-25 M ⊕ in the interior of TOI-4201b, depending on the metal distribution in the interior, its history, and its energy state.  \nAccording to the standard core accretion models (e.g., Fortney et al. 2007; Baraffe et al. 2008), the mean density of TOI-4201b is consistent with pure H, He composition and lighter than that, i.e. inflated radius. We perform dedicated structure-evolution models of TOI4201b under various conditions, based on Vazan et al.  \nFigure 11. Planet mass and radius diagram. The gray dots are giant planets around FGK stars. The colored dots are all transiting gas giants around M dwarfs with mass measurements. Different colors represent different equilibrium temperatures. Three constant density curves ( ρ =0.5, 1.0, 2.0 g cm -3 ) are shown as dashed lines for reference. With a density around 1.8 g cm -3 , TOI-4201b is one of the densest gas giants transiting M dwarfs.  \n<!-- image -->  \n(2013, 2015) in order to fit its radius-mass-age parameters. Since atmospheric properties of TOI-4201b are not available, we consider albedo values of 0-0.3, and atmospheric opacity metallicity between solar and ten times solar (Freedman et al. 2014). Evolution models covered a wide range of parameters of cold and hot start conditions (e.g., Marley et al. 2007), and various metal distributions in the interior. The models suggest a core mass of 0-5 Earth masses and a metal-poor envelope. Alternative scenarios that involve post-formation giant impact allow larger metal content. In such cases TOI-4201b can contain up to 25 Earth masses of metals, gradually distributed in the deep interior, surrounded by a metalpoor envelope.  \nIn Figure 12, we show the radius evolution of three TOI-4201b mass planets. Under standard evolution conditions TOI-4201b is an inflated hot Jupiter, since its observed radius-age (grey box) is larger than any model with efficient (adiabatic) cooling, even if the model is of a metal-free planet (dashed blue). The hot diluted core model (red) simulates post energetic giant impact scenario that eroded parts of the core, injecting energy and slowing the cooling of the planet. Yet, the maximum metal mass in this scenario cannot exceed 25 M ⊕ , and the model meets the observed values only at a limited range. A classical core-envelope model with 5 M ⊕ core and high atmospheric opacity of ten times solar (green) barely fits the observed properties. These models indi-  \nFigure 12. Radius evolution of three TOI-4201b like planets. The observed values of radius-age appear in grey box. The evolution models are of a metal-free planet (dashed blue), hot diluted core of 25 M ⊕ (red), and of 5 M ⊕ core with high atmospheric opacity (green). Models emphasize the requirement of additional energy sources to explain the observed properties of TOI-4201b.  \n<!-- image -->  \na difficulty to fit observations with high accuracy without additional energy sources.  \nBased on the power-law relation from Thorngren et al. (2016), we compute the total heavy element mass required for TOI-4201b, and we obtain a high value around 100 M ⊕ under the assumption that the planet is not affected by inflation mechanisms. Apart from a typical 10 M ⊕ core occupied in Thorngren et al. (2016), the rest of the metals with a mass of about 90 M ⊕ are supposed to go directly into the gas envelope. However, the resulting low metal content (up to 3%) of TOI-4201b from our interior structure modeling makes this planet challenging in terms of planet formation theory. In this system, despite the high metallicity of TOI-4201, the giant planet appears to be a metal-poor planet and its bulk metallicity is at most identical to its star metallicity. This finding challenges the classical core-accretion trend of the stellar-planet metallicity (Thorngren et al. 2016). The formation path of TOI-4201b - a massive metal-poor planet around a metal-rich low-mass star remains unclear. Given the low metallicity of the planet, planet formation by gravitational instability cannot be ruled out (Boss 2011).', '5.3. Potential heating mechanisms': 'The discrepancy between the host star iron abundance and planet metallicity ( Z planet < 0 . 03 and Z star = 0 . 067 ± 0 . 012) could be slightly mitigated if additional heating processes have an effect. Such processes may  \ninject energy and slow down the cooling, leading to an inflated planet radius. Although the low equilibrium temperature ( < 1000 K) makes Ohmic dissipation less probable (Batygin & Stevenson 2010), we propose three possible pathways of planet heating.', '5.3.1. Tidal heating': "The short distance from the star suggests that tidal heating may introduce non-negligible energy into the planet's interior, as it went through the circularization of its orbit (Leconte et al. 2010). The low orbital eccentricity of 0 . 041 ± 0 . 015 also suggests that the tidal effect has an influence on the planet. We calculate the tidal circularization timescale of TOI-4201b using the equation from Goldreich & Soter (1966), which yields a τ e = 0 . 68 Gyr, assuming a typical planet quality factor Q p ∼ 10 6 of hot Jupiters. This is very close to the lower bound age estimation 0 . 7 -2 . 0 Gyr for the host star. If this is the case, it favors the younger age limit for the system and the tidal heating probably took place until recently.", '5.3.2. Gas giant merger': "Another possible scenario is a head-on collision of two less massive gas giants (Li et al. 2010; Liu et al. 2015), which will result in a coalescence of solid cores and gaseous envelopes without substantial material loss. Due to the high energy generated during the impact, the final merger remnant might be puffed up. This scenario also accounts for the high mass of the planet.  \nPlanet collision should be more likely around metalrich stars. With high solid density protoplanetary disks, they are likely to have multiple giant planets at the beginning. Such a configuration has also been found around M dwarfs, for example GJ 317, a metalrich ([Fe/H]=0 . 30 ± 0 . 08 dex) M star hosting two cold Jupiters (Anglada-Escud'e et al. 2012; Feng et al. 2020). Violent planet-planet interactions will probably then occur through gravitational perturbations (Rasio & Ford 1996). Under certain conditions, one of these formed giant planets will be scattered inward to large eccentricities (Dawson & Murray-Clay 2013). While in some extreme cases, the instability may cause close encounters if planets have sufficiently small spacings, which will result in a direct collision (Li et al. 2021).  \nAs a violent process, one would naturally expect to see a misalignment between the spin axis of the host star and the orbital angular momentum axis of the planet, induced by planet-planet scattering (Liu & Ji 2020). Observations have shown that several hot Jupiters have high obliquities (projected spin-orbit angle), which may be left during the interactions between planets (Winn &  \nFabrycky 2015). If TOI-4201b did undergo a dynamically hot formation history, we may be able to find clues through studying the stellar obliquity (Stef'ansson et al. 2022).", '5.3.3. Embryo capture by a gas giant': 'An alternative route to inject energy is through embryo capture. At first, there probably exists a proto gas giant outside the snow line accompanied with several rocky embryos such as super Earths in the inner disk. During the inward migration of the cold giant planet (Lin et al. 1996), these interior embryos would be trapped by the mean motion resonance (MMR) of the gas giant, leading to an orbital decay along the path. The planetary resonant perturbation probably then excites the orbits of these unstable embryos (Zhou & Lin 2007; Wang et al. 2021). If their eccentricity damping rate is fast, we will see a short-period giant accompanied by close-in super Earths. However, the orbits of rocky embryos may crisscross that of the migrating gas giant if the damping timescale is long. They will pass through the envelope of the proto gas giant and collides into its center core. Similar to the aforementioned pathway, the energy produced by the impact may inflate the gas giant.  \nA planetary system with both outer gas giants and inner small planets is not unusual. Huang et al. (2016) found that nearly half of warm Jupiters co-exist with low-mass planets. In addition, recent work from Wu et al. (2023) also reported that at least 12 ± 6% hot Jupiters as well as ≥ 70 ± 16% warm Jupiters have nearby planetary companions by investigating the transit timing variation signals, which probably support this hypothesis. Although the number is still limited, such a planetary configuration has also been found around M dwarfs through RV surveys (e.g., GJ 876; Rivera et al. 2010). Since the disk migration and embryo capture are quiescent, unlike the previously mentioned process, here we expect a well-aligned system with a low stellar obliquity. Such gentle gas disk migrations are also unlikely to produce large eccentricities of giant planets (Dunhill et al. 2013), like what is expected for TOI-4201b.', '5.4. Tidally induced evolution': "As TOI-4201b is a new addition to the growing number of giant planets orbiting M dwarfs (see AlvaradoMontes 2022), studying its tidal evolution could allow us to shed some light on the energy dissipation of M-dwarf hosts and put some constraints on stellar and planetary interiors. We study the tidal evolution of TOI-4201b under the formalism used by Alvarado-Montes (2022) to understand the orbital decay of giant planets orbiting M dwarfs due to the dissipation of inertial waves (IWs) in  \nFigure 13. Top panel: Semi-major axis as a function of time. Bottom panel: advance in mid-transit for TOI-4201b for an observation baseline of six years measured from today. For both panels, the colors represent different stellar rotation periods ( P rot ) from 10 to 40 days, as expected for M dwarfs (McQuillan et al. 2014). The orange line stands for the measured P rot in this work.  \n<!-- image -->  \nconvective envelopes and internal gravity waves (IGWs) in stellar radiative regions. For this particular system, IWs were not excited in the convective envelope due to the short orbital period of the planet and the slow rotation of the host star (i.e. P ≱ P rot 2 ; see Barker 2020).  \nIGWs in the stellar radiative regions are not expected to be damped due to wave breaking, as a planetary mass threshold (i.e. critical mass) M crit > 10 3 M J would be necessary for such a mechanism to occur for the given stellar age and stellar mass (see M crit values in Figure 9 in Barker 2020). However, if IGWs are somehow excited (e.g., via radiative diffusion), we predict a stellar quality factor Q ' = Q ' IGW ≈ 2 . 6 × 10 7 . This could also be the  \ncase in later stages of stellar evolution as M crit becomes smaller, thus producing fully damped IGWs due to wave breaking. For this tidal quality factor, we found that TOI-4201b would undergo orbital decay from its initial position to the Roche limit in a timescale of ∼ 8 Gyr for a 17.3-day stellar rotation period. The orbital decay for other stellar rotation periods is also shown in Figure 13 (top panel), and the migration timescales resemble those of other short-period giant planets undergoing orbital decay for P rot < 40 days (see e.g., Brown et al. 2011; Alvarado-Montes & Garc'ıa-Carmona 2019).  \nUsing the results from the tidal simulations performed for TOI-4201b (top panel in Figure 13), we computed the rate of change of the orbital period of TOI-4201b as ˙ P = 0 . 03185 +0 . 00049 -0 . 00044 ms yr -1 . Using these results, we calculated the advance in the mid-transit time of TOI-4201b (bottom panel in Figure 13) where it can be seen that an advance of ∼ 60 ms is expected for the next six years of observations, being three orders of magnitude smaller than those calculated for ultrashort-period Jupiters (see e.g., McCormac et al. 2019; Alvarado-Montes et al. 2021). It is worth noting that we studied the tidal effects for TOI-4201b assuming a co-planar system, so future observations (see next subsection) can help us refine our calculations by using a formalism where the evolution of the projected obliquity λ would also be included.", '5.5. Prospects for future observations': 'Based on the mass correlation from Thorngren et al. (2016), TOI-4201b is supposed to contain a large amount of metal in the envelope but we would expect to see a smaller planet radius in that case. Our interior structure model, instead, suggests the envelope is metalpoor. It motivates future atmospheric characterization using space telescopes like JWST to study the chemical composition in the planet atmosphere. Therefore, we compute the Transmission Spectroscopy Metric (TSM; Kempton et al. 2018), which is a quality factor of an object suitable for atmospheric composition studies. We find a TSM = 24 ± 6, much smaller than the threshold of 90 for high-quality Jovian-like planets recommended by Kempton et al. (2018), indicating that TOI-4201b is a challenging case but still possible with multi-visit observations.  \nIn addition, studying stellar obliquity may help us trace back the planet dynamic history. The host star shows a clear rotation signal around 17.3 days, corresponding to a stellar rotation speed of v sin i ∼ 1 . 8 km s -1 . To probe the opportunity of detecting the Rossiter-McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924) and measure the projected spin or-  \nt angle, we estimate the RM signal amplitude to a first order using  \nA RM ∼ 2 3 ( R p /R ∗ ) 2 √ 1 -b 2 × v sin i, (1)  \nwhere b is the impact parameter. The moderate impact parameter also eliminates the usual covariance between v sin i and the projected obliquity λ . Coupled with results from the joint fit, we find an A RM about 44 m/s. Although the host star is faint and the transit duration is short, the RM observation is still possible with ground spectroscopic facilities on large telescopes like MAROON-X and ESPRESSO.', '5.6. Metallicity preference of hot and warm Jupiters around M dwarfs': 'In order to compare the planet-metallicity dependence to probe the similarity of the planet formation history, we next investigate the metallicity distribution of four planet groups: HJG, WJG, HJM and WJM, corresponding to hot and warm Jupiters around G (0 . 90 ≤ M ∗ ≤ 1 . 06 M ⊙ ) and M dwarfs based on the discoveries retrieved from the NASA Exoplanet Archive (Akeson et al. 2013). To take the planet insolation into consideration, here we define the hot and warm Jupiter groups with a scaled semi-major axis ( a/R ∗ ) and planet mass boundary cut. We designate hot Jupiters as planets having mass above 0.3 M J and a/R ∗ < 20 while cold Jupiters as those with a/R ∗ > 20. Under this definition, our sample contains 103 HJG, 70 WJG, 9 HJM and 12 WJM. First, both hot and warm Jupiters are likely to form around metal-rich G-type stars. These two samples have a similar preference for stellar metallicity, both of which have a median [Fe/H] around 0.14 dex as shown in Figure 14. Interestingly, most hot Jupiters are orbiting M dwarfs with even higher super-solar metallicity ( ∼ 0.28 dex). Such a strong dependence between giant planets and stellar metallicity favors the core accretion formation mechanism (Santos et al. 2004; Fischer & Valenti 2005; Sousa et al. 2011; Wang et al. 2018). Meanwhile, the iron abundances of M dwarfs that host warm Jupiters have a much wider distribution, from -0 . 3 to 0.3 dex with a median of 0.06 dex. Therefore, the formation of warm Jupiters around M dwarfs is probably less sensitive to stellar metallicity, indicating a different formation pathway that is less dependent on metallicity. But it is still possible that this feature is biased due to the limited M dwarf sample size or different facilities and methods used to measure the iron abundance. Future observations and homogeneous spectroscopic analysis may help draw a firm conclusion.', '6. CONCLUSIONS': 'Figure 14. Planet mass ( M p or M p sin i ) versus host star metallicity. The red and blue dots are hot ( a/R ∗ < 20) and warm ( a/R ∗ ≥ 20) Jupiters around M dwarfs. The vertical red and blue dashed lines are the median metallicity of two samples: 0.28 and 0.06 dex. The background gray dots and triangles are hot and warm Jupiters around G-type stars under the same definition. Both of them have a median iron abundance of 0.14 dex, shown as a vertical black line. The top two panels are the histogram of these four samples.  \n<!-- image -->  \nIn this paper, we report the discovery and characterization of TOI-4201b, a massive and dense hot Jupiter transiting an early-M dwarf. TOI-4201b was first alerted as a planet candidate based on the TESS data, and the planetary nature was then confirmed through groundbased photometric, spectroscopic and imaging observations. TOI-4201b has a radius of 1 . 22 ± 0 . 04 R J with a mass of 2 . 48 ± 0 . 09 M J . It orbits the host M star every 3.58 days on a nearly circular orbit ( e = 0 . 041 ± 0 . 015). The bulk density ρ p = 1 . 82 ± 0 . 19 g cm -3 of TOI4201b makes it one of the most massive and densest hot Jupiters around M dwarfs. Although the host star has super-solar metallicity, interior structure modeling suggests that the TOI-4201b is metal-poor, which differs from the classical picture of positive stellar-planet metallicity correlation. Planet formation scenarios more typical for the core accretion model can be envisioned if additional energy sources, like tidal heating or giant impacts, operate to inflate the planet radius. Future  \nstudies on the planet atmosphere and stellar obliquity may shed light on its formation and evolution.  \nWe also compare the stellar metallicity distributions of four planet groups: hot and warm Jupiters around G and M stars. We find that the M dwarfs hosting hot Jupiters have a higher mean metallicity than G dwarfs, favoring the core accretion scenario. Warm Jupiters around M dwarfs, instead, show a weak dependence on stellar iron abundance, which perhaps indicates a different formation story.', '7. ACKNOWLEDGMENTS': "We thank Beibei Liu and Haochang Jiang for the useful discussions on planet formation and scattering. We are grateful to Coel Hellier for the insights regarding the WASP data.  \nThis work is partly supported by the National Science Foundation of China (Grant No. 12133005). This research uses data obtained through the Telescope Access Program (TAP), which has been funded by the TAP member institutes. The authors acknowledge the Tsinghua Astrophysics High-Performance Computing platform at Tsinghua University for providing computational and data storage resources that have contributed to the research results reported within this paper.  \nBased on observations obtained at the CanadaFrance-Hawaii Telescope (CFHT) which is operated from the summit of Maunakea by the National Research Council of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. The observations at the Canada-France-Hawaii Telescope were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site. Based on observations obtained with SPIRou, an international project led by Institut de Recherche en Astrophysique et Plan'etologie, Toulouse, France.  \nThis work is partly supported by the Natural Science and Engineering Research Council of Canada and the Institute for Research on Exoplanets through the Trottier Family Foundation. This work makes use of observations from the Las Cumbres Observatory global telescope network and it is partly supported by JSPS KAKENHI Grant Numbers JP17H04574, JP18H05439, and JST CREST Grant Number JPMJCR1761.  \nSome of the observations in this paper made use of the High-Resolution Imaging instrument Zorro and were obtained under Gemini LLP Proposal Number: GN/S2021A-LP-105. Zorro was funded by the NASA Exoplanet Exploration Program and built at the NASA Ames Research Center by Steve B. Howell, Nic Scott, Elliott P. Horch, and Emmett Quigley. Zorro was  \nmounted on the Gemini South telescope of the international Gemini Observatory, a program of NSF's OIR Lab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. on behalf of the Gemini partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigaci'on y Desarrollo (Chile), Ministerio de Ciencia, Tecnolog'ıa e Innovaci'on (Argentina), Minist'erio da Ciˆencia, Tecnologia, Inova¸c˜oes e Comunica¸c˜oes (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).  \nBased on observations obtained with the Samuel Oschin 48-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grant No. AST-1440341 and a collaboration including Caltech, IPAC, the Weizmann Institute for Science, the Oskar Klein Center at Stockholm University, the University of Maryland, the University of Washington, Deutsches Elektronen-Synchrotron and Humboldt University, Los Alamos National Laboratories, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.  \nThis work makes use of observations from the LCOGT network. Part of the LCOGT telescope time was granted by NOIRLab through the Mid-Scale Innovations Program (MSIP). MSIP is funded by NSF.  \nThis paper includes data gathered with the 6.5 meter Magellan Telescopes located at Las Campanas Observatory, Chile.  \nThe ULiege's contribution to SPECULOOS has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013) (grant Agreement n · 336480/SPECULOOS), from the Balzan Prize and Francqui Foundations, from the Belgian Scientific Research Foundation (F.R.S.-FNRS; grant n · T.0109.20), from the University of Liege, and from the ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. SPECULOOS-North has received financial support from the Heising-Simons  \nFoundation and from Dr. and Mrs. Colin Masson and Dr. Peter A. Gilman.  \nThe postdoctoral fellowship of KB is funded by F.R.S.-FNRS grant T.0109.20 and by the Francqui Foundation. This publication benefits from the support of the French Community of Belgium in the context of the FRIA Doctoral Grant awarded to MT. J.d.W. and MIT gratefully acknowledge financial support from the Heising-Simons Foundation, Dr. and Mrs. Colin Masson and Dr. Peter A. Gilman for Artemis, the first telescope of the SPECULOOS network situated in Tenerife, Spain. BVR thanks the Heising-Simons Foundation for support. MG is FNRS-F.R.S Research Director. KAC and SQ acknowledges support from the TESS mission via subaward s3449 from MIT. SQ acknowledges support from the TESS GI Program under award 80NSSC21K1056.  \nFunding for the TESS mission is provided by NASA's Science Mission Directorate. This research has made use of the Exoplanet Follow-up Observation Program website, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This paper includes data collected by the TESS mission, which are publicly available from the Mikulski Archive for Space Telescopes (MAST). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www. cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www. cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This work made use of tpfplotter by J. Lillo-Box (publicly available in www.github.com/jlillo/tpfplotter), which also made use of the python packages astropy , lightkurve , matplotlib and numpy .  \nFacilities: TESS, Gaia, CFHT/SPIRou, Magellan/MagE, Gemini/Zorro, LCOGT, MuSCAT, SPECULOOS, ZTF  \nSoftware: lightkurve (Lightkurve Collaboration et al. 2018), astropy (Astropy Collaboration et al. 2013, 2018), AstroImageJ (Collins et al. 2017), juliet (Espinoza et al. 2019), batman (Kreidberg 2015), radvel (Fulton et al. 2018), tpfplotter (Aller et al. 2020)", 'REFERENCES': "Akeson, R. L., Chen, X., Ciardi, D., et al. 2013, PASP, 125,  \n- Schlegel, D. J., Finkbeiner, D. P., & Davis, M. 1998, ApJ, 500, 525, doi: 10.1086/305772  \n```\nZechmeister, M., & Kurster, M. 2009, A&A, 496, 577, doi: 10.1051/0004-6361:200811296 Zhou, G., Huang, C. X., Bakos, G. ' A., et al. 2019, AJ, 158, 141, doi: 10.3847/1538-3881/ab36b5\n```  \n```\nZhou, J.-L., & Lin, D. N. C. 2007, ApJ, 666, 447, doi: 10.1086/520043\n```"}