TKS X: Confirmation of TOI-1444b and a Comparative Analysis of the Ultra-short-period Planets with Hot Neptunes

The TESS team reported one object of interest, TOI-1444.01, with an orbital period of 0.47 day. Beginning with the transit parameters reported on the ExoFOP website, ³⁶ we removed the data points taken within the transit window of TOI-1444.01. Using the resultant, out-of-transit light curve, we tried to measure the stellar rotation period of TOI-1444 using both a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982) and an auto-correlation function (ACF, McQuillan et al. 2014). These two methods basically look for any quasi-periodic flux modulations that may be attributed to stellar rotation coupled with surface magnetic activity. Neither the periodogram nor the ACF detected a signal above a 1% false positive rate. Moreover, the highest peaks reported by these two methods did not agree. Therefore, we could not measure the rotation period of the host from the flux variation seen in the TESS light curve. The weak flux variation (standard deviation of about 1200 ppm for a 2 minute cadence over a 300 day baseline) may be the result of a slow rotation and/or low stellar activity. Rotational broadening of TOI-1444 was not detectable in our spectrum above other broadening effects. Following Petigura (2015), we set an upper limit of v sini <2 km s⁻¹. The proxy for chromospheric activity in the Ca II H&K lines S_HK is 0.145 ± 0.004. This is slightly lower than the median S_HK of stars of similar B–V color (0.146, Isaacson & Fischer 2010).

TOI-1444b was initially detected by the TESS Science Processing Operations Center (SPOC; Jenkins et al. 2016) in a transit search of Sectors 15 and 16 that occurred. The 0.47 day signal was detected at a 7.6σ level with an adaptive, noise-compensating, matched filter (Jenkins et al. 2010); it passed all the diagnostic tests performed and published in the resulting data validation reports. The vetting included tests for eclipsing binaries, such as an odd and even depth variation test, a secondary eclipse test, and a ghost diagnostic test. The difference imaging centroid test showed that the source of the transit signature was consistent with the target star TIC 258514800, but could not exclude nearby stars in the TIC catalog. The signal strength grew as additional observations were collected by TESS and the difference imaging centroid test located the source within 36 once the full set of observations were completed in Sector 26.

We searched the TESS light curve for any additional transit signals, particularly around the 16 day periodicity of the candidate planetary signal seen in the RV data set (Section 5). We first removed the data points within the transit window of TOI-1444.01. We then fitted a cubic spline of length 1.5 days to detrend any long-term stellar or instrumental flux variation. We applied the box-least-squares algorithm (BLS, Kovács et al. 2002) to the resultant light curve. Our BLS pipeline is implemented in C++ and has yielded a number of planet discoveries including other ultra-short-period planets, e.g., K2-131b (Dai et al. 2017) and K2-141b (Barragán et al. 2018). We followed the suggested improvement of BLS as outlined in Ofir (2014). This involves using a nonlinear frequency grid given the theoretical scaling of transit duration with an orbital period for stars of a certain mean density. We also adopted the signal detection efficiency (SDE) defined in Ofir (2014) to quantify the significance of a BLS signal. In short, SDE is the local variation of the BLS spectrum normalized by the local standard deviation. It helps to remove any period-related systematics.

We recovered TOI-1444.01 with an SDE = 21. However, we did not detect any additional transiting signal with an SDE larger than 5. Visual inspection of the phase-folded light curve shows that none of the top candidate signals has a transit-like shape. We also did not detect any convincing single-transit events visually. Notably, there is no transiting signal consistent with the 16 day orbital period of the candidate planetary signal in the RV data set. No additional transiting signal was found by the SPOC pipeline either.

We modeled the transit light curve of TOI-1444.01 or TOI-1444b to constrain its transit parameters. We started from the transit ephemeris reported by the TESS team. We used the Python package Batman (Kreidberg 2015) to generate the model transit light curves. We also imposed a prior on the host-star mean density using the result in Table 1. This precise prior on mean stellar density helps to break some of the degeneracy in modeling transit morphology and often leads to improved transit parameters (Seager & Mallén-Ornelas 2003). We adopted a quadratic limb-darkening profile. We imposed Gaussian priors (width of 0.3) on the limb-darkening coefficients u₁ and u₂. We queried EXOFAST ³⁷ (Eastman et al. 2013) for the theoretical limb darkening given the spectroscopic parameters of TOI-1444 (u₁ = 0.48, u₂ = 0.22). We also adopted the efficient sampling reparameterization of limb-darkening coefficients proposed by Kipping (2013). The other parameters in our transit model are the orbital period P_orb, the time of conjunction T_c, the planet-to-star radius ratio R_p/R_⋆, the scaled orbital distance a/R_⋆ and the impact parameter $b\equiv a\cos i/{R}_{\star }$.

We first fit all transits of TOI-1444b globally. We found the best-fit solution using the Levenberg–Marquardt method, specifically using that implemented in Python package lmfit. Using this global fit as a template, we then fit each individual transit, allowing the mid-transit time to vary freely. The resultant transit times do not show quasi-sinusoidal variation as one would expect in the case of transit-timing variations. We did not detect any prominent periodicity; instead, it is consistent with a linear ephemeris, which we assume in subsequent analysis.

We sampled the posterior distribution of transit parameters with the affine-invariant Markov Chain Monte Carlo method as implemented in Python package emcee (Foreman-Mackey et al. 2013). We used 128 walkers, starting from the maximum likelihood solution found by lmfit. With 50,000 links, the Gelman–Rubin potential scale reduction factor reduced to below 1.01, suggesting convergence of the sampling process. We summarize the resultant posterior distribution in Table 2. Figure 2 shows the phase-folded and binned TESS light curve and the best-fit transit model.

Zoom In Zoom Out Reset image size

Table 2. Planetary Parameters of TOI-1444

Parameter	Symbol	Posterior Distribution
Planet b
Planet/Star Radius Ratio	Rp / R⋆	${0.01410}_{-0.00058}^{+0.00060}$
Time of Conjunction (BJD-2457000)	tc	1711.3676 ± 0.0016
Impact Parameter	b	0.38 ± 0.14
Scaled Semimajor Axis	a/R⋆	${2.73}_{-0.074}^{+0.078}$
Orbital Inclination (deg)	i	${82}_{-2}^{+3}$
Orbital Eccentricity	e	0 (fixed)
Orbital Period (days)	Porb	0.4702694 ± 0.0000044
Semiamplitude (m/s)	K	${3.30}_{-0.59}^{+0.58}$
Planetary Radius (R⊕)	Rp	1.397±0.064
Planetary Mass (M⊕)	Mp	3.87±0.71
Secondary Eclipse Depth (ppm)	δsec	27 ± 12
Time of Secondary Eclipse (days)	tsec	0.236 ± 0.019
Amplitude of Illumination Effect (ppm)	Aill	24 ± 10
Phase Offset of Illumination Effect (°)	θill	16 ± 27
Planet c
Time of Conjunction (BJD-2457000)	tc	713.0 ± 2.0
Orbital Period (days)	Porb	16.066 ± 0.024
Orbital Eccentricity	e	0 (fixed)
Semiamplitude (m s−1)	K	${3.12}_{-0.76}^{+0.75}$
Projected Planetary Mass (M⊕)	${M}_{{\rm{p}}}\sin i$	11.8 ± 2.9
RV Jitter (m s−1)	σ	3.2 ± 0.4

Download table as:  ASCIITypeset image

We looked for any phase-curve variation and secondary eclipse of TOI-1444b in the TESS band. Again, we removed data taken within the transit window of TOI-1444b. We then detrended any long-term stellar variation or systematic effect using the procedure outlined in Sanchis-Ojeda et al. (2013a). In short, a linear function in time whose width is at least one orbital period is fitted to remove local flux variation longer than the orbital period. In practice, we tried 1×, 2×, or 3 × P_orb and found nearly identical results. We also compared PDC_SAP and SAP versions of TESS photometry and found no substantial difference. Figure 3 shows the phase-curve variation from all TESS sectors of PDC_SAP photometry.

Zoom In Zoom Out Reset image size

We then modeled the phase-curve variation and secondary eclipse simultaneously. We modeled the secondary eclipse using Batman. We fixed the transit parameters using the best-fit primary transit solution and changed the limb-darkening coefficients to 0. The phase curve variation was modeled as a Lambertian disk (e.g., Demory et al. 2016). The parameters in this model include the depth of secondary eclipse δ_sec, the amplitude of illumination effect A_ill, the time of secondary eclipse t_sec, and the phase offset of the illumination effect θ_ill. To constrain the posterior distribution, we similarly ran an MCMC analysis with emcee (Foreman-Mackey et al. 2013). We found that both δ_sec = 27 ± 12 ppm and A_ill = 24 ± 10 ppm showed marginal detection of no more than 2.5σ. The fact that δ_sec are similar to A_ill in amplitude even though they are fitted separately boosts our confidence in these detections. t_sec centers around the half-way from mid-transit although with substantial error bar t_sec = 0.236 ± 0.019 day or 0.502 ± 0.040 in the orbital phase. Again, this is expected given the very short tidal circularization timescale of a planet on such a short-period orbit. The phase curve offset θ_ill = 16 ± 27° shows a very weak preference for an eastward offset.

The TESS passband (600–1000 nm) extends substantially into the infrared so that one may expect to see thermal emission from the planet as well as reflected starlight. However, with only one band, we could not break the degeneracy between reflected stellar light and thermal emission from the planet. Figure 4 captures this degeneracy in the Bond albedo versus TESS-band geometric albedo plane. The blue shaded region shows the 68% confidence interval. Given how wide the confidence interval is, we refrain from making a strong interpretation of the marginal 2.5σ phase-curve detection in the TESS band. Phase-curve observation in the infrared would consolidate the phase-curve variation and may reveal the surface properties of TOI-1444b.

Zoom In Zoom Out Reset image size

We observed TOI-1444 on 2020 May 02 during the transit window of planet b as predicted by the SPOC pipeline analysis of TESS Sectors 15 and 16. The observation was carried out in the Pan-STARSS z-short band from the Las Cumbres Observatory Global Telescope (LCOGT; Brown et al. 2013) 1.0 m network node at the McDonald Observatory. The 4096 × 4096 LCO SINISTRO cameras have an image scale of 0389 per pixel, providing a field of view of 26' × 26'. The standard LCOGT BANZAI pipeline (McCully et al. 2018) was used to calibrate the images and the photometric data were extracted with the AstroImage (AIJ) software package (Collins et al. 2017).

Since the ∼200 ppm event detected by the SPOC pipeline is generally too shallow to detect with ground-based observations, we checked for a faint nearby eclipsing binary (NEB) that could be contaminating the SPOC photometric aperture. To account for possible contamination from the wings of neighboring star PSFs, we searched for NEBs at the positions of Gaia DR2 stars out to 25 from the target star. If fully blended in the SPOC aperture, a neighboring star that is fainter than the target star by 9.3 magnitudes in the TESS band could produce the SPOC-reported flux deficit at mid-transit (assuming a 100% eclipse). To account for possible delta-magnitude differences between the TESS band and Pan-STARSS z-short band, we included an extra 0.5 magnitude fainter (down to TESS-band magnitude 20.0). We visually compared the light curves of the 85 nearby stars that meet our search criteria with models that indicate the timing and depth needed to produce the ∼200 ppm event in the SPOC photometric aperture. We found no evidence of an NEB that might be responsible for the SPOC detection. By a process of elimination, we conclude that the transit is likely occurring on TOI-1444.

We first applied the Python package RVSearch ³⁸ to determine the number of planetary signals in our RV data set. In short, RVSearch sequentially looks for peaks in the LS periodogram of the RV data set after removing the best-fit Keplerian model of planetary signals detected in previous iterations. The code stops when the Bayesian Information Criterion (Δ BIC) no longer favors the addition of another planetary signal. RVSearch has been widely tested on known planetary systems. For more detail on RVSearch, we refer interested readers to Rosenthal et al. (2021).

Applying RVSearch to TOI-1444, a strong peak at the transiting period of TOI-1444b is recovered (Figure 5). The next strongest peak is a 38 day signal that is very close to a 36.5 day alias of the annual variation and was disregarded as a possible planetary signal. RVSearch detects a candidate 16 day signal. No corresponding peak is seen in the RV window function or the activity index S_HK. Moreover, the strength of the 16 day periodicity steadily increased as we included more and more RV data in the periodogram. In contrast, a signal caused by stellar activity loses coherence over time because the underlying surface magnetic activity typically emerges and subsides on a weeks to months timescale for solar-like stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

RVSearch can also perform injection and recovery tests of planetary signals in the M_p sini-P_orb plane. The result is a completeness contour showing the sensitivity of the RV data set at hand to planetary signals of different strength and periodicity in addtion to the detected planets (e.g., Figure 7). Given the current HIRES RV data set for TOI-1444, we were unable to identify a third planetary signal. The completeness contour for that third planet is shown in Figure 7. With the current HIRES data set, other sub-Neptune planets (<10M_⊕) within 1 au of TOI-1444 can easily remain hidden from detection.

Zoom In Zoom Out Reset image size

We first modeled the RV data set using Keplerian models. We used the publicly available Python package RadVel (Fulton et al. 2018). Each planetary signal is described by its orbital period P_orb, time of inferior conjunction T_c , eccentricity e, argument of periastron ω, and the RV semiamplitude K. We allowed an RV offset γ, a linear RV trend $\dot{\gamma }$, and we included a jitter term to encapsulate any additional astrophysical or instrumental noise. We reparameterized e and ω into $\sqrt{e}$ cosω, $\sqrt{e}$ sinω. Since the ephemeris is much better constrained with the transit data, we imposed Gaussian priors on P_orb and T_c using the results from Section 4.3. We imposed uniform priors on the RV semiamplitude K, the jitter σ_jit, $\sqrt{e}$ cosω (with range [−1, 1]), $\sqrt{e}$ sinω ([−1, 1]), γ and $\dot{\gamma }$. Our likelihood function is as follows:

$$\begin{eqnarray}\begin{array}{rcl}{ \mathcal L } & = & \prod _{i}\left(\displaystyle \frac{1}{\sqrt{2\pi ({\sigma }_{i}^{2}+{\sigma }_{\mathrm{jit}}^{2})}}\right.\\ & & \times \exp \left.\left[-\displaystyle \frac{{[\mathrm{RV}({t}_{i})-{ \mathcal M }({t}_{i})]}^{2}}{2({\sigma }_{i}^{2}+{\sigma }_{\mathrm{jit}}^{2})}\right]\right)\end{array}\end{eqnarray} \tag{ 1 }$$

where RV(t_i ) is measured RV, ${ \mathcal M }({t}_{i})$ is the sum of Keplerian planetary signals, a constant RV offset γ, and a linear RV trend $\dot{\gamma }{t}_{i};$ and σ_i is the internal uncertainty.

We varied the number of planets in our model and we tested if the current data set supports nonzero eccentricity and $\dot{\gamma }$. We selected the best model using both the Bayesian information criterion (Δ BIC) and Akaike information criterion (Δ AIC). Our best-fit model (lowest BIC and AIC) contains the RV signals of planet b and a nontransiting planet with a 16 day orbit identified by RVSearch. Nonzero eccentricity for either planet was not preferred by the current data set. A linear RV trend of $\dot{\gamma }=-0.008\pm 0.012$ m s⁻¹ is marginally disfavored by the current data set with ΔBIC ≈ − 2. We then performed an MCMC analysis to sample a posterior distribution of the various parameters. The sampling procedure is similar to that described in Section 4.3. We summarize the posterior distribution in Table 2.

The RV residuals after removing the best two-planet Keplerian model still show a root-mean-square variation (rms) of 3.4 m s⁻¹ which is substantially larger than the internal uncertainties in Table 3. Moreover, the RV residuals visually display some correlated noise component in Figure 8. We investigated if these residual noises can be tidied up by a Gaussian process (GP) model (e.g., Haywood et al. 2014; Grunblatt et al. 2015) and disentangle the planetary signal from the stellar activity.

Zoom In Zoom Out Reset image size

We adopted the GP model used in our previous works (Dai et al. 2017, 2019). Stellar surface magnetic activity rotating in and out of view of the observer is chiefly responsible for the correlated stellar noise in RV measurements, the quasi-sinusoidal flux variation in the light curve, and the variation of the chromospheric activity index S_HK (e.g., Aigrain et al. 2012). In other words, these effects stem from the same underlying physical process. We model them with a common GP model. The plan of attack is to train a quasi-periodic GP model on the out-of-transit TESS light curves and the HIRES S_HK index before applying it to the RV data set. Our kernel has the following form:

$$\begin{eqnarray}\begin{array}{rcl}{C}_{i,j} & = & {h}^{2}\exp \left[-\displaystyle \frac{{\left({t}_{i}-{t}_{j}\right)}^{2}}{2{\tau }^{2}}-{\rm{\Gamma }}{\sin }^{2}\displaystyle \frac{\pi ({t}_{i}-{t}_{j})}{T}\right]\\ & & +\left[{\sigma }_{i}^{2}+{\sigma }_{\mathrm{jit}}^{2}\right]{\delta }_{i,j}\end{array}\end{eqnarray} \tag{ 2 }$$

where C_i,j represents the covariance matrix, δ_i,j is the Kronecker delta function, h is the covariance amplitude, t_i is the time of ith observation, τ is the correlation timescale, Γ is the relative importance between the squared exponential and periodic parts of the kernel, T is the period of the covariance, σ_i is the reported uncertainty, and σ_jit is a jitter term. Our likelihood function is as follows:

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }=-\displaystyle \frac{N}{2}\mathrm{log}2\pi -\displaystyle \frac{1}{2}\mathrm{log}| {\boldsymbol{C}}| -\displaystyle \frac{1}{2}{{\boldsymbol{r}}}^{{\rm{T}}}{{\boldsymbol{C}}}^{-1}{\boldsymbol{r}}\end{eqnarray} \tag{ 3 }$$

where ${ \mathcal L }$ is the likelihood function, N is the total number of RV data points, C represents the covariance matrix, and r is the residual the observed RV minus the model RV.

We ran an MCMC analysis (similar to that described in Section 4.3) on the TESS light curves to constrain the posterior distribution of the various hyperparameters. As we mentioned earlier, we could not robustly detect the stellar rotation period in the TESS light curve from a periodogram analysis. This is may be due to the low activity of the host star. Correspondingly, our GP model of the TESS photometry yielded a very broad posterior distribution on various hyperparameters. We then tested whether the addition of the HIRES S_HK index gave better constraints on the hyperparameters. S_HK was sparsely sampled and only varied very mildly with an rms of about 0.004. Consequently, S_HK did not significantly improve the constraints on the GP hyperparameters.

Therefore, we chose to let the hyperparameters float freely in our final GP analysis of the RV data set. The exception is that we did limit T the covariance period (a proxy for stellar rotation period) between 1 and 200 days to avoid inference with the 0.47 day planet b. Understandably, without an imposed prior on the various hyperparameters, the flexibility of this GP model subsumed the signal of the candidate 16 day planet c. Planet c only has an upper limit of K_c < 4.3m s⁻¹ at a 95% confidence level in our GP model. In fact, the lowest BIC GP model only included the signal from planet b (Figure 9). An MCMC analysis showed that the posterior distribution of the semiamplitude of planet b is K_b = 3.59 ± 0.75m s⁻¹ which agrees with the ${K}_{b}={3.30}_{-0.59}^{+0.58}$ m s⁻¹ found by RadVel. Since the GP model has far more complexity than the RadVel model and the K value came out less well constrained, we adopted the results from RadVel for further analysis in this work.

Zoom In Zoom Out Reset image size

Before any discussion, we would like to loosely define USPs as terrestrial planets ( < 2R_⊕) with an orbital period of less than one day (the prevailing definition in the literature) as well as those with an insolation level larger than 650F_⊕ (see Table 4 for the complete list). 650F_⊕ is an empirical boundary, identified by Lundkvist et al. (2016) and Sanchis-Ojeda et al. (2014), beyond which photoevaporation is so strong that any Neptune-sized planets are quickly stripped down to their rocky cores by the strong stellar irradiation (Figure 10), thereby creating a "hot Neptune desert."

Zoom In Zoom Out Reset image size

Table 4. Additional Planetary Companions in USP and Hot-Neptune Systems Ranked by Planet Mass

System	Additional Planets Detected?	Comments	Reference
Kepler-78	N	activity-induced RV variation (50 m s−1 peak to peak, 12.5 day periodicity) prevents detection of additional planets	Howard et al. (2013)
GJ 1252	Y	15 m s−1 drift over 10 days likely due to an additional planet	Shporer et al. (2020)
K2-229	Y	two additional transiting planets on 8 and 31 day orbit	Santerne et al. (2018)
LTT 3780	Y	one additional transiting planet on 12 day orbit	Cloutier et al. (2020)
K2-36	Y	one additional transiting planet on 5 day orbit	Damasso et al. (2019)
Kepler-10	Y	one additional transiting planet on 45 day orbit	Dumusque et al. (2014)
TOI-1444	Y	one additional nontransiting planet on 16 day orbit	This work
CoRoT-7	Y	one additional nontransiting planet on 3.7 day orbit	Haywood et al. (2014)
K2-141	Y	one additional nontransiting planet on 7 day orbit	Malavolta et al. (2018)
HD 3167	Y	one additional nontransiting planet on 8.5 day orbit and 1 transiting planet on 29 day orbit	Christiansen et al. (2017)
HD 80653	Y	RV drift suggests one additional planet of $0.37\pm 0.08{M}_{\mathrm{Jup}}{\left(a/{au}\right)}^{2}$	Frustagli et al. (2020)
K2-131	N	activity-induced RV variation (60 m s−1 peak to peak, 3 day periodicity) prevents detection of additional planets	Dai et al. (2017)
K2-291	N	activity-induced RV variation (30 m s−1 peak to peak, 19 day periodicity) prevents detection of additional planets	Kosiarek et al. (2019)
WASP-47	Y	two additional transiting planets with a 4 day giant planet; 1 nontransiting giant planet on 600 day orbit	Vanderburg et al. (2017)
K2-106	Y	one additional nontransiting planet on 13 day orbit	Guenther et al. (2017)
55 Cnc	Y	four additional nontransiting planets between 14 and 5000 days including a close-in giant planet	Bourrier et al. (2018)
HD 213885	Y	one additional nontransiting planet on 5 day orbit	Espinoza et al. (2020)
K2-100	N	activity-induced RV variation (100 m s−1 peak to peak, 4.3 day periodicity) prevents detection of additional planets	Barragán et al. (2019)
LTT 9779	N	20 days of RV baseline	Jenkins et al. (2020)
TOI-849	N	>400 days of RV baseline	Armstrong et al. (2020)

Download table as:  ASCIITypeset image

Dai et al. (2019) performed a uniform analysis of all transiting USPs that also have RV mass constraints. In particular, they used Gaia parallax information to better constrain the host stellar properties. The increased precision on the stellar radius translated to increased precision on planetary radii. Moreover, the mean stellar density from Gaia and spectroscopy further disentangled degeneracies in transit modeling. The results significantly improved constraints on planetary radii. For example, the radius constraint of Kepler-78 b improved from R_p = 1.20 ± 0.09R_⊕ in Howard et al. (2013) to ${R}_{p}={1.228}_{-0.019}^{+0.018}{R}_{\oplus }$ with Gaia. Dai et al. (2019) also applied a Gaussian process model uniformly to all USP planets that required mitigation of stellar activity contamination in RV data sets. The improved mass and radius constraints on the sample of 10 USPs revealed a prevalence of 35%Fe-65%MgSiO₃ Earth-like composition.

In this work, we applied the same set of analysis to TOI-1444b and interpret the resultant mass and radius constraint along with other USP planets. The interest in USP planets has boomed in recent years; the number of USP RV mass measurements increased from 10 in Dai et al. (2019) to 17 at the time of writing of this work. Moreover, the inclusion of Gaia parallax information and Gaussian process modeling has become a standard practice in these new USP papers (Table 4). This puts the new USPs reported by different groups on relatively equal footing and ready for comparison. In the upper panel of Figure 11, we show the mass and radius of all USP planets with RV mass measurements. At a glance, USPs seem to cluster around an Earth-like composition. To quantify the compositions, we adopted a simple two-layer model where planets have an iron core and a MgSiO₃ mantle. We used the procedure described by Zeng et al. (2016) to convert the measured mass and radius of the planet to a Fe-MgSiO₃ ratio. For an individual planet, the confidence interval is relatively wide, e.g., TOI-1444b can have 43 ± 30% of its mass in iron given our mass and radius measurement. However, as an ensemble, the USP planets generally cluster around the Earth-like composition with a weighted mean of 32 ± 4% mass in an iron core. This is consistent with the general picture that planets at the lower peak of the bimodal radius distribution are predominantly rocky (Rogers 2015; Dressing et al. 2015; Fulton et al. 2017).

Zoom In Zoom Out Reset image size

We note that we have excluded two USPs from the averaging of the iron core mass fraction. The mass of TOI-561 reported by Weiss et al. (2021) and Lacedelli et al. (2021) differ by almost a factor of two. More RV data are being taken to resolve this discrepancy (C. Brinkman et al, 2021 in preparation). By focusing on the USP planets, i.e., planets that are most strongly irradidated, we hope to probe the exposed planetary cores directly without worrying about the degeneracy introduced by a planetary atmosphere. This assumption has held up well for most USPs in our sample. We examined the composition of USPs against the insolation level. If a substantial atmosphere were present on a USP, one may expect the scale height to vary strongly with insolation level. The scale-height variation would have translated to a correlation between insolation level and the inferred planetary composition. However, no correlation between insolation and iron core mass fraction was found (Figure 11). That said, we did exclude 55 Cnc e, one of the largest and coolest USPs, from our analysis. The strong phase offset seen in Spitzer observations of 55 Cnc e (Demory et al. 2016) demands a strong heat circulation between the day and the night side of the planet that, as Angelo & Hu (2017) argue, requires the presence of an atmosphere on 55 Cnc e. The larger planet mass and the cooler equilibrium temperature of 55 Cnc e may help retain its atmosphere compared to other USPs in the current sample.

What is the relation between USPs and the recently reported ultrahot Neptunes K2-100 b (Mann et al. 2017), LTT 9779 b (Jenkins et al. 2020), and TOI-849 b (Armstrong et al. 2020)? Are they the same group of planets that only differ in size? In this section, we will show that USPs and ultrahot Neptunes differ in planetary/stellar properties and system architecture, suggesting that they are probably two distinct groups of planets.

We identified a sample of three very hot Neptunes between 3-5 R_⊕ and with an insolation level >2000F_⊕ (Figure 10). Again >650F_⊕ is the empirical boundary of the "hot Neptune desert" proposed by Lundkvist et al. (2016). >2000F_⊕ securely put these planets within the "desert". Currently, there are three confirmed hot Neptunes straddling this so-called "hot Neptune desert": LTT 9779 b(Jenkins et al. 2020), TOI-849 b(Armstrong et al. 2020), and K2-100b (Mann et al. 2017) which have well-measured mass and radii. We put these hot Neptunes in the same mass–radius diagram with the USPs (Figure 11 lower panel). Clearly, the USPs and the hot Neptunes occupy different parts of the mass–radius parameter space. The USPs are all below 10M_⊕, even though RV mass measurements bias toward the detection of heavier planets (a point we will return to later). The USPs also cluster around an Earth-like composition of iron-rock mixture. However, the three close-in hot Neptunes are near or above 100% water composition line in the mass–radius diagram. This indicates the presence of a substantial H/He atmosphere or a water envelope for these planets. If we assume that they have Earth-like rocky cores, a 1%–5% mass fraction H/He envelope is needed to reproduce the observed mass and radius. Recent Spitzer phase-curve and secondary eclipse observations of LTT 9779 b pointed to the presence of a heavy molecular weight (e.g., CO) atmosphere (Crossfield et al. 2020; Dragomir et al. 2020).

Winn et al. (2017) showed that the metallicity of USP host stars is solar-like which is similar to the longer-period sub-Neptune planets commonly discovered by Kepler (blue histogram in Figure 12). By contrast, hot Jupiters are known to preferentially occur around metal-rich systems (orange histogram in Figure 12). Interestingly, the three hottest Neptunes all have similarly metal-rich host stars: LTT 9779 with [Fe/H] = 0.27 ± 0.03 (Jenkins et al. 2020), K2-100 with [Fe/H] = 0.22 ± 0.09 (Mann et al. 2017), and TOI-849 with [Fe/H] = 0.19 ± 0.03 (Armstrong et al. 2020). Even though the sample size is only three, the distribution of hot Neptunes is beginning to show a stark distinction from the solar-like distribution of USPs, and bears more resemblance to the metal-rich environments of hot Jupiters. This result is reminiscent of the previous work by Dong et al. (2018) who showed an even more statistically robust preference for hot Neptunes to occur around metal-rich host stars if we consider orbital periods as long as 10 days.

Zoom In Zoom Out Reset image size

Another feature that differentiates USPs from hot Neptunes is the presence of additional planetary companions. Out of the 17 well-characterized USPs, 14 have other planetary companions in the same system discovered by transits or radial-velocity follow up (Table 4). We further note that the RV data sets of the three systems where the USPs remain the only detected planets (K2-131, Kepler-78, and K2-291) are heavily contaminated by activity-induced RV variation that amounts to tens of m s⁻¹. We applied RVSearch to these systems to constrain the sensitivity of the existing RV data set to additional planetary signals (similar to Figure 7). RVSearch showed that the activity-induced RV variation could easily inhibit the detection of injected planetary signals with <15M_⊕ within 1 au. Turning our attention to the hot Neptunes, among the three hottest Neptunes, none of them have additional detected planets even though some of RV monitoring baselines extended more than 100 days (Table 4). This hints that hot Neptunes tend to be "lonely" similar to the hot Jupiters.

To sum up, USPs tend to be rocky in composition while hot Neptunes possess substantial H/He atmospheres or other volatile envelopes. USPs occur in stars with solar-like abundances, while hot Neptunes are preferentially found in metal-rich systems. USPs tend to have other close-in planetary companions while hot Neptunes tend to be lonely. These properties distinguish the two groups. Moreover, the properties of the hottest Neptunes are quite similar to their bigger cousins hot Jupiters. This similarity actually extends to many Neptune-sized planets between 2–6 R_⊕ and orbital periods of 1–10 days as pointed out by Dong et al. (2018). Finally, we close the section by noting that USPs, hot Jupiters, and <10 day hot Neptunes all occur at about a 1% level around Sun-like stars (e.g., Sanchis-Ojeda et al. 2014; Cumming et al. 2008; Dong et al. 2018). However, the hot Neptunes on sub-day orbits seem rarer. The Kepler USPs from Sanchis-Ojeda et al. (2014) represent a uniformly studied sample amenable to occurrence rate studies. Among the 67 confirmed/validated USP candidates, only 1 (KOI-3913) has a radius in the Neptune regime >3R_⊕. The rest are all below 2R_⊕ (Figure 12). Even though this is still small-number statistics, the Kepler USP sample suggests that hot Neptunes on sub-day orbits are probably as rare as 0.1%–0.01% around Sun-like stars. Future TESS occurrence rate studies will firm up this number. Another system worth mentioning is NGTS-4b (West et al. 2019), a 3.2R_⊕ hot Neptune on a 1.3 day orbit around a K star. It is less strongly irradiated compared to K2-100b, TOI-849 b, and LTT 9779 b, right on the edge of the "hot Neptune desert". It shows many similarities with K2-100b, TOI-849 b, and LTT 9779 b: likely enshrouded by an atmosphere, having no other planetary companions. However, the metallicity of its host star was reported to be surprisingly low: reported in units of total metal [M/H]=-0.28±0.10. We encourage future observation to confirm the reported low metallicity.

The top contenders for USP formation theory invoke the secular interaction between USPs and longer-period planets (Petrovich et al. 2019; Pu & Lai 2019). Consider a USP that formed initially on longer orbital periods in the 2-10 days range. This is beyond the dust sublimation radius and where many Kepler planets are found today. If there are other planets in the same system with enough angular momentum deficit (AMD), i.e., the angular momentum difference between a system with coplanar, circular orbits and a system with eccentric, mutually inclined orbits, secular interaction between the planets could shuffle the AMD around. Since the progenitors of USPs are the innermost planet of their system, they have the lowest angular moment per unit mass. The same AMD could thus induce a significantly eccentric and inclined orbit. The augmented tidal interaction with the host star at periastron could then shrink the orbit of the USPs. This secular scenario has gained observational support as the USPs are indeed observed on more inclined orbits compared to other Kepler planets and USPs often have a longer orbital period ratio relative to their neighbors—likely a consequence of orbital decay (Dai et al. 2018; Steffen & Coughlin 2016). Here, we point out further observational support of the secular theory, i.e., a USP must have additional planets to initiate the secular interaction and provide enough AMD. In Table 4 we showed that USPs almost always have additional planetary companions (14/17). For 55 Cnc, the existing RV data set is big enough to reveal orbital eccentricities of the various planets, thereby allowing us to gauge whether the amount of AMD in a system could alter the USP orbit significantly. Indeed, Hansen & Zink (2015) showed that the architecture of 55 Cnc does contain the requisite AMD and is in general consitent with the secular formation scenario. Future extensive RV follow up of USPs may extend this AMD test to other systems.

Hot Neptunes, including those on sub-day orbits, show striking similarities with hot Jupiters: (1) they favor metal-rich environments, (2) they tend to be the only planet in a system, and (3) their occurrence rate sums up to about a 1% occurrence rate within 10 days around Sun-like stars (Dong et al. 2018). Another similarity between hot Jupiters and hot Neptunes is that many hot Neptunes were also observed to be on misaligned orbits characterized by large stellar obliquities, e.g., Kepler-63 (Sanchis-Ojeda et al. 2013b), HAT-P-11 (Sanchis-Ojeda & Winn 2011), and WASP-107 (Dai & Winn 2017; Rubenzahl et al. 2021). For hot Jupiters, a large spin–orbit angle has traditionally been interpreted as a signpost of a dynamically hot formation scenario that tilted planet orbit while also triggering orbital decay (see the view by Dawson & Johnson, 2018). The large obliquity of many hot Neptunes is suggestive of a formation channel similar to that of the hot Jupiters. It will be instructive to extend obliquity measurements to ultrahot Neptunes such as K2-100, LTT 9779, and TOI-849, although this is technically challenging with current generation spectrographs. Berger et al. (2020) also presented evidence that hot Neptunes are preferentially found around evolved stars. This could be an indication that the hot Neptunes migrated as a result of host stellar evolution.

Two recent planet formation theoretical works may also shed light on the distinction between USPs and hot Neptunes. Adams et al. (2020) performed a simple analysis that optimized the total energy of a pair of forming planets assuming the conservation of angular momentum, constant total mass reservoir and fixed orbital spacings. They found that when the total available mass is low (≲40M_⊕, a number that depends on stellar mass, semimajor axis, etc.), the energy-optimized state is an equal partition of mass between two competitively growing planets. This thus tends to create multi-planet systems with intra-system uniformity that was seen in many observed sub-Neptune multi-planet systems (Millholland et al. 2017; Wang 2017; Weiss et al. 2018). However, when the total mass is high enough, the system switches to a different optimized state in which the mass is concentrated in one of the planets, thereby creating a dominant, possibly lonely planet that may undergo runaway accretion. Whether a system does go into the runaway accretion state, as pointed out by Lee (2019) and Chachan et al. (2021), further depends on the local hydrodynamic flow conditions, the local opacity, and the time of core emergence, i.e., whether gas is still present in the disk. These theories put the USPs and hot Neptunes into perspective: in a metal-rich (high [Fe/H]) disk, the planetesimals assemble quickly and grow large enough to accrete the remaining gas in the disk (see also Dawson et al. 2015). If the total planetesimal mass is above some threshold, one particular planet embryo grows to be the dominant planet in a system as Adams et al. (2020) would predict. In short, metal-rich systems tend to breed hot Neptunes and hot Jupiters. By contrast, the rocky USP planets could be the product of late assembly in a gas-depleted disk in which core assembly proceeds oligarchically and slowly. This explains why their occurrence does not correlate strongly with host-star metallicity. Adams et al. (2020) would further predict that USP planets come in multi-planet systems where the planets are similar in size.

Table 4 shows that there are three USPs that also have detected giant planet companions. They are 55 Cnc with a 14.6 day giant planet (Butler et al. 1997), WASP-47 with 4.2 and 600 day giant planets (Becker et al. 2015; Neveu-VanMalle et al. 2016), and HD 80653 with a strong RV trend likely indicative of a giant planet $0.37\pm 0.08{M}_{\mathrm{Jup}}{\left(a/{au}\right)}^{2}$ (Frustagli et al. 2020). The attentive reader might have already guessed that these three systems are also the most metal-rich among the USP sample with [Fe/H] = 0.31 ± 0.04, 0.38 ± 0.05, and 0.28 ± 0.05, respectively (red points in Figure 12). These three USPs also happen to be the more massive ones among the USP sample (Figure 12). One is tempted to ask: are USPs in metal-rich environments also more massive ones? We use the ${M}_{\star }^{* }{10}^{[\mathrm{Fe}/{\rm{H}}]}$ as a proxy for the surface density of solid material available to USP growth. Dai et al. (2020) applied a minimum-mass extrasolar nebula framework to Kepler sub-Neptune planets which suggeested that the disk solid mass displayed a linear scaling with host-star mass even within the innermost 1 au. ${M}_{\star }^{* }{10}^{[\mathrm{Fe}/{\rm{H}}]}$ may be a reasonable proxy for the solid-disk density. We found ${M}_{\star }^{* }{10}^{[\mathrm{Fe}/{\rm{H}}]}$ and USP mass do show a positive correlation with a p value of 0.04 in a Pearson test (Figure 12). This correlation may suggest that the surface density of solids in the disk directly controls the overall availability of solids, the assembly rate, and eventually manifests as the size of planetary cores that emerge out of the disks.

As we argued above, USPs are probably a distinct group compared to hot Jupiters and hot Neptunes given the differences in host-star metallicity and planet architecture. USPs are most likely cores of planets that escaped runaway accretion in the first place, i.e., super Earths or sub-Neptunes. The rate of atmospheric erosion, particularly photoevaporation, shows a steep dependence on planet mass. It is suppressed by orders of magnitude for giant planets (e.g., Murray-Clay et al. 2009). Briefly, this is because the gravitational potential well deepens with planet mass, but the heating efficiency of XUV irradiation does not. For example, the H/He envelopes on a 5 M_⊕ planet can be easily stripped in 100 Myr at 0.1 au around a G-type star, but the mass loss timescale quickly shoots up to more than a Hubble time for planetary cores with >10–15 M_⊕ depending on the insolation the planet receives (Wang & Dai 2019). In summary, planets heavier than Neptune can only lose a small fraction of their envelope over their lifetime; thus they are unlikely to be stripped down to a rocky core that is observed as a USP planet today.

Therefore, the sample of USP planets help to place an upper limit on the threshold of runaway accretion. RV follow-up observations are heavily biased toward the detection of more massive planets. However, this bias works in our favor as we are probing an upper limit in mass. If one assumes a radiative outer layer, the critical mass for runaway accretion only has a logarithmic dependence on the local disk properties (e.g., Rafikov 2006). In other words, it is not sensitive to the location where the planets formed and is estimated to be about 10 M_⊕ (Rafikov 2006). If one more carefully accounts for the envelope opacity, the threshold mass may be more variable between 2 and 8 M_⊕ (Lee & Chiang 2016). The lower the local gas opacity, the lower the mass threshold for runaway accretion (see also Chachan et al. 2021). In this framework, the super-puff, low density, gas-rich planets with <5M_⊕, >5 R_⊕, likely formed beyond the snowline where the opacity is low (e.g., Kepler-51 Libby-Roberts et al. 2020). Coming back to the current sample of 17 USPs, they are mostly consistent with exposed rocky cores. They seem to probe the upper limit of the threshold for runaway accretion. This seems to be at the predicted 8 M_⊕ for the high-opacity inner disk (∼0.1 au) where USPs have likely formed. We note that the only USP that is suspected to have some atmosphere is 55 Cnc e (Angelo & Hu 2017) right at ∼8M_⊕. This putative current-day atmosphere on 55 Cnc e was also suggested to be secondary in nature. 55 Cnc e was neither detected in Lyα (Ehrenreich et al. 2012) or metastable He (Zhang et al. 2021) transit observations.