A Massive Hot Jupiter Orbiting a Metal-rich Early M Star Discovered in the TESS Full-frame Images

Tianjun Gan; Charles Cadieux; Farbod Jahandar; Allona Vazan; Sharon X. Wang; Shude Mao; Jaime A. Alvarado-Montes; D. N. C. Lin; Étienne Artigau; Neil J. Cook; René Doyon; Andrew W. Mann; Keivan G. Stassun; Adam J. Burgasser; Benjamin V. Rackham; Steve B. Howell; Karen A. Collins; Khalid Barkaoui; Avi Shporer; Jerome de Leon; Luc Arnold; George R. Ricker; Roland Vanderspek; David W. Latham; Sara Seager; Joshua N. Winn; Jon M. Jenkins; Artem Burdanov; David Charbonneau; Georgina Dransfield; Akihiko Fukui; Elise Furlan; Michaël Gillon; Matthew J. Hooton; Hannah M. Lewis; Colin Littlefield; Ismael Mireles; Norio Narita; Chris W. Ormel; Samuel N. Quinn; Ramotholo Sefako; Mathilde Timmermans; Michael Vezie; Julien de Wit

doi:10.3847/1538-3881/acf56d

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 frequently than Sun-like 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).

The 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 before 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 goes up from 0.5 toward 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 owing to the cutoff of the inward flow of solids (Mulders et al. 2021).

On 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 mid- to 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 example, 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).

The Transiting Exoplanet Survey Satellite (TESS; 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 past 3 yr (e.g., HATS-71b, Bakos et al. 2020; HATS-74Ab and HATS-75b, Jordán et al. 2022; TOI-3629b and TOI-3714b, Cañas 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ñas et al. 2020). All of these findings gradually enable glimpses of their formation and evolution history.

In this paper, 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, and the comparison of stellar metallicity distribution of hot/warm Jupiters around G/M stars. We conclude with our findings in Section 6.

2. Observations

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 revisited 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, 2020b; Fausnaugh et al. 2020). Due to the faintness of TOI-4201 (Tmag = 13.5), the light curve was not examined initially because the TESS vetting team only inspects targets with Tmag ≤ 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 ≤ Tmag ≤ 13.5 that have QLP light curves.

However, 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 pixel 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 × 15 pixel 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.

**Figure 1.** Top panels: the 15 × 15 pixel TESS FFI from Sectors 6 and 33. The black shaded region represents the 3 × 3 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 with blue tickmarks.
Download figure:
Standard image High-resolution image

Combining 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érn-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 binaries and refine the transit ephemeris and 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 the details of all observations in Table 1 and describe each of them below.

Table 1. Ground-based Photometric Follow-up Observations for TOI-4201

Obs Date	Telescope	Filter	Pixel Scale (arcsec)	FWHM (arcsec)	No. of Observations	Coverage
2021-09-26	LCO-SAAO	$g^{\prime}$ and $i^{\prime}$	0.39	5.1	23 & 24	Ingress
2021-10-04	LCO-CTIO	$g^{\prime}$ and $i^{\prime}$	0.39	3.4	26 & 25	Full
2022-01-30	MuSCAT	$g^{\prime}$ and $r^{\prime}$ and $z^{\prime}$	0.36	3.0	124 & 301 & 299	Full
2023-02-04	SPECULOOS-North	$z^{\prime}$	0.35	1.0	519	Full

Download table as: ASCII Typeset image

2.2.1. LCOGT

We acquired two Sloan Digital Sky Survey (SDSS) $g^{\prime}$ - and $i^{\prime}$ -band alternating time-series observations using 1 m telescopes from the Las Cumbres Observatory Global Telescope network (LCOGT; Brown et al. 2013). The first partial transit (ingress-only) observation was done at the South African Astronomical Observatory (SAAO) on 2021 September 26, while the second one was obtained at the Cerro Tololo Inter-American 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^{\prime} \times 26^{\prime}$ and a plate scale of 0 farcs 389 pixel⁻¹. The exposure times of $g^{\prime}$ - and $i^{\prime}$ -band observations are 300 and 180 s, respectively. After the raw images were calibrated 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 farcs 5) and 16-pixel (6 farcs 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 January 30 with MuSCAT using exposure times of 120, 50, and 50 s for g, r, and z_s bands, respectively. MuSCAT is a multiband 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 a $6\buildrel{\,\prime}\over{.} 1\times 6\buildrel{\,\prime}\over{.} 1$ field of view, 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. Five comparison stars and aperture radii of 12 pixels (4 farcs 2) yield the optimum light curves for all bands.

2.2.3. SPECULOOS-North

We obtained a full transit of TOI-4201b by the SPECULOOS-North telescope on 2023 February 4 in the Sloan $z^{\prime}$ filter with an exposure time of 20 s. SPECULOOS-North is a 1.0 m Ritchey−Chrétien telescope equipped with a thermoelectrically cooled 2K × 2K Andor iKon-L BEX2-DD CCD camera with a pixel scale of 0 farcs 35 and a field of view of $12^{\prime} \times 12^{\prime}$ (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ártir in Baja California, México (Demory et al. 2020). Data reduction and photometric measurements were performed using the PROSE ³⁶ pipeline (Garcia et al. 2021). The resulting light curve is presented in Figure 2.

**Figure 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 in the lower left corner of each panel.
Download figure:
Standard image High-resolution image

2.3. Spectroscopic Follow-up

2.3.1. MagE

We 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 farcs 7, and we observed the target during nautical twilight. We used the 0 farcs 7 × 10'' slit to obtain resolution λ/Δλ ≈ 6000 over 3200–10000 Å. We collected two 290 s exposures at an air mass 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 and ThAr arc lamp and Xe flash and incandescent flat-field 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. 2020a, 2020b) using standard settings. The final calibrated spectrum has a median signal-to-noise ratio (S/N) of ≈100 in the 8000 Å region.

Figure 3 compares the MagE spectrum of TOI-4201 to the best-match spectral template constructed from SDSS 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épine 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 shows Hα in absorption with an equivalent width of 0.426 ± 0.017 Å and no detection of Li i, ruling out a young (≲10 Myr) source. Using the calibration of Lépine 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 SpectroPolarimètre InfraROUge (SPIRou; Donati et al. 2020) between 2022 November 6 and 2023 January 1. SPIRou is a near-infrared high-resolution (R ≈ 75,000) fiber-fed echelle spectrograph installed on the 3.6 m Canada–France–Hawaii Telescope (CFHT), covering a wavelength range from 0.98 to 2.5 μm. Since TOI-4201 is a faint target (Hmag > 11), instead of taking simultaneous drift calibration with the thermalized Farby–Pérot (FP) etalon, we opted to use the Dark mode to avoid contaminating science spectra. We conducted all of our observations under an environment of air mass around 1.3 and seeing about 0 farcs 6. Every night, we took two continuous sequences of TOI-4201 at a random phase, with a 1200 s exposure time for each one. For our 28 spectra, the median S/N per pixel in the middle of H band is about 40.

The 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 preprocessed 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 (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-step 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₂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 air mass and water column. The last step mitigates finite-resolution effects, i.e., the fact that deconvolution between stellar+telluric spectrum and 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 detail in E. Artigau et al. (2023, in preparation).

The RV measurements were obtained from the telluric-corrected spectra using the line-by-line (LBL) method of Artigau et al. (2022). In this framework, Doppler shifts are measured on small chunks of the spectrum called "lines" (∼16,000 lines for an M dwarf 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-1759b, Martioli et al. 2022; TOI-2136b, Gan et al. 2022b; TOI-1452b, Cadieux et al. 2022; TOI-1695b, Kiefer et al. 2023). After subtracting the systemic velocity of about 42.1 km s⁻¹, we list the LBL RVs of TOI-4201 in Table 2, which have a median uncertainty of 24 m s⁻¹.

Table 2. SPIRou Radial Velocities for TOI-4201

BJD	RV (m s⁻¹)	σ_RV (m s⁻¹)
2,459,889.94627	−343.4	29.7
2,459,891.05138	−225.4	22.6
2,459,891.06557	−177.0	22.8
2,459,891.93548	362.6	28.9
2,459,892.96892	42.1	28.2
2,459,892.98312	−42.1	25.4
2,459,894.00997	−516.7	24.5
2,459,894.02416	−526.7	25.4
2,459,895.02960	272.1	25.2
2,459,895.04379	306.3	25.1
2,459,895.96208	335.8	23.3
2,459,895.97627	362.0	23.1
2,459,899.05405	441.1	24.0
2,459,899.06824	415.7	24.1
2,459,900.00510	119.1	24.5
2,459,900.01929	89.8	26.3
2,459,900.93532	−341.4	32.9
2,459,904.99832	−327.4	30.0
2,459,915.92131	−263.0	23.6
2,459,915.93550	−289.3	23.2
2,459,924.91005	166.7	29.1
2,459,924.92423	160.6	29.3
2,459,943.93428	−467.7	24.0
2,459,943.94847	−422.0	23.6
2,459,944.84612	−79.2	23.0
2,459,944.86031	−69.8	22.9
2,459,945.91587	467.0	21.8
2,459,945.93007	467.7	22.0

Download table as: ASCII Typeset image

We then performed a generalized Lomb−Scargle (GLS) periodogram (Zechmeister & Kürster 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 frequencies of f = 1 + 1/3.58 day⁻¹ and f = 1 − 1/3.58 day⁻¹. After subtracting the best Keplerian model from the joint fit (see Section 4), we find that they disappear and no other significant peaks, which may point to the existence of a secondary planet, show up with false-alarm probability (FAP) below 0.1%.

**Figure 5.** Left panel: the time-series SPIRou RVs 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 the joint fit. The orange shaded region denotes the 1σ credible intervals of the RV model. Residuals are plotted below.
Download figure:
Standard image High-resolution image

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 owing 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 nondetection of small planets residing within 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).

TOI-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 (562 and 832 nm) with output data products including a reconstructed image with robust contrast limits on companion detections. Nine sets of 1000 × 0.06 s 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 mag below that of the star itself from the 8 m telescope diffraction limit (20 mas) out to 1 farcs 2. At the distance of TOI-4201 (d = 189 pc), these angular limits correspond to spatial limits of 3.8−227 au.

**Figure 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.
Download figure:
Standard image High-resolution image

3. Stellar Characterization

3.1. Stellar Parameters

We apply three different methods to determine the stellar properties of TOI-4201: SED fitting, spectroscopic analysis, and empirical relations. We take the weighted stellar radius, mass, and effective temperature and list them in Table 3.

Table 3. Summary of Stellar Parameters for TOI-4201

Parameter	Value	Ref.
Main identifiers:
TIC	95057860	TIC V8^[1]
Gaia ID	2997312063605005056	Gaia DR3^[2]
Equatorial coordinates:
α_J2015.5	06:01:53.93	TIC V8
δ_J2015.5	−13:27:40.93	TIC V8
Photometric properties:
TESS (mag)	13.501 ± 0.007	TIC V8
B (mag)	16.696 ± 0.139	APASS^[3]
V (mag)	15.297 ± 0.024	APASS
G (mag)	14.480 ± 0.003	Gaia DR3
G_BP (mag)	15.472 ± 0.003	Gaia DR3
G_RP (mag)	13.495 ± 0.004	Gaia DR3
J (mag)	12.258 ± 0.021	2MASS^[4]
H (mag)	11.564 ± 0.024	2MASS
K (mag)	11.368 ± 0.025	2MASS
W1 (mag)	11.272 ± 0.024	WISE^[5]
W2 (mag)	11.301 ± 0.021	WISE
W3 (mag)	11.283 ± 0.155	WISE
Astrometric properties:
ϖ (mas)	5.291 ± 0.018	Gaia DR3
μ_α (mas yr⁻¹)	11.731 ± 0.017	Gaia DR3
μ_δ (mas yr⁻¹)	6.052 ± 0.018	Gaia DR3
Stellar parameters:
Spectral type	M1.0 ± 0.5	This work
RV (km s⁻¹)	42.1 ± 0.3	This work
Distance (pc)	189.0 ± 0.6	This work
U_LSR (km s⁻¹)	−23.9 ± 0.2	This work
V_LSR (km s⁻¹)	−17.6 ± 0.2	This work
W_LSR (km s⁻¹)	6.0 ± 0.1	This work
M_* (M_⊙)	0.61 ± 0.02	This work
R_* (R_⊙)	0.63 ± 0.02	This work
ρ_* (g cm⁻³)	3.44 ± 0.35	This work
$\mathrm{log}{g}_{* }\ (\mathrm{cgs})$	4.64 ± 0.03	This work
L_* (L_⊙)	0.076 ± 0.004	This work
T_eff (K)	3794 ± 79	This work
[Fe/H]	0.30 ± 0.10	This work^[6]
[Fe/H]	0.22 ± 0.20	This work^[7]
[Fe/H]	0.52 ± 0.08	This work^[8]
P_rot (days)	17.3 ± 0.4	This work
Age (Gyr)	0.7 − 2.0	This work

Note. [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.

Download table as: ASCII Typeset image

3.1.1. SED Fitting

We 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. (2017), and Stassun & Torres (2018). We pulled the JHK_S magnitudes from the Two Micron All Sky Survey (2MASS; Cutri et al. 2003), the W1–W3 magnitudes from the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), the zy magnitudes from Pan-STARRS (Chambers et al. 2016), and the G, G_BP, and 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).

**Figure 7.** SED 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 best-fit PHOENIX atmosphere model (black). The Gaia spectrum is shown as a gray swathe; the inset focuses on the Gaia spectrum.
Download figure:
Standard image High-resolution image

We 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 χ² of 1.4. The metallicity in particular is further constrained by the Gaia spectrum (see gray swathe in Figure 7 and inset). Integrating the (unreddened) model SED gives the bolometric flux at Earth, F_bol = (6.80 ± 0.16) × 10⁻¹¹ erg s⁻¹ cm⁻². Taking the F_bol and T_eff together with the Gaia parallax gives the stellar radius, R_* = 0.63 ± 0.03 R_⊙.

We 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 use 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. 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 and then estimate R_* from the Stefan−Boltzmann relation and 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.

The planet host sits high on the color–magnitude diagram for a main-sequence M dwarf (Mann et al. 2015), suggesting that 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 [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 F. Jahandar et al. (2023, in preparation). 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. 2012b; 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. 2012a; Allard et al. 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 S/N 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 the 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 $\sqrt{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 but a 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.

**Figure 8.** SPIRou observation of OH and Fe i lines of TOI-4201 (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.
Download figure:
Standard image High-resolution image

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 TOI-4201 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).

Next, 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).

Finally, we obtain the stellar mass M_* using Equation (2) in Mann et al. (2019) according to the M_*–M_K relation.³⁷ 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 Equation (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 GLS periodogram (Zechmeister & Kürster 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- and 20.4-day signals have a false-alarm probability (FAP) below 0.1%.

**Figure 9.** The GLS 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-day rotational signal of TOI-4201. A zoomed-in plot is shown in the upper right corner.
Download figure:
Standard image High-resolution image

Due 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 the 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 deg², ZTF scans the entire Northern sky every 2 days, providing a large amount of data for time-domain science. TOI-4201 is observed by ZTF using two different CCDs (ID = 2 and 16) in both g and r bands. All of these observations were done with an exposure time of 30 s, 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 March 27 and 2022 December 27, (ii) 184 measurements from CCD-16 in g band between 2018 March 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 data sets, and we show the results in Figure 9. Except for the first one that shows no significant signals with FAP below 10%, all three of the other light curves have a peak around 17.34 days with an 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.

Based 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.

Building 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 of 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 three-dimensional velocity 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 from the SPIRou measurement. Taking the solar velocity components relative to the LSR (U_⊙, V_⊙, W_⊙) = (9.58, 10.52, 7.01) km s⁻¹ 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⁻¹. We further make use of the BANYAN Σ algorithm (Gagné 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 shows no evidence of cluster membership. We finally follow the procedure in Bensby et al. (2003, 2014) and estimate a low probability ratio of 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 RVs 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).

We 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 parameterization 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₁ and r₂ instead and allow them to 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³⁸ 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 jitter term for every photometric data set to account for additional white noise.

In terms of the RV modeling, we fit a standard Keplerian orbit, leaving $e\sin \omega$ and $e\cos \omega$ as free parameters. We place wide informative uniform priors on the RV semiamplitude K, the systematic velocity μ_SPIRou, and 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.

Table 4. Parameter Priors and Best-fit Values Along with the 68% Credibility Intervals in the Final Joint Fit for TOI-4201

Parameter	Prior	Value	Description
Orbit parameters:
P (days)	${ \mathcal N }$ (3.6, 0.1)	3.5819194 ± 0.0000011	Orbital period
T_c (BJD −2,457,000)	${ \mathcal N }$ (1470.9716, 0.1)	1470.9619 ± 0.0004	Mid-transit time
r₁	${ \mathcal U }$ (0, 1)	0.6662 ± 0.0109	Parameterization for p and b
r₂	${ \mathcal U }$ (0, 1)	0.1955 ± 0.0008	Parameterization for p and b
$e\sin \omega$	${ \mathcal U }$ (−1, 1)	−0.0402 ± 0.0155	Parameterization for e and ω
$e\cos \omega$	${ \mathcal U }$ (−1, 1)	0.0006 ± 0.0106	Parameterization for e and ω
K (m s⁻¹)	${ \mathcal U }$ (0, 600)	456.3 ± 12.5	RV semiamplitude
μ_SPIRou (m s⁻¹)	${ \mathcal U }$ (−500, 500)	−11.2 ± 8.8	Systemic velocity
σ_SPIRou (m s⁻¹)	${ \mathcal U }$ (0, 100)	45.2 ± 7.1	RV jitter

Stellar parameter:
ρ_* (kg m⁻³)	${ \mathcal L }{ \mathcal U }$ (10², 10⁵)	3786 ± 146	Stellar density

Dilution factors:
D_{TESS S06}	${ \mathcal T }{ \mathcal N }$ (0.9, 0.1², 0, 1)	0.85 ± 0.01
D_{TESS S33}	${ \mathcal T }{ \mathcal N }$ (0.9, 0.1², 0, 1)	0.85 ± 0.01
D_ground	1 (Fixed)	⋯

Limb-darkening coefficients:
q_{1,TESS S06}	${ \mathcal U }$ (0, 1)	0.16 ± 0.06
q_{2,TESS S06}	${ \mathcal U }$ (0, 1)	0.65 ± 0.23
q_{1,TESS S33}	${ \mathcal U }$ (0, 1)	0.32 ± 0.12
q_{2,TESS S33}	${ \mathcal U }$ (0, 1)	0.19 ± 0.12
q_LCO,SAAO,g	${ \mathcal U }$ (0, 1)	0.55 ± 0.27
q_LCO,SAAO,i	${ \mathcal U }$ (0, 1)	0.19 ± 0.12
q_LCO,CTIO,g	${ \mathcal U }$ (0, 1)	0.75 ± 0.09
q_LCO,CTIO,i	${ \mathcal U }$ (0, 1)	0.30 ± 0.08
q_MuSCAT,g	${ \mathcal U }$ (0, 1)	0.54 ± 0.06
q_MuSCAT,r	${ \mathcal U }$ (0, 1)	0.62 ± 0.03
q_MuSCAT,z	${ \mathcal U }$ (0, 1)	0.31 ± 0.04
q_SPECULOOS,z	${ \mathcal U }$ (0, 1)	0.29 ± 0.04

Photometric jitter:
σ_{TESS S06} (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${8.3}_{-8.3}^{+54.2}$
σ_{TESS S33} (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${0.1}_{-0.1}^{+32.7}$
σ_LCO,SAAO,g (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${15.7}_{-14.9}^{+43.6}$
σ_LCO,SAAO,i (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${31.8}_{-30.9}^{+99.5}$
σ_LCO,CTIO,g (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${2663.7}_{-476.9}^{+599.2}$
σ_LCO,CTIO,i (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${2000.5}_{-310.7}^{+374.3}$
σ_MuSCAT,g (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${12.4}_{-11.8}^{+19.6}$
σ_MuSCAT,r (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${38.8}_{-37.6}^{+48.3}$
σ_MuSCAT,z (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${2000.5}_{-12.5}^{+19.4}$
σ_SPECULOOS,z (ppm)	${ \mathcal L }{ \mathcal U }$ (10⁻⁶, 10⁵)	${2309.0}_{-169.4}^{+172.8}$

Derived parameters:
R_p/R_*	⋯	0.1955 ± 0.0008	Scaled planet radius
a/R_*	⋯	13.694 ± 0.151	Scaled semimajor axis
a (au)	⋯	0.040 ± 0.001	Semimajor axis
b	⋯	0.499 ± 0.016	Impact parameter
i (deg)	⋯	88.0 ± 0.1	Orbital inclination
e	⋯	0.041 ± 0.015	Orbital eccentricity
ω (deg)	⋯	130.0 ± 44.5	Argument of periapsis

Planetary physical parameters:
R_p (R_J)	⋯	1.22 ± 0.04	Planet radius
M_p (M_J)	⋯	2.48 ± 0.09	Planet mass
ρ_p (g cm⁻³)	⋯	1.82 ± 0.19	Planet density
${T}_{\mathrm{eq}}^{[1]}$ (K)	⋯	725 ± 20	Equilibrium temperature
S (S_⊕)	⋯	45.2 ± 3.2	Planetary insolation

Note. ${ \mathcal N }$ (μ, σ²) means a normal prior with mean μ and standard deviation σ. ${ \mathcal U }$ (a, b) stands for a uniform prior between a and b. ${ \mathcal L }{ \mathcal U }$ (a, b) is a log-uniform prior between a and b. ${ \mathcal T }{ \mathcal N }$ (μ, σ², a, b) represents a truncated normal prior ranging from a to b. [1] We do not consider heat distribution between the dayside and nightside here and assume albedo A_B = 0.

Download table as: ASCII Typeset image

5. Discussions

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⁻³ 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án 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⁻³ and a ≤ 0.1 au), an unpopulated area proposed by Gan et al. (2022a), making it an interesting object for theories of planet formation.

Figure 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⁻³, significantly higher than other observed gas giants around M dwarfs, as is shown in Figure 11. The high 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 and 25 M_⊕ in the interior of TOI-4201b, depending on the metal distribution in the interior, its history, and its energy state.

According 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 TOI-4201b under various conditions, based on Vazan et al. (2013) and Vazan et al. (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 10 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 metal-poor envelope.

In 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 (gray 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 the 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 10 times solar (green) barely fits the observed properties. These models indicate a difficulty to fit observations with high accuracy without additional energy sources.

Based 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 inject 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 nonnegligible 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⁶ of hot Jupiters. This is very close to the lower bound age estimation of 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.

Planet collision should be more likely around metal-rich 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 metal-rich ([Fe/H] = 0.30 ± 0.08 dex) M star hosting two cold Jupiters (Anglada-Escudé 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), whereas 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).

As 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 & Fabrycky 2015). If TOI-4201b did undergo a dynamically hot formation history, we may be able to find clues through studying the stellar obliquity (Stefànsson 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 collide into its center core. Similar to the aforementioned pathway, the energy produced by the impact may inflate the gas giant.

A 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 coexist with low-mass planets. In addition, recent work from Wu et al. (2023) also reported that at least 12% ± 6% of hot Jupiters and ≥70% ± 16% of 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 Alvarado-Montes 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 convective envelopes and internal gravity waves (IGWs) in stellar radiative regions. For this particular system, IWs were not excited in the convective envelope owing to the short orbital period of the planet and the slow rotation of the host star (i.e., $P\ngeqq \tfrac{{P}_{\mathrm{rot}}}{2};$ see Barker 2020).

IGWs in the stellar radiative regions are not expected to be damped owing to wave breaking, as a planetary mass threshold (i.e., critical mass) M_crit > 10³ 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^{\prime} =Q{{\prime} }_{\mathrm{IGW}}\approx 2.6\times {10}^{7}$ . This could also be the case 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).

Using 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 $\dot{P}\,={0.03185}_{-0.00044}^{+0.00049}$ ms yr⁻¹. Using these results, we calculated the advance in the mid-transit time of TOI-4201b (bottom panel of Figure 13), where it can be seen that an advance of ∼60 ms is expected for the next 6 yr of observations, being three orders of magnitude smaller than those calculated for ultra−short-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 coplanar 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 that the envelope is metal-poor. 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 multivisit observations.

In 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\sim 1.8$ km s⁻¹. To probe the opportunity of detecting the Rossiter−McLaughlin (RM) effect (Rossiter 1924; McLaughlin 1924) and measure the projected spin–orbit angle, we estimate the RM signal amplitude to a first order using

$\begin{eqnarray}&&{A}_{\mathrm{RM}}\sim \displaystyle \frac{2}{3}{\left({R}_{p}/{R}_{* }\right)}^{2}\sqrt{1-{b}^{2}}\times v\sin i,\end{eqnarray} \tag{ 1 }$

where 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 of 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_⊙ ≤ 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 semimajor 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 and 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 supersolar 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 owing 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.

**Figure 14.** Planet mass (M_p or ${M}_{p}\sin i$ ) vs. host star metallicity. The red and blue circles are hot (a/R_* < 20) and warm (a/R_* ≥ 20) Jupiters around M dwarfs, respectively. The vertical red and blue dashed lines are the median metallicity of two samples: 0.28 and 0.06 dex, respectively. The background gray circles and triangles are hot and warm Jupiters around G-type stars, respectively, 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 histograms of these four samples.
Download figure:
Standard image High-resolution image

**Figure 14.** Planet mass (M_p or ${M}_{p}\sin i$ ) vs. host star metallicity. The red and blue circles are hot (a/R_* < 20) and warm (a/R_* ≥ 20) Jupiters around M dwarfs, respectively. The vertical red and blue dashed lines are the median metallicity of two samples: 0.28 and 0.06 dex, respectively. The background gray circles and triangles are hot and warm Jupiters around G-type stars, respectively, 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 histograms of these four samples.
Download figure:
Standard image High-resolution image

6. Conclusions

In 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 ground-based 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⁻³ of TOI-4201b makes it one of the most massive and densest hot Jupiters around M dwarfs. Although the host star has supersolar 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 studies on the planet atmosphere and stellar obliquity may shed light on its formation and evolution.

We 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.

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.

This 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.

Based on observations obtained at the Canada–France–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étologie, Toulouse, France.

This 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 is partly supported by JSPS KAKENHI grant Nos. JP17H04574 and JP18H05439 and JST CREST grant No. JPMJCR1761.

Some of the observations in this paper made use of the High-Resolution Imaging instrument Zorro and were obtained under Gemini LLP Proposal No. GN/S-2021A-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 mounted 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ón y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).

Based 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−Milwaukee, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW.

This 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.

This paper includes data gathered with the 6.5 m Magellan Telescopes located at Las Campanas Observatory, Chile.

The 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 No. 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 Foundation and from Dr. and Mrs. Colin Masson and Dr. Peter A. Gilman.

The postdoctoral fellowship of K.B. 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 M.T. J.d.W. and M.I.T. 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. B.V.R. thanks the Heising-Simons Foundation for support. M.G. is FNRS-F.R.S Research Director. K.A.C. and S.Q. acknowledge support from the TESS mission via subaward s3449 from MIT. S.Q. acknowledges support from the TESS GI Program under award 80NSSC21K1056.

Funding 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.

Facilities: TESS - , Gaia - , CFHT - , Magellan:Baade - , Gemini:Gillett - , LCOGT - Las Cumbres Observatory Global Telescope, MuSCAT - , SPECULOOS - , ZTF - .

Software: 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).

A Massive Hot Jupiter Orbiting a Metal-rich Early M Star Discovered in the TESS Full-frame Images

Article metrics

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations

2.1. TESS

2.2. Ground-based Photometry

2.2.1. LCOGT

2.2.2. MuSCAT

2.2.3. SPECULOOS-North

2.3. Spectroscopic Follow-up

2.3.1. MagE

2.3.2. SPIRou

2.4. High Angular Resolution Imaging

3. Stellar Characterization

3.1. Stellar Parameters

3.1.1. SED Fitting

3.1.2. Spectroscopic Analysis

3.1.3. Empirical Relation

3.2. Stellar Rotation and Age Estimation

3.3. Stellar Kinematic Properties

4. Joint Fit

5. Discussions

5.1. TOI-4201b: A Massive and Dense Giant Planet

5.2. Interior Structure of TOI-4201b

5.3. Potential Heating Mechanisms

5.3.1. Tidal Heating

5.3.2. Gas Giant Merger

5.3.3. Embryo Capture by a Gas Giant

5.4. Tidally Induced Evolution

5.5. Prospects for Future Observations

5.6. Metallicity Preference of Hot and Warm Jupiters around M Dwarfs

6. Conclusions

Acknowledgments

Footnotes

A Massive Hot Jupiter Orbiting a Metal-rich Early M Star Discovered in the TESS Full-frame Images

Article metrics

Share this article

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction

2. Observations

2.1. TESS

2.2. Ground-based Photometry

2.2.1. LCOGT

2.2.2. MuSCAT

2.2.3. SPECULOOS-North

2.3. Spectroscopic Follow-up

2.3.1. MagE

2.3.2. SPIRou

2.4. High Angular Resolution Imaging

3. Stellar Characterization

3.1. Stellar Parameters

3.1.1. SED Fitting

3.1.2. Spectroscopic Analysis

3.1.3. Empirical Relation

3.2. Stellar Rotation and Age Estimation

3.3. Stellar Kinematic Properties

4. Joint Fit

5. Discussions

5.1. TOI-4201b: A Massive and Dense Giant Planet

5.2. Interior Structure of TOI-4201b

5.3. Potential Heating Mechanisms

5.3.1. Tidal Heating

5.3.2. Gas Giant Merger

5.3.3. Embryo Capture by a Gas Giant

5.4. Tidally Induced Evolution

5.5. Prospects for Future Observations

5.6. Metallicity Preference of Hot and Warm Jupiters around M Dwarfs

6. Conclusions

Acknowledgments

Footnotes