The Multiplanet System TOI-421: A Warm Neptune and a Super Puffy Mini-Neptune Transiting a G9 V Star in a Visual Binary^*

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2014) mission's primary scientific driver is to measure masses for transiting planets smaller than 4 R_⊕ around bright stars, in order to explore the transition from sub-Neptunes (with extended envelopes) to rocky planets (with compact atmospheres), that occurs at about 1.8 R_⊕. The Kepler mission revealed that small planets (especially in the super-Earth and sub-Neptune regime), in compact coplanar multiplanet systems, are very common (Latham et al. 2011; Lissauer et al. 2011, 2014; Rowe et al. 2014). While most of the Kepler stars are distant and faint, making radial velocity (RV) follow-up very difficult, the TESS mission is focused on the nearest and brightest stars, so that an intensive follow-up, as well as atmospheric characterization, can more easily be achieved. Those prospects are rather important, for example, for future space-based observations with James Webb Space Telescope (JWST).

Since 2018 July, TESS has been scanning the sky and performing a photometric search for planets transiting bright stars. In its primary mission, this survey will cover 26 sectors, each of them monitored for ∼27 days, with candidate alerts released almost every month. TESS is expected to detect ∼10,000 transiting exoplanets (Sullivan et al. 2015; Barclay et al. 2018; Huang et al. 2018b). More than 1000 planet candidates have been revealed, with dozens of confirmed planets so far (e.g., Brahm et al. 2019; Esposito et al. 2019; Díaz et al. 2020; Lendl et al. 2020), some of which are multiplanet systems (e.g., Gandolfi et al. 2018, 2019; Huang et al. 2018a; Dragomir et al. 2019; Günther et al. 2019; Quinn et al. 2019).

Multiplanet systems are prime targets for testing planetary formation and evolution theories. Orbiting the same star, they offer an opportunity to simplify the assumptions of initial conditions and compare planets with different sizes and compositions in the same system. Such systems are also interesting for transmission spectroscopy, which allows us to characterize planetary atmospheres and compare them at different levels of incident stellar flux.

In the present paper, we report on the discovery and characterization of a sub-Neptune (TOI-421 b) and a Neptune (TOI-421 c) transiting the bright star BD-14 1137 (TOI-421, TIC 94986319), observed by TESS in Sectors 5 and 6. The work presented here is part of the RV follow-up project carried out by the KESPRINT collaboration.⁷⁶ (e.g., Grziwa et al. 2016; Van Eylen et al. 2016; Gandolfi et al. 2017; Nowak et al. 2017; Barragán et al. 2018; Korth et al. 2019; Persson et al. 2019), which aims to confirm and characterize planet candidates from the K2 and TESS space missions.

This paper is organized as follows. In Section 2, we describe the observations carried out with different space- and ground-based facilities. These include TESS photometry, AO imaging, ground-based photometry, and high-resolution spectroscopy. In Section 3, we derive the stellar parameters of TOI-421. Section 4 reports the TESS photometry analysis with the detection of the transit signals. Section 6 presents the High Accuracy RV Planet Searcher (HARPS) frequency analysis to confirm the planetary nature of the transiting companions and investigate additional signals associated with the stellar activity. Sections 7 and 8 present a preliminary RV modeling and a joint analysis of TESS light-curve and RV data, respectively. We discuss our results in Section 9, including the simulations of possible atmospheric signals in different wavelength bands and of the atmospheric evolution of both planets. Finally, we draw our conclusion in Section 10.

2.1. TESS Photometry

TOI-421 (TIC 94986319)—whose identifiers, coordinates, proper motion, optical and infrared magnitudes, and fundamental parameters are listed in Table 1—was monitored in Sectors 5 and 6 of the TESS mission between 2018 November 15 and 2019 January 7. The target was imaged on CCD 3 of camera 2 in Sector 5, and on CCD 4 of camera 2 in Sector 6. A total of 34,622 photometric data points were collected, each with an exposure time of 2 minutes. The target light curves were processed by the Science Processing Operations Center (SPOC; Jenkins et al. 2016) data reduction pipeline. The Simple Aperture Photometry (SAP) was used in the SPOC pipeline to produce the light curves (Twicken et al. 2010; Morris et al. 2017), and the Presearch Data Conditioning (PDCSAP) algorithm was used to remove known instrumental systematics (Smith et al. 2012; Stumpe et al. 2012, 2014). For the transit detection and analysis presented in this work (Section 4), we downloaded TOI-421's PDCSAP light curves, which are publicly available at the Mikulski Archive for Space Telescopes web-page.⁷⁷

Table 1.  Main Identifiers, Equatorial Coordinates, Proper Motion, Parallax, Optical and Infrared Magnitudes, and Fundamental Parameters of TOI-421

Parameter	Value	Source
Main identifiers
TIC	94986319	ExoFOPa
BD	−14 1137	ExoFOP
TYC	5344-01206-1	ExoFOP
2MASS	J05272482-1416370	ExoFOP
Gaia DR2	2984582227215748864	Gaia DR2b
Equatorial coordinates, parallax, and proper motion
R.A. (J2000.0)	05h27m2483	Gaia DR2
Dec. (J2000.0)	−14°16'3705	Gaia DR2
π (mas)	13.3407 ± 0.0361	Gaia DR2
μα (mas yr−1)	−35.687 ± 0.046	Gaia DR2
μδ (mas yr−1)	50.450 ± 0.064	Gaia DR2
Optical and near-infrared photometry
TESS	9.2711 ± 0.006	TIC v8c
G	9.7778 ± 0.0002	Gaia DR2
Bp	10.2034 ± 0.0012	Gaia DR2
Rp	9.2265 ± 0.0012	Gaia DR2
B	10.735 ± 0.076	TIC v8
V	9.931 ± 0.006	TIC v8
J	8.547 ± 0.020	2MASSd
H	8.219 ± 0.033	2MASS
Ks	8.071 ± 0.018	2MASS
W1	8.058 ± 0.023	AllWISEe
W2	8.110 ± 0.020	AllWISE
W3	8.060 ± 0.021	AllWISE
W4	7.809 ± 0.168	AllWISE
Fundamental parameters
v sin i⋆( km s−1)	1.8 ± 1.0	This work
Teff (K)	${5325}_{-58}^{+78}$	This work
log g⋆ (cgs)	${4.486}_{-0.018}^{+0.025}$	This work
[Fe/H] (dex)	−0.02 ± 0.05	This work
M⋆ (M⊙)	${0.852}_{-0.021}^{+0.029}$	This work
R⋆ (R⊙)	0.871 ± 0.012	This work
Age (Gyr)	${9.4}_{-3.1}^{+2.4}$	This work
Distance (pc)	74.94 ± 0.58	This work
AV (mag)	${0.11}_{-0.08}^{+0.12}$	This work

Notes.

^a https://exofop.ipac.caltech.edu/. ^bGaia Collaboration et al. (2018). ^cStassun et al. (2018). ^dCutri et al. (2003). ^eCutri (2013).

Download table as:  ASCIITypeset image

2.2. Ground-based Photometry

2.3. AO Imaging

2.4. Ground-based Spectroscopy

We carried out spectroscopic follow-up observations of TOI-421 using different facilities—as described in the subsections below—to spectroscopically confirm the planetary nature of the transit signals detected in the TESS light curve and determine the masses of the two transiting planets.

2.4.1. FIES

We started the RV follow-up of TOI-421 in 2019 February using the FIbre-fed Echellé Spectrograph (FIES; Frandsen & Lindberg 1999; Telting et al. 2014) mounted at the 2.56 m Nordic Optical Telescope (NOT) of Roque de los Muchachos Observatory (La Palma, Spain). FIES is mounted inside an insulated building where the temperature is kept constant within 0.02 °C and fed with octagonal fibres to improve the RV stability of the spectrograph (Stürmer et al. 2018). We employed the intermediate-resolution fiber, which provides a resolving power of R = 45,000 over the wavelength range 3660–9275 Å, and acquired 10 spectra between 2019 February 2 and March 13. We used the same observing strategy as in Buchhave et al. (2010) and Gandolfi et al. (2013), i.e., we split the observations into three subexposures to remove cosmic-ray hits, and bracketed the three exposures with long-exposed (T_exp ≈ 80 s) ThAr spectra to trace the instrument drift and improve the wavelength solution. We reduced the data using IRAF and IDL standard procedures and extracted relative RV measurements by cross-correlating the observed Echellé spectra with the first epoch spectrum (Table 2).

Table 2.  FIES RV Measurements of TOI-421

BJDTDBa	RV	σ	Texp	S/Nb
−2450000	( km s−1)	( km s−1)	(s)
8517.532619	0.0000	0.0064	1800	65
8522.486513	−0.0070	0.0042	2400	88
8523.431757	−0.0078	0.0046	2400	84
8524.489067	−0.0126	0.0046	2400	82
8539.420470	−0.0054	0.0042	2400	89
8540.437340	−0.0146	0.0028	2700	120
8541.423092	−0.0140	0.0035	2400	94
8554.407973	−0.0064	0.0043	2400	88
8556.406509	−0.0152	0.0042	2700	87

Notes.

^aBarycentric Julian dates are given in barycentric dynamical time. ^bS/N ratio per pixel at 550 nm.

Download table as:  ASCIITypeset image

2.4.2. CORALIE

We also took seven high-resolution spectra (R ≈ 60,000) of TOI-421 using the CORALIE Echellé spectrograph on the Swiss 1.2 m Euler telescope at La Silla Observatories, Chile (Queloz et al. 2001). We extracted the RV measurements by cross-correlating the CORALIE Echellé spectra with a binary G2 mask (Pepe et al. 2002). The Doppler measurements show no significant RV variation at a level of ∼8 m s⁻¹ and exclude that TOI-421 is an eclipsing binary mimicking planetary transits.

2.4.3. HARPS

We acquired 105 high-resolution (R = 115,000) spectra of TOI-421 using the HARPS (Mayor et al. 2003) spectrograph mounted at the ESO-3.6 m telescope of La Silla Observatory, Chile. Installed in a pressure- and temperature-controlled enclosure and fed with octagonal fibres, HARPS has demonstrated a long-term precision at the 1 m s⁻¹ level and below (Lovis et al. 2006). Our HARPS observations were performed over two observing seasons between 2019 February and 2020 January, as part of the observing programs 1102.C-0923, 0103.C-0874, 0103.C-0759, 0103.C-0442, and 60.A-9709. We used the second fiber of the instrument to monitor the sky background and set the exposure time to 900–2100 s, depending on sky conditions and constraints of the observing schedule. We reduced the data using the dedicated HARPS Data Reduction Software (DRS) and computed the RVs by cross-correlating the Échelle spectra with a G2 numerical mask (Baranne et al. 1996; Pepe et al. 2002; Lovis & Pepe 2007). We also used the DRS to extract the line profile asymmetry indicators—namely, the FWHM and the bisector inverse slope (BIS) of the cross-correlation function (CCF), and the Ca ii H and K lines activity indicator⁷⁸ $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$. We report the HARPS RV measurements and their uncertainties, along with the FWHM, BIS, $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$, exposure time, and signal-to-noise (S/N) ratio per pixel at 550 nm in Table 3.

Table 3.  HARPS RV Measurements of TOI-421

BJDTDBa	RV	σRV	BIS	FWHM	$\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$	σ $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$	Texp	S/Nb
−2450000	( km s−1)	( km s−1)	( km s−1)	( km s−1)			(s)
8528.596454	79.5473	0.0011	−0.0181	6.7449	−4.948	0.019	1200	69
8529.643890	79.5417	0.0013	−0.0198	6.7420	−4.842	0.018	1200	61
8530.589677	79.5415	0.0011	−0.0257	6.7420	−4.978	0.019	1200	73
8540.601897	79.5363	0.0010	−0.0224	6.7489	−5.001	0.022	1200	83
8578.531222	79.5516	0.0023	−0.0220	6.7481	−5.016	0.061	900	41
8580.540429	79.5490	0.0014	−0.0238	6.7417	−4.942	0.033	900	64
8581.569558	79.5492	0.0015	−0.0219	6.7370	−5.045	0.049	900	62
8709.911463	79.5529	0.0007	−0.0221	6.7564	−4.879	0.007	1800	104
8711.894944	79.5463	0.0010	−0.0279	6.7561	−4.920	0.011	1800	83
8713.921013	79.5440	0.0010	−0.0258	6.7547	−4.893	0.011	1700	81
8714.911893	79.5458	0.0010	−0.0191	6.7591	−4.883	0.011	1800	82
8715.921283	79.5436	0.0012	−0.0228	6.7568	−4.870	0.013	1800	67
8716.915892	79.5401	0.0011	−0.0198	6.7504	−4.877	0.012	1600	72
8717.927169	79.5408	0.0008	−0.0235	6.7513	−4.871	0.008	1800	97
8718.883686	79.5417	0.0009	−0.0262	6.7578	−4.912	0.010	1800	87
8721.904488	79.5465	0.0010	−0.0214	6.7569	−4.911	0.012	1800	81
8723.878421	79.5460	0.0008	−0.0231	6.7494	−4.916	0.009	1800	92
8724.913234	79.5496	0.0014	−0.0303	6.7434	−4.956	0.021	1500	58
8725.875540	79.5514	0.0014	−0.0309	6.7440	−4.919	0.019	1800	58
8726.862398	79.5498	0.0012	−0.0285	6.7390	−4.970	0.018	1200	67
8734.857658	79.5470	0.0009	−0.0235	6.7571	−4.868	0.010	2100	84
8736.885722	79.5499	0.0010	−0.0246	6.7596	−4.892	0.011	2100	81
8737.873716	79.5403	0.0013	−0.0185	6.7641	−4.954	0.021	1800	61
8738.871421	79.5412	0.0011	−0.0172	6.7524	−4.897	0.013	1800	74
8740.862468	79.5464	0.0012	−0.0214	6.7562	−4.889	0.016	1800	66
8741.837012	79.5468	0.0009	−0.0200	6.7482	−4.895	0.011	1800	84
8744.821435	79.5414	0.0010	−0.0254	6.7488	−4.909	0.012	1800	79
8745.817358	79.5393	0.0016	−0.0231	6.7471	−5.004	0.031	1800	52
8746.804032	79.5414	0.0012	−0.0261	6.7438	−4.990	0.020	1800	67
8747.835235	79.5389	0.0013	−0.0233	6.7452	−4.916	0.018	1800	63
8748.860025	79.5376	0.0009	−0.0264	6.7440	−4.928	0.011	1800	84
8750.820767	79.5414	0.0010	−0.0311	6.7503	−4.924	0.012	1800	76
8752.867278	79.5418	0.0011	−0.0203	6.7542	−4.920	0.013	1800	73
8753.779744	79.5395	0.0018	−0.0233	6.7559	−5.010	0.034	1800	49
8754.812264	79.5417	0.0011	−0.0278	6.7501	−4.892	0.011	1800	75
8755.830524	79.5492	0.0008	−0.0256	6.7512	−4.908	0.008	1620	101
8756.882374	79.5524	0.0008	−0.0249	6.7509	−4.900	0.009	1800	94
8757.807038	79.5465	0.0019	−0.0316	6.7564	−4.968	0.036	1800	46
8760.807608	79.5425	0.0013	−0.0232	6.7515	−4.931	0.018	1800	62
8761.806384	79.5452	0.0009	−0.0256	6.7525	−4.884	0.009	1800	90
8762.847560	79.5450	0.0009	−0.0279	6.7565	−4.895	0.011	1800	84
8763.827086	79.5402	0.0009	−0.0232	6.7547	−4.903	0.010	1800	84
8766.857081	79.5479	0.0009	−0.0195	6.7589	−4.895	0.011	1800	85
8767.723508	79.5481	0.0011	−0.0254	6.7614	−4.889	0.014	1800	72
8767.743659	79.5491	0.0012	−0.0248	6.7619	−4.864	0.013	1800	69
8780.805190	79.5336	0.0008	−0.0271	6.7448	−4.924	0.010	1800	96
8781.847532	79.5365	0.0009	−0.0299	6.7484	−4.942	0.012	1800	90
8782.800562	79.5381	0.0017	−0.0305	6.7489	−4.925	0.023	1800	51
8784.840817	79.5395	0.0013	−0.0303	6.7384	−5.043	0.023	1800	61
8785.732151	79.5386	0.0014	−0.0322	6.7431	−4.947	0.019	1800	57
8785.863860	79.5420	0.0016	−0.0373	6.7482	−4.916	0.021	1800	50
8791.866359	79.5440	0.0010	−0.0291	6.7541	−4.908	0.013	1800	79
8792.789899	79.5466	0.0009	−0.0231	6.7469	−4.920	0.010	1800	90
8793.813558	79.5425	0.0014	−0.0241	6.7486	−4.965	0.019	1800	60
8794.779193	79.5371	0.0013	−0.0228	6.7494	−4.934	0.018	1800	61
8795.703042	79.5318	0.0014	−0.0265	6.7443	−4.943	0.019	1800	58
8796.819006	79.5386	0.0017	−0.0264	6.7501	−4.992	0.029	1800	50
8797.761303	79.5437	0.0019	−0.0324	6.7361	−4.961	0.029	1800	46
8798.828187	79.5355	0.0021	−0.0316	6.7468	−4.967	0.033	1800	42
8798.849647	79.5430	0.0020	−0.0302	6.7529	−4.982	0.033	1800	45
8799.804743	79.5415	0.0016	−0.0213	6.7485	−4.935	0.022	1800	52
8802.684679	79.5488	0.0007	−0.0260	6.7498	−4.880	0.006	1800	111
8803.723755	79.5537	0.0008	−0.0257	6.7537	−4.873	0.007	1800	101
8805.801094	79.5514	0.0009	−0.0266	6.7587	−4.890	0.011	1800	87
8806.812098	79.5532	0.0009	−0.0205	6.7557	−4.897	0.010	1800	92
8820.826561	79.5478	0.0008	−0.0262	6.7460	−4.954	0.012	1800	102
8821.721745	79.5424	0.0007	−0.0243	6.7479	−4.964	0.010	1800	110
8823.814612	79.5494	0.0008	−0.0267	6.7466	−5.000	0.016	1800	100
8824.842097	79.5463	0.0010	−0.0259	6.7576	−4.948	0.015	1800	86
8826.728971	79.5352	0.0007	−0.0247	6.7491	−4.936	0.008	1800	123
8830.721095	79.5404	0.0008	−0.0259	6.7451	−4.924	0.010	1800	99
8831.736359	79.5368	0.0009	−0.0244	6.7450	−4.943	0.012	1800	90
8832.741251	79.5387	0.0011	−0.0260	6.7508	−4.894	0.014	1800	74
8833.766142	79.5419	0.0009	−0.0248	6.7506	−4.923	0.011	1800	95
8834.754203	79.5460	0.0008	−0.0280	6.7481	−4.896	0.010	1800	97
8835.678212	79.5481	0.0009	−0.0247	6.7447	−4.920	0.012	1800	87
8835.804299	79.5456	0.0007	−0.0245	6.7512	−4.914	0.009	1800	119
8841.709874	79.5455	0.0017	−0.0244	6.7541	−4.941	0.026	1500	50
8843.694063	79.5398	0.0013	−0.0252	6.7475	−4.906	0.018	1500	62
8845.701891	79.5434	0.0009	−0.0212	6.7444	−4.916	0.013	1500	90
8847.726300	79.5392	0.0008	−0.0224	6.7489	−4.923	0.012	1500	97
8856.746648	79.5476	0.0012	−0.0297	6.7456	−4.967	0.022	1800	70
8856.768244	79.5443	0.0014	−0.0242	6.7371	−4.969	0.027	1800	64
8857.687817	79.5417	0.0008	−0.0245	6.7430	−4.947	0.011	1800	107
8857.709216	79.5415	0.0007	−0.0254	6.7424	−4.958	0.011	1800	110
8858.570410	79.5403	0.0010	−0.0234	6.7414	−4.956	0.014	1800	77
8858.689630	79.5424	0.0009	−0.0248	6.7450	−4.933	0.013	1800	93
8859.582685	79.5422	0.0018	−0.0185	6.7456	−4.962	0.033	1800	49
8859.675147	79.5433	0.0011	−0.0281	6.7448	−4.966	0.016	1800	72
8860.603211	79.5408	0.0010	−0.0237	6.7412	−4.948	0.013	1800	85
8860.695071	79.5420	0.0010	−0.0277	6.7480	−4.979	0.017	1800	80
8861.558044	79.5435	0.0011	−0.0266	6.7424	−4.961	0.016	1800	73
8861.720386	79.5431	0.0009	−0.0250	6.7457	−4.939	0.014	1800	90
8862.575813	79.5388	0.0008	−0.0257	6.7449	−4.955	0.011	1800	94
8862.688991	79.5341	0.0009	−0.0294	6.7392	−4.953	0.013	1800	94
8863.592677	79.5362	0.0010	−0.0253	6.7433	−4.959	0.015	1800	87
8863.684930	79.5385	0.0008	−0.0243	6.7416	−4.936	0.012	1800	109
8864.587945	79.5413	0.0010	−0.0257	6.7473	−4.931	0.013	1800	80
8864.681911	79.5418	0.0009	−0.0235	6.7408	−4.926	0.012	1800	90
8865.614228	79.5448	0.0010	−0.0248	6.7429	−4.928	0.013	1800	81
8865.717198	79.5444	0.0010	−0.0180	6.7384	−4.954	0.015	1800	86
8866.667986	79.5441	0.0008	−0.0222	6.7457	−4.935	0.012	1800	100
8868.665966	79.5436	0.0009	−0.0225	6.7426	−4.949	0.013	1650	93
8869.699467	79.5441	0.0007	−0.0246	6.7449	−4.990	0.012	1660	109
8871.661192	79.5486	0.0010	−0.0245	6.7466	−4.967	0.016	1800	78

Notes.

^aBarycentric Julian dates are given in barycentric dynamical time. ^bS/N ratio per pixel at 550 nm.

Download table as:  ASCIITypeset images: 1 2

2.4.4. HIRES

Between 2019 September 17 and 2020 March 3, we obtained 28 spectra (R = 55,000) of TOI-421 over 27 nights using the High Resolution Échelle Spectrometer (HIRES; Vogt et al. 1994) on the 10 m Keck I telescope. The spectra were collected using an iodine cell for wavelength reference. The median exposure time was 680 s, which allowed us to achieve a S/N ratio of 200 per reduced pixel at 5500 Å. We also obtained a high S/N ratio template spectrum without the iodine cell, which was used as input for the standard forward-modeling procedures of the California Planet Search (Howard et al. 2010). We report the RVs and uncertainties based on the weighted mean and weighted error in the mean of the ∼700 individual spectral chunks in Table 4.

Table 4.  HIRES RV Measurements of TOI-421

BJDTDBa	RV	σ	Texp	S/Nb
−2450000	( km s−1)	( km s−1)	(s)
8744.056124	−0.0047	0.0011	770	218
8777.027485	0.0047	0.0011	592	218
8788.071425	0.0061	0.0011	680	221
8794.976150	−0.0076	0.0009	716	220
8796.022921	−0.0090	0.0011	774	221
8797.052515	−0.0061	0.0010	900	212
8798.104180	−0.0030	0.0010	706	218
8798.917940	−0.0062	0.0013	900	186
8802.899539	0.0038	0.0010	537	219
8809.050906	0.0086	0.0012	650	218
8815.908118	−0.0055	0.0012	808	216
8819.967016	0.0069	0.0010	900	214
8827.942590	−0.0138	0.0011	595	218
8832.964946	−0.0009	0.0012	877	221
8833.929016	0.0011	0.0012	582	220
8844.885674	−0.0048	0.0011	634	219
8845.918812	−0.0050	0.0011	651	217
8852.868230	0.0035	0.0012	633	220
8855.804440	0.0072	0.0012	900	197
8856.831220	0.0043	0.0011	900	182
8857.840778	−0.0011	0.0011	824	220
8869.855111	0.0026	0.0011	774	219
8870.880081	0.0055	0.0011	762	219
8878.833553	−0.0047	0.0011	567	219
8879.783275	−0.0060	0.0011	865	220
8880.810664	−0.0012	0.0011	582	219
8884.813904	0.0038	0.0011	512	219
8885.935404	0.0021	0.0015	899	164
8903.828998	0.0020	0.0012	900	174
8905.785063	−0.0007	0.0010	631	218
8906.756423	−0.0017	0.0011	670	218
8907.751496	−0.0039	0.0012	900	187
8911.799720	−0.0028	0.0010	668	216

Notes.

^aBarycentric Julian dates are given in barycentric dynamical time. ^bS/N ratio per pixel at 550 nm.

Download table as:  ASCIITypeset image

2.4.5. PFS

TOI-421 was also selected as a high-priority target by the Warm gIaNts with tEss (WINE) collaboration, which aims at systematically characterizing long-period transiting giant and Neptune-size planets from the TESS mission (see, e.g., Brahm et al. 2019; Jordán et al. 2020). In this context, we monitored TOI-421 with the Planet Finder Spectrograph (PFS; Crane et al. 2006, 2008, 2010) mounted on the 6.5 m Magellan II Clay Telescope at Las Campanas Observatory (LCO) in Chile. The observations were performed on nine different nights, between 2019 February 18 and October 11. For these observations, we used the 03 × 25 slit, which delivers high-resolution spectra with a resolving power of R = 130,000. The observing strategy consisted of obtaining two consecutive 1200 s spectra per night of TOI-421 in order to increase the total S/N ratio per epoch and also to average out the stellar and instrumental jitter. These observations were performed with the use of an iodine cell for determining the RV of the star. The PFS data were processed with a custom IDL pipeline that is capable of delivering RVs with a precision less than 1 m s⁻¹ (Butler et al. 1996). Additionally, three consecutive 1200 s iodine-free exposures of TOI-421 were obtained to construct a stellar spectral template for computing the RVs. The PFS RVs are listed in Table 5.

Table 5.  PFS RV Measurements of TOI-421

BJDTDBa	RV	σ	Texp	S/Nb
−2450000	( km s−1)	( km s−1)	(s)
8592.508995	−0.0006	0.0009	2400	73
8708.926166	0.0059	0.0009	1200	61
8717.901009	−0.0017	0.0010	1200	50
8738.854892	−0.0002	0.0012	1200	51
8739.858020	−0.0007	0.0010	1200	57
8763.842388	−0.0023	0.0008	1200	66
8764.859431	0.0011	0.0007	1200	71
8767.887884	0.0057	0.0007	1200	77
8768.852571	0.0042	0.0009	900	55

Notes.

^aBarycentric Julian dates are given in barycentric dynamical time. ^bS/N ratio per pixel at 550 nm.

Download table as:  ASCIITypeset image

The fundamental stellar parameters of TOI-421 were independently determined by three different methods—one based on spectral synthesis, one on template matching, and one using different sets of isochrones, as described in the paragraphs below. For this purpose, we utilized the co-added HARPS spectrum, which has an S/N ratio of ∼700 per pixel at 5500 Å.

Method 1. We derived the effective temperature T_eff, surface gravity log g_⋆, iron abundance [Fe/H], and projected rotational velocity v sin i_⋆ with the spectral analysis package Spectroscopy Made Easy (SME,version 5.2.2; Valenti & Piskunov 1996; Piskunov & Valenti 2017). We used the ATLAS12 model spectra (Kurucz 2013) and the line lists from the Vienna Atomic Line Database.⁷⁹ (Ryabchikova et al. 2015) to model the co-added HARPS spectrum. Following the modeling procedure detailed in Fridlund et al. (2017) and Persson et al. (2018), we measured T_eff from the wings of H-α line, and log g_⋆ from the line wings of the Ca and Mg triplets around 6100 Å and 5100 Å, respectively. We derived v sin i_⋆ and [Fe/H] from narrow unblended iron lines between 6000 and 6600 Å. The micro- and macro-turbulent velocities, v_mic and v_mac, were kept fixed using the Bruntt et al. (2010) and Doyle et al. (2014) calibration equations for Sun-like stars to 0.9 km s⁻¹ and 2.5 km s⁻¹, respectively. The final best-fitting model spectrum was checked using the Na doublet at 5888 and 5895 Å. We found T_eff = 5194 ± 60 K, log g_⋆ = 4.45 ± 0.05(cgs), [Fe/H] = −0.04 ± 0.06 (dex), and v sin i_⋆ = 1.8 ± 1.0 km s⁻¹. The uncertainties are internal error bars that do not account for the systematic uncertainties of the atmospheric models.

We measured the visual interstellar extinction (A_V) along the line of sight to TOI-421 following the method described in Gandolfi et al. (2008). Briefly, we built the spectral energy distribution (SED) of the star from the broadband photometry listed in Table 1, and fitted the SED using synthetic magnitudes computed from a low-resolution BT-NextGen model spectrum (Allard et al. 2012) with the same spectroscopic parameter as the star. We adopted the Cardelli et al. (1989) extinction law and assumed a total-to-selective extinction ratio of R_V = (A_V)/E(B − V) = 3.1. We found that the interstellar extinction along the line of sight is consistent with zero (A_V = 0.03 ± 0.03), as expected given the relatively short distance to the star (∼75 pc).

In order to compute the stellar mass, radius, and age we employed the Bayesian PARAM 1.3.⁸⁰ model tool tracks (da Silva et al. 2006) with the PARSEC isochrones (Bressan et al. 2012). Input parameters are T_eff and [Fe/H] from SME, the apparent visual magnitude, and the parallax (Table 1). We added 0.1 mas in quadrature to the parallax uncertainty of 0.0361 mas to account for systematic errors of Gaia's astrometry (Luri et al. 2018). The resulting stellar parameters are M_⋆ = 0.86 ± 0.02 M_⊙, R_⋆ = 0.86 ± 0.02 R_⊙, log g_⋆ = 4.48 ± 0.02 (cgs), and an age of 9.2 ± 2.3 Gyr. We note that the derived log g_⋆ is in very good agreement with the spectroscopic value.

Method 2. We used SpecMatch-Emp (Yee et al. 2017) to analyze the co-added HARPS spectrum via comparison with a spectroscopic library of well-characterized stars, enabling empirical estimates of the effective temperature, stellar radius, and photospheric iron content. Following the procedure described in Hirano et al. (2018) to adapt the HARPS spectrum for use with SpecMatch-Emp, we obtained T_eff = 5337 ± 110 K, R_⋆ = 0.972 ± 0.100 R_⊙, and [Fe/H] = 0.05 ± 0.09 dex.

Method 3. We also computed a set of uniformly inferred stellar parameters using isochrones (Morton 2015) and MIST (Choi et al. 2016) to fit Two Micron All Sky Survey (2MASS) JHKs photometry (Skrutskie et al. 2006) and Gaia Data Release 2 (DR2) parallax (Gaia Collaboration et al. 2016, 2018), adding 0.1 mas of uncertainty in quadrature to account for systematics (Luri et al. 2018). We assumed priors on T_eff and [Fe/H] based on the spectroscopic results from SME, using a more conservative prior width of 100 K for T_eff to account for systematic errors, and we sampled the posterior using MultiNest (Feroz et al. 2019). We obtained the following parameter estimates: T_eff = ${5325}_{-58}^{+78}$ K, log g_⋆ = ${4.486}_{-0.018}^{+0.025}$, [Fe/H] = −0.02 ± 0.05 dex, M_⋆ = ${0.852}_{-0.021}^{+0.029}$ ${M}_{\odot }$, R_⋆ = 0.871 ± 0.012 R_⊙, age = ${9.4}_{-3.1}^{+2.4}$ Gyr, distance = 74.94 ± 0.58 pc, and A_V = ${0.11}_{-0.08}^{+0.12}$.

We note that the stellar parameter estimates obtained using these three independent methods agree well within 1–1.5σ, and they serve to ensure the accuracy of our results.

Since the derived parameters we find using the Method 3 agree more with those found in Method 2, we adopted the results of Method 2 and list the final adopted parameter estimates in Table 1. We do so because Method 1 relies on the profile of the Balmer Hα in order to determine the T_eff, and this method become less accurate for cooler temperatures. The quoted projected rotational velocity v sin i_⋆ has, however, been derived following Method 1, because Method 2 does not provide it. The spectroscopic parameters of TOI-421 translate into a spectral type and luminosity class of G9 V (Pecaut & Mamajek 2013).

We also looked for solar-like oscillations (García & Ballot 2019) using an optimized aperture for asteroseismology (L. González-Cuesta et al. 2020, in preparation) for both sectors. While the probability of detection computed following Schofield et al. (2019) is very low for this star, we applied the A2Z pipeline (Mathur et al. 2010) on the concatenated light curves of Sectors 5 and 6. Some excess of power and periodicity of the modes (also known as the large frequency separation) was found around 2000 and 3000 μHz but with a very low confidence level. The expected frequency of maximum power being around 3500 μHz and the S/N still being very low, the seismic analysis is not conclusive.

The detection of a 16 day transit signal was issued by the TESS Science Office QLP pipeline in Sector 5, and was subsequently identified in the SPOC pipeline (Twicken et al. 2018) in the Sector 6 data set. The SPOC Data Validation difference image centroid offsets for TOI-421 in the multi-sector run indicated that the source of the transit signature was within 14 of the proper motion–corrected location of the target star. The detection was then released as a planetary candidate via the TOI releases portal⁸¹ on 2019 February 8. We independently performed transit searches on the PDCSAP light curve using the DST algorithm (Cabrera et al. 2012). The variability in the light curve was filtered using the Savitzky–Golay method (Savitzky & Golay 1964; Press et al. 2002), and a transit model described by a parabolic function was used for transit searches. The algorithm recovered the detection of TOI-421 c where the transit signal has a period of 16.069 ± 0.002 days, a transit depth of 2735.90 ± 76.04 ppm, and a duration of 2.72 ± 0.05 hr.

The detection of a 5.2 day signal in the follow-up RV data (see Section 6) prompted further analysis of the TESS light curve. The DST algorithm further detected a transit signal with a period of 5.197 ± 0.003 days, a depth of 654.21 ± 88.70 ppm, and a duration of 1.23 ± 0.14 hr. The detection of both transit signals were also confirmed with the software package EXOTRANS (Grziwa & Pätzold 2016), which applies the box-fitting least squares algorithm (Kovács et al. 2002) for transit searches. Figure 1 shows the PDCSAP light curve of TOI-421 along with the detection of the two planets.

Zoom In Zoom Out Reset image size

In order to check whether the 5.2 and 16 day transit signals could arise from another source, as well as to assess the dilution level of TOI-421, we visually inspected archival images and compared the positions of Gaia DR2 (Gaia Collaboration et al. 2018) sources with the SPOC photometric apertures from Sectors 5 and 6. We used the coordinates of TOI-421 from the TESS Input Catalog⁸² (Stassun et al. 2018) to retrieve Gaia DR2 sources and a 3' × 3' image from DSS2.⁸³ Following the procedures described in Gandolfi et al. (2019), we computed a photometric dilution level of 1.8 ± 0.4% for TOI-421, which is a little smaller than the SPOC contamination ratios of 0.028 and 0.024 in Sectors 5 and 6, respectively. This small difference is most likely due to the fact that the SPOC target star is affected by an artifact in the 2MASS catalog caused by a diffraction spike from the 2MASS telescope secondary spider. In our specific case, this creates a TICv8 neighbor 69 north of the target star, which does not exist. We note that our dilution calculation is not affected by this issue because it is based on Gaia DR2.

We discovered that the nearby fainter star, spatially located at ∼294 NW of TOI-421 (Gaia ID 2984582227215748224, ΔG = 4.8, indicated by an orange circle in Figure 2), has a parallax and Gaia G-band extinction that are consistent within the error bars with those of TOI-421, and that the two stars have similar proper motions. We concluded that the pair very likely form a visual binary and TOI-421 is the primary component. According to Gaia DR2 effective temperature (T_eff = ${3676}_{-385}^{+376}$ K), the secondary component is an M dwarf. The angular separation and parallax imply a separation between the two stars of about 2200 au.

Zoom In Zoom Out Reset image size

The Keck AO observations show no additional stellar companions were detected to within a resolution ∼0053 FWHM (Figure 3). The sensitivities of the final combined AO image were determined by injecting simulated sources azimuthally around the primary target every 45° at separations of integer multiples of the central source's FWHM (Furlan et al. 2017). The brightness of each injected source was scaled until standard aperture photometry detected it with 5σ significance. The resulting brightness of the injected sources relative to the target set the contrast limits at that injection location. The final 5σ limit at each separation was determined from the average of all of the determined limits at that separation; the uncertainty on the 5σ limit was set by the rms dispersion of the azimuthal slices at a given radial distance. The sensitivity curve is shown in Figure 3 along with an inset image zoomed to the primary target, showing no other companion stars.

Zoom In Zoom Out Reset image size

We also fitted the FIES, HARPS, HIRES, and PFS measurements using the same RV models presented in Table 6 (Section 7), including an RV linear trend to account for the presence of a potential outer companion. We found an acceleration consistent with zero within less than 1σ.

Table 6.  RV Model Selection

Model	AICc	BIC	Nfree	Ndata	rmsa	$\mathrm{ln}{ \mathcal L }$ b
2eS	587.83	633.88	19	123	2.25	−260.96
2cS	591.48	629.18	15	123	2.38	−268.33
2eGP	596.78	644.79	20	123	1.57	−267.24
2cGP	601.06	640.92	16	123	1.69	−275.19
2c	630.25	661.16	12	123	2.78	−282.70
2e	634.49	674.36	16	123	2.71	−279.66

Notes.

^aRoot mean square of the data minus the model. ^bLog-likelihood of the data given the model.

Download table as:  ASCIITypeset image

In order to search for the Doppler reflex motion induced by the 16 day transiting planet candidate and unveil the presence of possible additional signals induced by other orbiting planets and/or stellar activity, we performed a frequency analysis of the RV measurements, as well as of the Ca ii activity index ($\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$) and CCF asymmetry indicators (FWHM and BIS). Toward this aim, we used only the HARPS measurements, as they form the largest homogeneous data set among our spectroscopic data, and we analyzed only the 98 HARPS measurements collected between 2019 August and 2020 January, in order to avoid the presence of spurious peaks associated with the 1 yr sampling in the power spectrum of the HARPS time series.

The generalized Lomb–Scargle periodogram (Zechmeister & Kürster 2009) displays a significant peak at the frequency of the transit signal reported by the TESS QLP pipeline (f_c ≈ 0.062 days⁻¹, P_c ≈ 16.1 days). Following the bootstrap method (Murdoch et al. 1993; Hatzes 2016), we assessed its false alarm probability (FAP) by computing the GLS periodogram of 10⁵ mock time-series obtained by randomly shuffling the Doppler measurements and their uncertainties, while keeping the time stamps fixed. We found that the peak at f_c has an FAP ≪ 0.1%. No significant peak at f_c (Figure 4, fourth, fifth, and sixth panels) appears in the periodogram of the Ca ii activity index $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$ nor those of the CCF asymmetry indicators (FWHM and BIS), providing strong evidence that this Doppler signal is due to the stellar reflex motion induced by the transiting planet TOI-421 c detected in the TESS light curve.

Zoom In Zoom Out Reset image size

The periodogram of the 98 HARPS RVs also shows a significant peak (FAP < 0.1%) at f_b ≈ 0.193 days⁻¹ (P_b = 5.2 days) whose power increases once the Doppler signal of TOI-421 c is removed⁸⁴ (Figure 4; first and second panel). This second peak has no counterpart in the periodograms of $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$, FWHM, and BIS, indicating that it is not due to stellar activity. A reanalysis of the TESS light curve unveils the presence of a transit signal at 5.2 days, as described in Section 4, confirming that the Doppler signal discovered in the RV time series is associated with a second transiting planet (TOI-421 b).

The periodogram of the RV residuals following the subtraction of the Doppler signal induced by TOI-421 c (Figure 4; second panel) shows also a second significant peak⁸⁵ (FAP < 0.1%) at f_rot = 0.024 ± 0.003 days⁻¹, corresponding to a period of P_rot = ${42}_{-5}^{+6}$ d and an RV semi-amplitude variation of ∼2.4 m s⁻¹, whose power becomes stronger once the Doppler signal of TOI-421 b is also removed (third panel). This peak is also significantly detected in the GLS periodogram of $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$ (FAP < 0.1%; fourth panel). It is also found in the periodogram of the FWHM (fifth panel), although with a higher FAP ≈ 0.2%. This suggests that the rotation period of the star is close to ∼42 days and that the third Doppler signal at f_rot is induced by intrinsic stellar variability associated with the presence of active regions rotating on and off the visible stellar disk. We note that the periodograms of $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$ and FWHM (Figure 4, fourth and fifth panels) also display a peak at the first harmonic of the rotation period (2f_rot), which is likely due to the presence of active regions at opposite stellar longitudes.

In order to further investigate the nature of the ∼42 day signal, we searched the WASP-South data (Section 2.2.2) for any rotational modulation using the methods from Maxted et al. (2011). We found no significant periodicity, with a 95% confidence upper limit on the amplitude of 1 mmag. We did find a significant periodicity compatible with the lunar cycle, but this was seen only in the data from the 85 mm lenses, which are more vulnerable to moonlight, and not in the 200 mm data; furthermore, in the 85 mm data, the same signal was also seen in adjacent field stars, so we concluded that it is not intrinsic to TOI-421.

With a mean $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$ of −4.93 ± 0.04 (Table 3), TOI-421 is a relatively quiet star. The lack of significant rotational modulation in the WASP-South light curve might be explained if the spot-induced variability of TOI-421 is too low to produce a photometric signal with an amplitude higher than the WASP-South photometric precision. We further investigated this scenario with the code SOAP2.0 (Dumusque et al. 2014), which simulates stellar activity using a fine grid to model the photosphere of a spotted rotating star. For each grid cell, SOAP2.0 simulates the local CCF using as a reference the solar HARPS CCF, and it accounts for the contribution of spots and plages using the HARPS CCF of a solar active region. The adoption of of the solar CCF is an advantage for our test, because TOI-421 is a G9 star, meaning that it is more likely to behave like our Sun. For a given set of stellar parameters and spots/plages distribution, size, and temperature, SOAP2.0 can estimate the photometric and Doppler signals induced by active regions, accounting for the inhibition of the convective blueshift and limb darkening/brightening effects.

For TOI-421, we used the effective temperature T_eff and stellar radius R_⋆ listed in Table 1, and we assumed a rotation period of P_rot = 42 days. To account for the wavelength range covered by the HARPS spectra, we adopted the linear and quadratic limb darkening coefficients of the Sun (q₁ = 0.29 and q₂ = 0.34; Claret & Bloemen 2011). Assuming a simplified model with one single spot, we found that the ∼2.4 m s⁻¹ RV semi-amplitude variation induced by stellar activity could be accounted for by a typical sunspot, with a temperature contrast with respect to the quiet photosphere of ΔT = 663 K (Meunier et al. 2010), a radius of 0.10 R_⋆, and placed at a latitude of 30°, which is the average active latitude for the Sun (Donati & Landstreet 2009; Strassmeier 2009). The corresponding photometric variation would have an amplitude of 5000 ppm, equivalent to ∼5 mmag, which is higher than the 1 mmag upper limit of the WASP-South time series of TOI-421. We performed a similar test by replacing the spot with a plage. We assumed a temperature contrast of ΔT = 251 K (typical plage contrast for the Sun; see Meunier et al. (2010)) and a radius of 0.12 R_⋆, while we kept latitude identical to that of the spot (30°). While the resulting RV semi-amplitude is ∼2.4 m s⁻¹, the photometric variation is 500 ppm, which corresponds to ∼0.5 mmag, i.e., lower than the 1 mmag upper limit on the amplitude observed in the WASP-South photometry. In Figure 5, we show the results of our simulations. The upper panel displays the effect of the simulated spot on both the RV and the photometric signal, whereas the bottom panel shows the same for a plage. The photometric effect induced by a plage is thus one order of magnitude smaller than that caused by a spot. These results agree with those reported in Dumusque et al. (2014) and more recently by Shapiro et al. (2016) and Milbourne et al. (2019)

Zoom In Zoom Out Reset image size

To understand if the Doppler signal of TOI-421 is spot- or plage-dominated, we produced an average contour map of the HARPS CCF residuals (after the division for the average CCF), plotted versus RV and stellar rotation phase. The results are displayed in Figure 6. Positive deviations (i.e., cool star spots) are shown in red, while negative deviations (hot spots) are shown in blue. These can account for the RV variation due to stellar activity, if we consider their associated perturbation to be ΔRV ≃ 2 × FWHM × ΔI × f ≃ 14 × f, where ΔI ∼ 0.004 is the intensity range and f ≤ 1 the filling factor (Carleo et al. 2020). The contours show that the activity of the star is dominated by plages, though some spots are also evident.

Zoom In Zoom Out Reset image size

We conclude that the activity-induced Doppler signal of TOI-421 can be explained by plages that would also account for the nondetection of any photometric variability in the WASP-South time series. Alternatively, the star was photometrically quieter at the time of the WASP-South observations (2008–2015) and more active during our spectroscopic follow-up (2019–2020). This is corroborated by the fact that the contour map of the HARPS CCF also shows the presence of spots.

We also analyzed the TESS data to look for surface rotation modulation in the light curve due to the passage of spots or active regions on the visible disk. We did the analysis with two different light curves. The first is the one described in Section 2.1, and the second is the one optimized for asteroseismology following González-Cuesta et al. (2020, in preparation). The optimized aperture was obtained by selecting larger and larger apertures with different thresholds in the flux, starting from the SPOC aperture (the smaller one) and ranging up to the larger one with a threshold of 10e-/s, with increments of 10e-/s. In the resultant light curve, only points with a quality flag equal to zero have been retained. Missing points have been interpolated using inpainting techniques as in García et al. (2014b). The light curve has also been corrected from outliers and stitched together following García et al. (2011). The optimization is done by comparing the signal measured at the expected region for the modes (around 3500 μHz, Section 3) and the high-frequency noise above 2000 μHz. The best aperture found was the one with a threshold of 80 e-/s for Sector 5, and the one with a threshold of 10 e-/s for Sector 6. For both light curves, we removed the transits and concatenated them. Then we applied our rotation pipeline following García et al. (2014a), Mathur et al. (2014), and Santos et al. (2019). Our methodology consists of performing a time–frequency analysis with wavelets, computing the autocorrelation function (ACF) and a composite spectrum (CS) that is a combination of the first two methods (see Ceillier et al. 2017). We found a signal at 45 ± 3.54 days via the three methods. The heights of the ACF and CS are above the thresholds defined in Ceillier et al. (2017) to reliably select a rotation period. However, we usually require that we observe three rotation periods in order to more reliably assess the rotation period. With only 58 days of observations, we cannot not fulfill this criteria. However, having this period obtained with these three different methods and using two different processing of the light curves (in terms of apertures) suggests that this period could be from stellar origin and it is independent of the processing of the light curve or the aperture selected. The HARPS RV analysis also finds a rotation period of ∼42 days, lending more weight to the likelihood of it being real, and an analysis of TESS data complements the spectroscopic analysis. Indeed, with the TESS photometry alone, we could not be confident enough about the period found, as it could still be a harmonic of a longer periodicity or even something of instrumental origin.

Motivated by the results of our frequency analysis, we performed a series of fits to the RV data in order to enable model selection and obtain system parameter estimates. Specifically, we used RadVel.⁸⁶ (Fulton et al. 2018) to test six different two-planet models: circular orbits ("2c"), eccentric orbits ("2e"), circular orbits with a Gaussian Process (GP) noise model ("2cGP"), eccentric orbits with a GP noise model ("2eGP"), circular orbits with an additional sine curve for the stellar activity ("2cS"), and eccentric orbits with an additional sine curve for the stellar activity ("2eS"). We used a GP model with a quasi-periodic kernel, which has been used extensively in the literature to model stellar RV signals (see, e.g., Haywood et al. 2014; Grunblatt et al. 2015; Dai et al. 2017); to avoid overfitting, we imposed wide Gaussian priors on the hyperparameters loosely informed by our frequency analysis.

To compare the quality of these models, we computed both the commonly used Bayesian Information Criterion (BIC) and the Akaike Information Criterion (AICc; corrected for small sample sizes), which is a second-order estimator of information loss. The results of these fits are presented in Table 6, sorted in ascending order of AICc (best to worst). The 2eS model is strongly favored over the other models by the AICc, suggesting that the orbits of the two transiting planets are significantly eccentric and that the stellar activity is described reasonably well by a sinusoid. We note that the BIC presents a slight preference for the 2cS model, but the AICc has been suggested to have practical advantages over the BIC (Burnham & Anderson 2004). We performed a full MCMC exploration of the parameter space of the 2eS model using RadVel, yielding semi-amplitudes of ${K}_{{\rm{b}}}={4.33}_{-0.35}^{+0.37}$ and ${K}_{{\rm{c}}}={3.05}_{-0.34}^{+0.35}$ m s⁻¹ for the planetary components, and 2.52 ± 0.36 m s⁻¹ for the stellar component; the eccentricities are significant at the ∼2σ level and constrained to ${e}_{{\rm{b}}}={0.147}_{-0.065}^{+0.069}$ and ${e}_{{\rm{c}}}={0.171}_{-0.086}^{+0.087}$.

We performed a global analysis of our RV and transit data with the open-source code pyaneti ⁸⁷ (Barragán et al. 2019). Briefly, pyaneti⁸⁸ creates posterior distributions of the fitted parameters using a Markov Chain Monte Carlo (MCMC) sampling approach.

Transits are modeled following a Mandel & Agol (2002) quadratic limb darkening model, implemented with the parameterization suggested by Kipping (2013). A preliminary fit of the light curve shows that the limb darkening coefficients are not well-constrained for the LCOGT data (we note that the limb darkening coefficients for the TESS data are constrained by data themselves). Therefore, we used the code LDTk (Husser et al. 2013; Parviainen & Aigrain 2015) to estimate the limb darkening coefficients corresponding to a star with the stellar parameters presented in Table 1, as well as the Pan-STARSS Y-band (940–1060 nm). We used a Gaussian prior with mean given by the LDTk estimation and a conservative standard deviation of 0.1. We also note that we sampled for the stellar density ρ_⋆ and recover a/R_⋆ for each planet using Kepler's third law; for more details, see, e.g., Winn (2010). We implemented a Gaussian prior on the stellar density using the stellar mass and radius provided in Section 3 to constrain better the orbit eccentricities (see, e.g., Van Eylen & Albrecht 2015). For the transit analysis, we did not include the potential contamination of nearby stars because the expected effect is smaller than the white noise of the light curve (see Section 5). We note that the combination of this transit analysis together with the RV data provides a stronger constraint on the orbital eccentricities. The final analysis supports the conclusions of Section 7 (see Table 7).

Table 7.  TOI-421 Parameters

Parameter	Priora	Valueb
Model Parameters for TOI-421 b
Orbital period Porb (days)	${ \mathcal U }[5.1917,5.2017]$	5.19672 ± 0.00049
Transit epoch T0 (BJD—2450000)	${ \mathcal U }[8441.2335,8441.3335]$	${8441.2847}_{-0.0018}^{+0.0020}$
$\sqrt{e}\sin {\omega }_{\star }$	${ \mathcal U }(-1,1)$	${0.27}_{-0.21}^{+0.15}$
$\sqrt{e}\cos {\omega }_{\star }$	${ \mathcal U }(-1,1)$	$-{0.268}_{-0.090}^{+0.117}$
Scaled planetary radius Rp/ R⋆	${ \mathcal U }[0,0.1]$	${0.0285}_{-0.0018}^{+0.0019}$
Impact parameter, b	${ \mathcal U }[0,1.1]$	${0.933}_{-0.024}^{+0.016}$
RV semi-amplitude variation K (m s−1)	${ \mathcal U }[0,50]$	2.97 ± 0.27
Model Parameters for TOI-421 c
Orbital period Porb (days)	${ \mathcal U }[16.0642,16.0741]$	16.06819 ± 0.00035
Transit epoch T0 (BJD—2450000)	${ \mathcal U }[8440.0804,8440.1804]$	${8440.13162}_{-0.00068}^{+0.00070}$
$\sqrt{e}\sin {\omega }_{\star }$	${ \mathcal U }(-1,1)$	${0.348}_{-0.086}^{+0.065}$
$\sqrt{e}\cos {\omega }_{\star }$	${ \mathcal U }(-1,1)$	$-{0.164}_{-0.078}^{+0.084}$
Scaled planetary radius Rp/ R⋆	${ \mathcal U }[0,0.1]$	${0.0542}_{-0.0010}^{+0.0011}$
Impact parameter, b	${ \mathcal U }[0,1.1]$	${0.738}_{-0.032}^{+0.029}$
RV semi-amplitude variation K (m s−1)	${ \mathcal U }[0,50]$	4.66 ± 0.29
Model Parameters of Activity-induced RV Sinusoidal Signal
Period Prot (days)	${ \mathcal U }[35,50]$	${43.24}_{-0.55}^{+0.57}$
Epoch T0 (BJD—2450000)	${ \mathcal U }[8412.2835,8452.2835]$	${8430.33}_{-4.98}^{+4.77}$
RV semi-amplitude variation K (m s−1)	${ \mathcal U }[0,50]$	${2.36}_{-0.3}^{+0.3}$
Other System Parameters
Stellar density ρ⋆	${ \mathcal N }[1.91,0.14]$	1.93 ± 0.13
Systemic velocity γ ( km s−1)c	${ \mathcal U }[79.0318,80.0537]$	${79.54382}_{-0.00025}^{+0.00024}$
Instrumental systemic velocity γ ( km s−1)c	${ \mathcal U }[-0.5,0.5]$	$-{0.0169}_{-0.0015}^{+0.0015}$
Instrumental systemic velocity γ ( km s−1)c	${ \mathcal U }[-0.5,0.5]$	${0.00222}_{-0.00093}^{+0.00091}$
Instrumental systemic velocity γ ( km s−1)c	${ \mathcal U }[-0.5,0.5]$	$-{0.00096}_{-0.00043}^{+0.00041}$
RV jitter term σHARPS ( m s−1)	${ \mathcal U }[0,100]$	${1.88}_{-0.18}^{+0.20}$
RV jitter term σFIES ( m s−1)	${ \mathcal U }[0,100]$	${0.85}_{-0.65}^{+1.41}$
RV jitter term σPFS ( m s−1)	${ \mathcal U }[0,100]$	${2.44}_{-0.64}^{+1.00}$
RV jitter term σHIRES ( m s−1)	${ \mathcal U }[0,100]$	${2.11}_{-0.34}^{+0.40}$
Limb darkening q1 TESS	${ \mathcal U }[0,1]$	${0.269}_{-0.083}^{+0.121}$
Limb darkening q2 TESS	${ \mathcal U }[0,1]$	${0.65}_{-0.35}^{+0.24}$
Limb darkening q1 Pan-STARSS Y-band	${ \mathcal N }[0.24,0.1]$	${0.164}_{-0.097}^{+0.101}$
Limb darkening q2 Pan-STARSS Y-band	${ \mathcal N }[0.36,0.1]$	${0.348}_{-0.100}^{+0.101}$
Derived Parameters for TOI-421 b
Planet mass (M⊕)	⋯	7.17 ± 0.66
Planet radius (R⊕)	⋯	${2.68}_{-0.18}^{+0.19}$
Planet density (g cm−3)	⋯	${2.05}_{-0.41}^{+0.52}$
Semimajor axis a (au)	⋯	0.0560 ± 0.0018
e	⋯	${0.163}_{-0.071}^{+0.082}$
ω⋆ (deg)	⋯	${128.9}_{-27.2}^{+24.9}$
Orbital inclination i (deg)	⋯	${85.68}_{-0.46}^{+0.36}$
Transit duration (hr)	⋯	${1.107}_{-0.063}^{+0.065}$
Equilibrium temperatured Teq (K)	⋯	${981.4}_{-15.8}^{+16.3}$
Insolation Fp (F⊕)	⋯	${154.57}_{-9.73}^{+10.53}$
Derived Parameters for TOI-421 c
Planet mass (M⊕)	⋯	${16.42}_{-1.04}^{+1.06}$
Planet radius (R⊕)	⋯	${5.09}_{-0.15}^{+0.16}$
Planet density (g cm−3)	⋯	${0.685}_{-0.072}^{+0.080}$
Semimajor axis a (au)	⋯	0.1189 ± 0.0039
e	⋯	0.152 ± 0.042
ω⋆ (deg)	⋯	${114.7}_{-13.3}^{+15.6}$
Orbital inclination i (deg)	⋯	${88.353}_{-0.084}^{+0.078}$
Transit duration (hr)	⋯	${2.71}_{-0.038}^{+0.043}$
Equilibrium temperatured Teq (K)	⋯	${673.6}_{-10.9}^{+11.2}$
Insolation Fp (F⊕)	⋯	${34.32}_{-2.16}^{+2.34}$

Notes.

^a ${ \mathcal U }[a,b]$ refers to uniform priors between a and b, ${ \mathcal N }[a,b]$ to Gaussian priors with median a and standard deviation b, and ${ \mathcal F }[a]$ to a fixed value a. ^bParameter estimates and corresponding uncertainties are defined as the median and 68.3% credible interval of the posterior distributions. ^cHARPS RVs are absolute measurements: the HARPS γ velocity corresponds to systemic velocity. FIES, PFS, and HIRES RVs are relative measurements. ^dAssuming albedo equal to zero.

Download table as:  ASCIITypeset images: 1 2

The RV model was chosen following the results presented in Section 7. We used two Keplerian orbits to model the Doppler reflex motion induced by the two transiting planets, as well as an extra sinusoidal to take into account the activity signal induced by the star at its rotation period.

We used 500 Markov chains to explore the parameter space. We stopped the sampling once all chains converged; following Gelman & Rubin (1992), we define convergence when R < 1.02 for all parameters. The posterior distribution was created using the last 2500 converged iterations and the 500 chains, leading to a posterior of 1,250,000 points per parameter.

Details on the fitted parameters, adopted priors, and parameter estimates are given in Table 7. Inferred parameters are defined as the median and 68% region of the credible interval of the posterior distributions for each fitted parameter. Figure 7 shows the RV model time series. The phase-folded RV and transit plots are shown in Figures 8 and 9, respectively. We note that there is an apparent shift of ∼20 minutes on the LCO transit in Figure 9. Given that the expected TTVs in the system are smaller (see Section 9.4), it is likely that this effect is caused by systematics in the ground-based data.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The innermost planet, TOI-421 b (P_b = 5.2 days), has a mass of M_b = ${16.42}_{-1.04}^{+1.06}$ M_⊕ and a radius of R_b = ${2.68}_{-0.18}^{+0.19}$ R_⊕, yielding a density of ρ_b = ${2.05}_{-0.41}^{+0.52}$ g cm⁻³. The outer transiting planet, TOI-421 c (P_c ≈ 16.1 days), has a mass of M_c = ${16.42}_{-1.04}^{+1.06}$ M_⊕ and a radius of R_c = ${5.09}_{-0.15}^{+0.16}$ R_⊕ resulting in a mean density of ρ_c = ${0.685}_{-0.072}^{+0.080}$ g cm⁻³. Figure 10 shows the position of TOI-421 b and c in the mass–radius diagram along with the sample of small planets (R_p ≤ 6 R_⊕) whose masses and radii have been measured with a precision better than 20%. Given their positions with respect to theoretical mass–radius relations, both planets are expected to host an atmosphere dominated by light elements, namely H and He. We performed a series of simulations—including hydrogen escape rate (Section 9.1), planetary atmospheric evolution (Section 9.2), and retrievals of the transmission spectrum (Section 9.3)—which indicate that TOI-421 b and c are very intriguing planets for atmospheric characterization. Finally, we compute the helium 10830 Å simulation following the approach in Oklopčić & Hirata (2018), finding no significant helium absorption (≲0.5% at line center) for either of the two planets.

Zoom In Zoom Out Reset image size

9.1. Hydrogen Escape

Atmospheric escape in close-in planets takes place when the high-energy (X-ray+EUV; hereafter XUV) stellar photons photoionize and heat up the planetary upper atmospheres (e.g., Murray-Clay et al. 2009). The close distances of TOI-421 b and c to the star suggest that their atmospheres may be significantly heated by the stellar high-energy emission and hence are undergoing escape. Therefore, these two planets are very appealing objects for studying the effects of mass loss. Interestingly, two other objects with similar bulk densities, namely K2-18 b (M_p = 8.63M_⊕, R_p = 2.6R_⊕, ρ_p = 2.67 g cm⁻³; see Benneke et al. (2019)) and GJ 3470b (M_p = 13.9M_⊕, R_p = 4.6R_⊕, ρ_p = 0.80 g cm⁻³; see Awiphan et al. (2016)), have shown the presence of atmospheric escape through the detection of Lyα absorption during transit (Bourrier et al. 2018; dos Santos et al. 2020).

Here, we use the 1D hydrodynamic atmospheric escape model presented in Allan & Vidotto (2019) to predict the behavior of the planetary atmospheres under the influence of the photoionizing flux of the host star. With this, we can infer the current properties of the planetary upper atmospheres, including the mass-loss rate. For the output of this 1D model, we study the atmospheric signatures of TOI-421 b and TOI-421 c for the neutral hydrogen Lyα and H-α lines during transit. Our model uses as input the XUV stellar luminosity, which was derived considering the median $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$ value (−4.93 ± 0.04) from Table 3, first by converting it into a Ca ii H and K chromospheric emission flux using the equations listed in Fossati et al. (2017b) and then by converting this emission in XUV flux using the scaling relations of Linsky et al. (2013) and Linsky et al. (2014). We find an XUV flux at a distance of 1 au of F_XUV = 23.12 erg cm⁻² s⁻¹, corresponding to an XUV luminosity of L_XUV = 6.5 × 10²⁸ erg s⁻¹. This implies XUV fluxes of F_XUV = 1654.8 erg cm² s⁻¹ and F_XUV = 7452.0 erg cm² s⁻¹ at the distances of planet b and c, respectively.

Figure 11 summarizes some key properties of the escaping atmosphere for both planets: radial component of the velocity (top), temperature (middle), and ionization fraction (bottom). For planet b (inner planet), we derive an atmospheric escape rate of 4.5 × 10¹⁰ g s⁻¹. We calculate the Roche lobe distance to be at 9.7 R_p; at this distance, the speed of the escaping material is about 30 km s⁻¹. Planet c has a similar escape rate of 4.4 × 10¹⁰ g s⁻¹, with material reaching a speed of 28 km s⁻¹ at the Roche lobe distance (14.3 R_p). The reason for the comparable escape rate is that, although planet c receives an XUV flux that is 4.5 times smaller than planet b (due to its larger orbital distance), it has a lower surface gravity (g_c ≃ 24% of Jupiter's gravity versus g_b ≃ 38% for planet b). It is more difficult for low-gravity planets to hold onto their atmospheres, thus the lower gravity of planet c compensates for its lower incident XUV flux, reaching comparable escape to planet b.

Zoom In Zoom Out Reset image size

Despite the similarities in the escape rates and velocities, we predict different Lyα transit absorptions for these two planets. Figure 12 shows the predicted light curves at Lyα line center, where we see that planet b (the inner planet) shows a maximum absorption of 35% and planet c (the outer planet) shows a maximum absorption of 53%. These different absorptions are caused by the different ionization fractions in each planet's atmosphere (see bottom panel of Figure 11), with planet c showing more neutral hydrogen in its atmosphere. The light curves presented in Figure 12 are symmetric about mid-transit. This is due to the one-dimensional geometry of the model, hence the assumption of spherically symmetric planetary atmospheres. However, we expect light curves to be asymmetric with respect to mid-transit, but these asymmetries can only be captured by 3D models that include interactions with the stellar wind (e.g., Villarreal D'Angelo et al. 2014, 2018), which we leave to a forthcoming paper.

Zoom In Zoom Out Reset image size

To encourage future observations in the UV, we made use of the absorption profile obtained from the upper atmosphere simulations and computed the expected Lyα profile as observed with the G140M grating of the STIS spectrograph on board the Hubble Space Telescope (HST). Figure 13 shows the resulting profiles at three different times, reproducing an out-of-transit observation (black line) and an observation at mid-transit of planet b (blue-dashed line) and planet c (red-dashed line). To compute the out-of-transit Lyα profile of TOI-421, we scaled the intrinsic Lyα profile of ξ Boo A (G8 V, d = 6.70 pc, R_⋆ = 0.78 R_⊙) derived by Wood et al. (2005). We assumed that the strength of the profiles of the two stars should be similar, as the stellar type and activity index for ξ Boo A ($\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$ ∼ −4.40; see Morgenthaler et al. (2012)) are close to those of TOI-421. Most of the predicted absorption will be hidden by the ISM absorption and geocoronal emission as shown in Figure 13. However, the ISM will affect the spectrum—absorbing a significant fraction of its flux beyond almost 100 pc, and probably most of it by 200 pc—while TOI-421 is 75 pc distant. Moreover, because of the interaction of the escaping planetary upper atmosphere with the stellar wind and the resolution of the G140M grating, significant planetary Lyα absorption will be visible in the line wings, making it detectable. The attenuation caused by the ISM over the intrinsic profile of TOI-421 is assumed to be the same as the one estimated for ξ Boo A, also computed in the work of Wood et al. (2005). This is a crude approximation, but a more adequate estimation would require at least a near-ultraviolet observation of the ISM absorption at the position of the Mgii h and k lines (Wood et al. 2005). Using HST, we expect to be able to detect the absorption signal caused by planet c during transit, as this is the largest absorption modeled, in accordance with what is shown in Figure 12.

Zoom In Zoom Out Reset image size

We also compute the H-α transit light curves following the approach outlined in Allan & Vidotto (2019). Although we predict a large absorption in Lyα, no appreciable absorption is expected to be detectable in H-α. This is because most of the hydrogen in the atmospheres of TOI-421 b and TOI-421 c is in the ground state—for H-α to be formed, it is required that some hydrogen be in the first excited state. The fact that we do not predict any H-α absorption in our model needs to be reassessed in 3D calculations. Recently, Villarreal D'Angelo et al. (2020) demonstrated that the interaction between the upper planetary atmosphere and the stellar wind, which is only captured in 3D models, could increase the atmospheric temperature in the interaction zone. As a result, this can increase the number of neutral hydrogen in the first excited state, possibly enhancing H-α absorption. We will further explore this in a future work.

9.2. Planetary Atmospheric Evolution

In addition to the hydrodynamic model presented in Section 9.1, as a cross-check, we also computed the mass-loss rates employing the interpolation routine of Kubyshkina et al. (2018), which is based on 1D hydrodynamic simulations, obtaining a value of 2.7 × 10¹⁰ g s⁻¹ for planet b (inner planet) and 1.6 × 10¹⁰ g s⁻¹ for planet c (outer planet), in agreement with the results obtained using the hydrodynamic atmospheric escape model of Allan & Vidotto (2019). We further notice the low Λ⁸⁹ value (Λ = 19.33) of planet b, which implies that the planetary gravity is hardly capable of holding a hydrogen-dominated atmosphere (Fossati et al. 2017a). In fact, this is remarkable because other planets with similarly low Λ values have an average density indicative of a mostly rocky composition (Gandolfi et al. 2017). This is confirmed by the mass-loss rates of planets b and c, which imply that they would have lost 20% and 5%, respectively, of their mass over 1 Gyr. While planet c should have an atmospheric mass fraction large enough to sustain such an intense mass loss over Gyr, the escape is probably too intense for planet b to be able to still retain a hydrogen-dominated atmosphere, as instead suggested by the average density.

We employed the tool presented by Kubyshkina et al. (2019a, 2019b) to constrain the atmospheric evolution of both planets, their initial atmospheric mass fractions, and the evolution of the rotation rate (and therefore also of the XUV emission) of the host star. In short, the framework mixes three ingredients: a model of the stellar XUV flux evolution (after Pizzolato et al. 2003; Mamajek & Hillenbrand 2008; Sanz-Forcada et al. 2011; Wright et al. 2011), a model relating planetary parameters and atmospheric mass (after Johnstone et al. 2015b; Stökl et al. 2015), and a model computing escape (after Kubyshkina et al. 2018). The framework also accounts for the evolution of the stellar bolometric luminosity, and hence planetary equilibrium temperature, using the MESA Isochrones and Stellar Tracks (MIST; Paxton et al. 2018) grid.

For a given core mass, the framework sets the core radius assuming an Earth-like density, and the atmospheric mass is considered to be negligible (Owen & Wu 2017). Then, starting at 5 Myr (the assumed age of the dispersal of the protoplanetary disk), at each time step the framework extracts the mass-loss rate from the grid based on the stellar flux and system parameters, using it to update the atmospheric mass fraction and the planetary radius. This procedure is then repeated until the age of the system is reached or the planetary atmosphere is completely escaped. The framework simulates the atmospheric evolution of both planets, simultaneously. The main framework's assumption is that the analyzed planets have (or had) a hydrogen-dominated atmosphere and that the planetary orbital separation does not change after the dispersal of the protoplanetary disk.

The input parameters of the framework are planetary masses, planetary radii, orbital separations, current stellar rotation rate, and stellar mass, while the free parameters are the index of the power law describing the evolution of the stellar XUV flux and the initial planetary radius (i.e., the initial atmospheric mass fraction at the time of the dispersal of the protoplanetary disk; f_at). The input parameters are set equal to the measurements with Gaussian priors having a width equal to the measurement uncertainties, while we take flat priors for the output parameters. The output parameters are constrained by implementing the atmospheric evolution algorithm in a Bayesian framework employing the MCMC tool of Cubillos et al. (2017b).

Figure 14 shows the main results of the simulation. The posteriors on the input parameters of the host star and of the outer planet (TOI-421 c) match the priors. The evolution of the stellar rotation rate (or XUV emission) is mostly unconstrained. For the outer planet, the analysis leads to a rather tight constraint on the initial atmospheric mass fraction of about 30%, and this result also holds when running the analysis solely on the outer planet, meaning that the anomaly found for the inner planet (see below) does not affect the other results. The outer planet could not have accreted an atmosphere much larger than f_at ≈ 30% while in the disk, because otherwise the stellar XUV emission would have not been able to remove enough of it to obtain the currently observed radius, even if the star was a fast rotator. Similarly, the planet could not have accreted an atmosphere much smaller than f_at ≈ 30% while in the disk, because otherwise the stellar XUV emission would have removed too much atmosphere, given the observed radius, even if the star was a slow rotator.

Zoom In Zoom Out Reset image size

The result obtained for the inner planet is extremely interesting. The framework is unable to find a configuration in which the planet is capable of retaining enough atmosphere to match the measured planetary mass, radius, and orbital separation. This is why the posterior of the planetary mass is slightly shifted toward higher masses compared to the prior (first top left panel of Figure 14). Moreover, the posterior of the planetary radius is significantly shifted toward smaller radii compared to the prior (third top panel of Figure 14). In other words, given the system parameters, the framework finds that the inner planet always loses its hydrogen atmosphere, regardless of the evolution of the stellar XUV emission. An inspection of the atmospheric evolutionary tracks indicates that the inner planet is expected to completely lose its atmosphere within 1 Gyr, while the estimated age of the system is significantly larger. We reran the simulation, looking for the planetary parameters that would enable the posteriors on mass and radius not to vary from the priors, obtaining either an orbital separation of about twice the measured one (keeping mass and radius equal to the measured values), a planetary mass of about 16 M_⊕ (keeping radius and orbital separation equal to the measured values), or a planetary radius of about 2 R_⊕ (keeping mass and orbital separation equal to the measured values). In the last two options, the planet would not host a hydrogen-dominated atmosphere.

CoRoT-24b was the first planet identified to have a low bulk density, compatible with the presence of a hydrogen-dominated atmosphere, but at the same time to also be subject to mass loss too extreme for hosting one (Lammer et al. 2016). Cubillos et al. (2017a) analyzed the upper atmospheric properties and high-energy irradiation of a large sample of mini-Neptunes, detected mostly by the Kepler satellite, finding that 15% of them share this same peculiar property. There is a range of possible solutions to this puzzle. One of the main assumptions in the atmospheric evolution framework is that orbital separations do not change with time following the dispersal of the protoplanetary nebula—which may not be the case, for example, if the system had a close enough encounter with another star in the past. It may also be that the hydrodynamic model overestimates the mass-loss rates, although past comparisons between observations and hydrodynamic models would tend to exclude order-of-magnitude errors in the computed rates. One further possibility is a bias in the measured planetary parameters (mass and/or radius), maybe caused by the presence of other undetected planets in the system biasing the planetary mass measurement, but this seems unlikely, given the quantity and quality of data. The mainstream explanation for this kind of planet is the presence of high-altitude aerosols that would lead overestimation of planetary radius (Lammer et al. 2016; Cubillos et al. 2017a; Gao & Zhang 2020). Future atmospheric characterization observations, particularly those at low resolution that are more sensitive to aerosols, will be able to identify whether this is the case or not (e.g., Libby-Roberts et al. 2020). There is also the additional possibility that the crust released a significant amount of light gases in the atmosphere, counteracting the effect of escape (e.g., Kite et al. 2019).

9.3. Simulated HST WFC3 Retrievals

With its large radius (∼5 R_⊕), TOI-421 c represents an excellent target for atmospheric characterization. TOI-421 b is somewhat more challenging; its scale height is comparable to TOI-421 c, but because c is smaller by a factor of two, the transit signal corresponding to one scale height will be as well, similarly to Lyα absorption. To assess how well the atmospheric properties of TOI-421 b and TOI-421 c could be derived from observations, we modeled the planetary atmospheres and derived transmission spectra with the open-source petitRADTRANS package (Mollière et al. 2019). The atmospheres were set up to be isothermal, at the equilibrium temperature of the planets. The absorber abundances were obtained from assuming chemical equilibrium, calculated with the chemistry module that is part of petitCODE (Mollière et al. 2017). We assumed a solar C/O ratio of 0.55, and two different metallicity values: 3 and 100× solar (Jupiter and Neptune-like, respectively). We also introduced a gray cloud deck, the position of which was varied between 100 and 10⁻⁵ bar, in 1 dex steps. We considered 100 bar to be the cloud-free model, as the atmosphere will become optically thick at lower pressures. The following gas opacities were included: the line opacities of H₂O, CH₄, CO, CO₂, Na, K; the Rayleigh scattering cross sections of H₂O, CH₄, CO, CO₂, H₂, and He; and H2–H2 and H2–He collision-induced absorption.

The atmospheric models described above were then retrieved with petitRADTRANS, using the PyMultiNest package (Buchner et al. 2014), which makes use of the nested sampling implementation MultiNest developed by Feroz et al. (2009). We created synthetic HST observations with the WFC3 instrument, assuming 12 wavelength points spaced equidistantly from 1.12 to 1.65 μm. The error on the flux decrease during transit was assumed to be 35 ppm for TOI-421 b and 33 ppm for TOI-421 c per channel, which we calculated for a single transit of the two planets, using the Pandexo_HST tool.⁹⁰ For reference, this is about 1/3 of the signal of one atmospheric scale height of TOI-421 c, when assuming a solar composition. For every input model, we wish to characterize whether or not the molecular features of H₂O or CH₄ can be identified in the spectra. To this end, we follow the technique outlined in Benneke & Seager (2013), i.e., we first retrieve the atmospheric temperature, reference pressure, cloud deck pressure, and vertically constant absorber abundances freely. This is called setup (i) in the following. We then remove the abundance parameter and opacity of either H₂O or CH₄ from the retrieval; these are setups (ii) and (iii). The Bayes factor B between model (i) and (ii) will constrain how confidently the atmospheric features of H₂O can be detected, while the B between model (i) and (iii) informs us about how reliably CH₄ can be detected.

It has recently been found that observational uncertainties in the range of 30 ppm can lead to significant differences in retrieved atmospheric abundances and temperatures when comparing the results of various retrieval codes (Barstow et al. 2020). These discrepancies most likely arise due to differences in the opacities, either because of the use of different line lists, or due to the choices made when converting line lists into opacities, such as the line broadening or cutoff. Because we use the same model to make the synthetic observations and the retrievals, our results can therefore be regarded as a limiting case, where the abovementioned issues are negligible. Moreover, because our discussion here focuses on the detection of molecular features, instead of constraining their abundances, we deem our approach acceptable for the exploratory study presented here. When running retrievals for real observations of similar data quality in order to constrain abundances and other atmospheric parameters, the use of multiple retrieval codes or opacity treatments (varying line lists, the broadening description, etc.) is recommended.

9.3.1. TOI-421 b

For the case of three times solar metallicity, we found that for cloud pressures larger than 10 mbar, substantial atmospheric features can be retrieved (we find B > 3; see Kass & Raftery (1995) for a definition of the B thresholds). For a metallicity of 100 times solar, this transition likewise occurs for pressures larger than 10 mbar. Due to the high equilibrium temperature of the planet, no CH₄ can be detected in the synthetic observations. This is because, at high temperatures, CO is chemically favored as the main C-bearing species instead of CH₄.

As an example, the left panel of Figure 15 shows the retrievals of the 100 times solar enrichment, clear atmosphere synthetic observation with the full model, and the models without CH₄ and H₂O. Because of the random noise instantiation, H₂O is only weakly detected in this example (B = 2.7), while running this test multiple times places the average B in favor of including H₂O at a value of B = 6.

Zoom In Zoom Out Reset image size

9.3.2. TOI-421 c

For the three times solar metallicity case, we found that for cloud pressures larger than 1 mbar, substantial atmospheric features can be retrieved. For a metallicity of 100 times solar, this transition occurs for pressures larger than 0.1 bar. Due to the smaller temperature and the assumption of equilibrium chemistry in the atmospheric model used to generate the observations, we find that the signal of CH₄ can be more confidently detected in the atmosphere of the planet than that of H₂O.

The right panel of Figure 15 shows the retrievals of the 100 times solar enrichment, clear atmosphere synthetic observation for the full retrieval model and the model without CH₄ or H₂O. The average B values for detecting CH₄ or H₂O are B = 4 and B = 3, respectively.

9.4. Prospects for Detecting Transit-timing Variation and the Rossiter–McLaughlin Effect

We have presented the discovery of a Neptune-sized planet and a sub-Neptune transiting TOI-421 (BD-14 1137, TIC 94986319), a G9 dwarf star observed by TESS. The host star is the primary component of a visual binary. Our RV follow-up observations led to the confirmation of the outer Neptune-sized planet (TOI-421 b) and the discovery of the second inner sub-Neptune (TOI-421 c), that we also found to transit its host star. We have determined both stellar and planetary parameters. We found that TOI-421 is a relatively quiet star with an activity index of $\mathrm{log}\,{R}_{\mathrm{HK}}^{{\prime} }$ = −4.93 ± 0.04. Based on our analysis of the HARPS and WASP-South data, we found that the intrinsic activity of TOI-421 can be explained mainly by plages.

Our TTV analysis shows that TOI-421 is an ideal target to compare planetary masses determined via TTV and Doppler techniques. We aim for future additional transit observations to explore this in more detail.

TOI-421 b and TOI-421 c are very appealing and suitable targets for atmospheric characterization. They are both expected to host extended atmospheres, showing significant signal in the Lyα line. Moreover, the atmospheric retrievals demonstrated that we can detect CH₄ in the atmosphere of the outer planet (TOI-421 c) if the atmosphere is in chemical equilibrium, and atmospheric evolution simulations showed that the inner planet (TOI-421 b) appears to be among the small sample of peculiar super puffy mini-Neptunes, making it also more intriguing for atmospheric studies and evolution theories. This multiplanet system, with its astonishing characteristics, would be a prime target for the upcoming JWST observations. Indeed, the two planets are among the first 30 targets with the highest expected S/Ns, as shown in Figure 16. Using the sample of exoplanets with R < 6 R_⊕, totaling more than 2000 exoplanets,⁹¹ TOI-421 b and TOI-421 c are within the top 30 most favorable targets for atmospheric characterization. This atmospheric characterization metric is based on a J-band, JWST-style observation, and is detailed in Niraula et al. (2017). It is particular noteworthy that this metric is scaled by the frequency of transits. This is motivated by the expectation that sensitive atmospheric observations will likely require many transits to build sufficient signal, and it may be prohibitive to accumulate the needed transits for longer-period exoplanets. Therefore, we used a metric that optimizes the S/N over a period of time rather than a per-transit metric. Regardless of the nuances of the metric, the TOI-421's planets are highly attractive targets for characterization of both their bound and extended atmospheres.

Zoom In Zoom Out Reset image size

This work was supported by the KESPRINT collaboration, an international consortium devoted to the characterization and research of exoplanets discovered with space-based missions.

We acknowledge the use of public TESS Alert data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center.

Funding for the TESS mission is provided by NASA's Science Mission directorate. We acknowledge the use of public TESS Alert data from pipelines at the TESS Science Office and at the TESS Science Processing Operations Center. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center for the production of the SPOC data products.

Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (NASA).

We are very grateful to the NOT and ESO staff members for their unique and superb support during the observations. A.A.V. and D.K. acknowledge funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement No 817540, ASTROFLOW). P.M. thanks Thomas Mikal-Evans and Laura Kreidberg for help computing predicted HST uncertainties for TOI-421. P.M. acknowledges support from the European Research Council under the European Union's Horizon 2020 research and innovation program under grant agreement No. 832428.

L.M.S. and D.G. gratefully acknowledge financial support from the CRT foundation under grant No. 2018.2323 "Gaseous or rocky? Unveiling the nature of small worlds." R.B. acknowledges support from FONDECYT Post-doctoral Fellowship Project 3180246. A.J. and R.B. acknowledge support from project IC120009 "Millennium Institute of Astrophysics (MAS)" of the Millennium Science Initiative, Chilean Ministry of Economy. A.J. acknowledges additional support from FONDECYT project 1171208. C.V.D. acknowledges the funding from the Irish Research Council through the postdoctoral fellowship (Project ID: GOIPD/2018/659). C.M.P., M.F., and I.G. gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19 and 136/13). L.S. acknowledges financial support from the Australian Research Council (Discovery Project 170100521). D.D. acknowledges support from NASA through Caltech/JPL grant RSA-1006130 and through the TESS Guest Investigator Program grant 80NSSC19K1727. D.H. acknowledges support from the Alfred P. Sloan Foundation, the National Aeronautics and Space Administration (80NSSC18K1585, 80NSSC19K0379), and the National Science Foundation (AST-1717000). S.M. acknowledges support from the Spanish Ministry with the Ramon y Cajal fellowship number RYC-2015-17697. L.G.C. thanks the support from grant FPI-SO from the Spanish Ministry of Economy and Competitiveness (MINECO), under research project SEV-2015-0548-17-2 and predoctoral contract BES-2017-082610. A.R.G.S. acknowledges support from NASA under grant No. NNX17AF27G. J.N.W. thanks the Heising-Simons foundation for support. R.A.G. acknowledges support from PLATO and GOLF CNES grants. M.R.D. acknowledges support from CONICYT-PFCHA/Doctorado Nacional-21140646, Chile. K.W.F.L., J.K., S.z.C.s., M.E., S.G., A.P.H., M.P., and H.R. acknowledge support by DFG grants PA525/18-1, PA525/19-1, PA525/20-1, HA3279/12-1, and RA714/14-1 within the DFG Schwerpunkt SPP 1992, "Exploring the Diversity of Extrasolar Planets" 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.