© The Authors 2026

1. Introduction

Soon after the first discoveries, a deficit of Neptunian planets at very short orbital periods was identified (e.g. Mazeh et al. 2005, 2016; Lecavelier Des Etangs 2007), a feature commonly referred to as the Neptunian desert. This deficit reflects a genuine scarcity of planets rather than an observational bias. Occurrence-rate analyses have shown that the desert remains highly underpopulated even when survey incompleteness is taken into account (Castro-González et al. 2024a; Cui et al. 2026). However, its origin remains an open problem (Matsakos & Königl 2016; Bourrier et al. 2018; Owen & Lai 2018).

Atmospheric escape can erode gas-rich planets at short orbital separations and has been invoked as a plausible contributor to the desert (Owen & Lai 2018). In parallel, high-eccentricity tidal migration (HEM; e.g. Ford & Rasio 2008; Correia et al. 2011) has been widely discussed as a mechanism shaping the desert. If outer perturbers drive planets to sufficiently small periastrons, tidal dissipation can lead either to disruption or to rapid circularisation onto compact final orbits, with characteristic periods comparable to those of the desert edge (e.g. Liu et al. 2013; Mazeh et al. 2016; Matsakos & Königl 2016). Today, it is increasingly thought that an interplay between these mechanisms shapes the properties of close-in Neptunes (e.g. Bourrier et al. 2018; Owen & Lai 2018), with contributions from tidal inflation and orbital decay (e.g. Hallatt & Millholland 2026a,b).

Recently, the observational picture of the ‘exo-Neptunian landscape’ has sharpened. Using occurrence-rate calculations to correct for survey biases, Castro-González et al. (2024a) redefined the boundaries of the Neptunian desert and show that it is bordered by an overdensity of planets concentrated at P_orb ≃ 3.2 − 5.7 d, which separates the desert from the more sparsely populated, longer-period Neptunian savanna (Bourrier et al. 2023). This feature, referred to as the Neptunian ridge, resembles the hot-Jupiter pile-up, suggesting that the closest-in Neptune- to Jupiter-size planets share related evolutionary pathways (e.g. Castro-González et al. 2024a). The origin of the pile-up remains debated, but HEM is widely considered to play a key role (see Dawson & Johnson 2018, and references therein). Recent observational evidence points to a similarly important contribution in the ridge (e.g. Correia et al. 2020; Bourrier et al. 2023, 2025; Doyle et al. 2025; Vissapragada & Behmard 2025; Handley et al. 2025; Yee et al. 2025; Handley & Batygin 2026).

In this work we revisit the HEM scenario in light of the updated, bias-corrected geometry of the Neptunian desert, ridge, and savanna inferred by Castro-González et al. (2024a).

2. Tidal survival and circularisation

In HEM scenarios, planets are driven onto highly eccentric orbits through mechanisms such as secular perturbations and planet–planet scattering. Depending on their periastron distance, they can be disrupted, collide with the star, be ejected, or survive and circularise. Their fate is governed by tidal interactions near periastron, which impose a fundamental constraint on the resulting population of close-in survivors.

Planets that approach their host stars too closely are tidally stripped and can be fully disrupted. This defines a minimum periastron distance set by the tidal radius (Roche 1849),

$$\begin{array}{c}\hfill r_{\mathrm{tide}}=\eta R_{\mathrm{p}}{\left(\frac{M_{★}}{M_{\mathrm{p}}}\right)}^{1/3},\end{array}$$(1)

where η is a dimensionless encounter parameter that depends on the planetary internal structure and on the hydrodynamics of the tidal interaction. This condition translates into a minimum circularised orbital period after tidal dissipation. Because the limit depends on the planet’s bulk density, it can be mapped onto the period–radius plane through a mass–radius relation, defining a tidal survival boundary, P_tide(R_p), for planets delivered through HEM (Matsakos & Königl 2016).

Survival alone does not guarantee compact post-HEM orbits. Tidal circularisation is highly sensitive to periastron distance, so planets that avoid very close passages may survive but fail to circularise on system timescales (e.g. Fabrycky & Tremaine 2007). Tidal evolution therefore operates most efficiently in a narrow periastron range near the disruption threshold.

Appendix A.1 gives the explicit expressions used to evaluate P_tide(R_p), and Appendices A.2 and A.3 show illustrative circularisation scalings for the Lidov-Kozai and secular-chaos channels. Other pathways, such as planet–planet scattering, are expected to show similar behaviour, since the outcome is set mainly by the tidal physics rather than by the specific excitation mechanism.

3. A tidal origin for the Neptunian ridge

We confronted the HEM tidal constraints with the distribution of short-period planets, focusing on the empirically determined, bias-corrected boundaries of the Neptunian desert and ridge (Fig. 1; Castro-González et al. 2024a). A key feature of the updated landscape is the flattening of the lower desert boundary at R_p ≃ 1.8 R_⊕ over P_orb ≃ 0.5 − 1.2 d, which we refer to as the ultra-short-period plateau. We restricted our analysis to R_p > 1.8 R_⊕ since rocky planets possess significant material strength. Moreover, the shortest-period planets systematically approach the disruption limit (Castro-González et al. 2024a, 2025), likely requiring orbital decay or alternative migration channels.

Fig. 1. Confrontation between tidal survival limits and the close-in planet population. (a) Period–radius plane. The solid black curve traces the Neptunian desert boundary, and the vertical dashed line marks the ridge limit (Castro-González et al. 2024a). Dashed red and blue curves show tidal survival boundaries Ptide(Rp; q) for the q = 0.84 and q = 0.16 mass–radius envelopes, adopting η = 4.5, and the shaded region indicates the survival band. Dotted curves show representative $P_{\mathrm{circ}}^{max}$ tracks with $Q{\prime }_{\mathrm{p}}={10}^5$ and M★ = 1 M⊙ (Owen & Lai 2018): orange and red correspond to tsys = 8 and 4 Gyr with a0 = 1.5 au along the q = 0.84 envelope, and blue to tsys = 8 Gyr with a0 = 1.0 au along the q = 0.16 envelope. Dash-dotted dark magenta and pink curves show illustrative secular-chaos circularisation curves (e.g. Matsakos & Königl 2016), computed using M(Rp; q = 0.5) and perturbing-companion masses of 0.7 and 0.4 MJ, respectively. (b) Mass–radius plane with the quantile envelopes M(Rp; q).

3.1. Mass–radius envelopes and tidal survival limits

To map the theoretical tidal limits onto the (P, R_p) plane, a mass–radius relation is required. We constructed this relation empirically using planets with measured masses from the NASA Exoplanet Archive, retaining only objects with sufficiently precise measurements (σ_Mp/M_p ≤ 20% and σ_Rp/R_p ≤ 10%) to ensure fractional density uncertainties ≲35% (for reference, the 1σ density scatter is ≃2.3 g cm⁻³ in the 10 − 30 M_⊕ regime).

At fixed radius, envelopes of log M are estimated using quantile statistics, yielding relations M_p(R_p; q), where q denotes the cumulative quantile level. We adopted the quantile pair q = 0.8413 and q = 0.1587, which correspond to the ±1σ levels of a normal distribution. These envelopes represent the high- and low-density ends of the gas-rich planet population. Together, they define a finite tidal survival band in the period–radius plane,

$$\begin{array}{c}\hfill P_{\mathrm{tide}}(q=0.8413)<P\le P_{\mathrm{tide}}(q=0.1587),\end{array}$$(2)

which captures the range of the shortest post-HEM orbits reachable for the bulk of systems without tidal disruption once the dispersion in density is taken into account¹.

3.2. A coupled HEM origin of the ridge and desert boundary

The slope of the desert boundary from the sub-Neptune to the super-Neptune/sub-Saturn regime (1.8 R_⊕ ≲ R_p ≲ 6 R_⊕) follows directly from the HEM tidal survival limit when evaluated along the high-density envelope, P_tide(q = 0.84). This agreement is not enforced by construction but instead emerges from the tidal survival criterion itself. Adopting a single representative value of the tidal encounter parameter, η = 4.5, which sets the overall period offset, the predicted boundary reproduces the bias-corrected desert edge over this full radius range. We note that η should not be interpreted here as a universal disruption constant but rather as an effective parameter for the onset of strong tidal disruption within the considered framework. In particular, the derived value is tied to the definition of the desert, the mass–radius relation, and the quantile envelopes. Still, the effective η inferred here is of the same order of magnitude as those reported in previous works, including the η = 2.7 from hydrodynamic calculations for Jupiter-like planets (Guillochon et al. 2011) and the η = 3.5 from a related HEM study (Matsakos & Königl 2016).

In the Jovian regime, additional effects such as radius inflation or post-circularisation orbital decay are likely required to account for residual deviations from the measured boundary, although the overall boundary shape and slope transitions remain broadly consistent with the tidal survival limit. Because a single value of η closely reproduces the overall desert geometry, we adopted it consistently to construct the tidal survival band. Remarkably, this band traces the period interval of the Neptunian ridge, 3 ≲ P_orb ≲ 6 d (Fig. 1). The ridge therefore emerges as the natural locus of the shortest-period HEM survivors once the observed density dispersion is taken into account. We note that this interpretation assumes that the density distribution relevant for HEM was broadly similar to that observed today. If some planets migrated at younger ages or were temporarily inflated during HEM, the low-density tail and thus the tidal survival band could have been broader, which would also affect the value of η.

Tidal dissipation increases very steeply towards small periastron distances, so efficient circularisation is concentrated among trajectories that pass close to the disruption threshold. As a consequence, successful HEM outcomes typically accumulate just beyond the survival limit (e.g. Fabrycky & Tremaine 2007; Giacalone et al. 2017; Owen & Lai 2018). This concentration is intrinsically stronger for lower-density planets, whose larger disruption radii shift the survival boundary to longer periods. Therefore, the combination of the steep tidal dependence on periastron distance and the finite density dispersion naturally confines surviving HEM outcomes to a narrow period interval that traces the Neptunian ridge. Figure 1 also shows representative maximum circularisation periods for the Lidov-Kozai channel and illustrative secular-chaos circularisation curves. While these first-order estimates are model-dependent and intended only as qualitative guides, their agreement with the desert–ridge structure supports the notion that HEM plays a major role in setting its geometry.

At larger radii, the same tidal framework predicts an overdensity of HEM-migrated planets in the hot-Jupiter pile-up (Fig. 1). At smaller radii, the HEM-accessible region is embedded within the high-occurrence sub-Neptune population, making any overdensity harder to isolate and requiring additional observables to determine a HEM contribution.

3.3. Tidal survival and clustering in the period–density plane

We further tested the HEM tidal survival formalism by projecting the close-in Neptunian population onto the period–density plane. Because the tidal survival condition depends on bulk density, it translates into a simple physical boundary in this space. For this analysis, we adopted η = 4.0, chosen so that the bulk of the planets lie outside the disruption region (Fig. 2). This value is slightly lower than the η = 4.5 adopted in Sect. 3.2, which is not unexpected since the two estimates rely on different assumptions regarding the desert boundary and the mass–radius distribution.

Fig. 2. Left: Period–density distribution of close-in Neptunian planets with precise masses and radii (4 ≤ Rp/R⊕ ≤ 8.5). The red curve shows the tidal survival limit adopting η = 4.0. Vertical dashed lines mark the desert–ridge and ridge–savanna boundaries (Castro-González et al. 2024a), and horizontal dotted blue and red lines indicate the low- and high-density quantiles and their intersections with the tidal limit. The grey-shaded region highlights the low-density depletion, which effectively limits typical HEM outcomes to the ridge interval. Right: Normalised density distribution, defining the savanna as Porb > 8 d.

The observed distribution shows two significant features previously reported in the literature: a group of ridge Neptunes are systematically denser than savanna planets, with a characteristic transition near 1 g cm⁻³ (Castro-González et al. 2024b), and the lower envelope of the density distribution rises towards shorter periods (the density ‘brink’; Bourrier et al. 2025). In Fig. 2, this lower envelope closely traces the predicted tidal survival limit. This agreement supports a scenario in which tidal disruption sets the minimum density required for planets to survive at a given orbital period under HEM. At fixed radius, higher-density planets can survive smaller periastron distances and reach shorter post-HEM periods, whereas lower-density planets are preferentially removed and are therefore absent at the shortest periods.

The same diagram also naturally explains the accumulation of HEM survivors in the ridge. The tidal survival curve intersects the low-density quantile of the Neptunian population at ≃6 d, closely matching the ridge–savanna boundary (Fig. 2). This result ultimately reflects the fact that the density distribution of Neptunian planets does not extend to sufficiently low values to effectively shift the minimum tidal survival period to orbital periods longer than the ridge interval. This depletion of very low-density planets is highlighted in Fig. 2. Overall, this density-limited survival window naturally explains both the confinement of HEM survivors to the ridge and the observed density gradient towards shorter periods.

A key open question is the origin of the dense Neptunes in the ridge. Castro-González et al. (2024b) suggest that the broad density spread reflects either intrinsically different internal structures, possibly linked to multiple formation pathways, or evolutionary effects such as atmospheric erosion. Within the HEM framework, another possibility is partial tidal stripping of initially low-density ridge planets during close periastron passages. Interestingly, in the right-hand panel of Fig. 2 we identify a non-Gaussian density distribution, with a main peak at ρ_p ≃ 0.6−0.8 g cm⁻³ and a secondary maximum near ρ_p ≃ 1.7 g cm⁻³. This feature is most evident in the ridge subsample and appears robust to the analyses presented in Appendix B. We therefore interpret it as a real density concentration associated with the ridge, whose origin deserves further investigation.

4. Discussion

The results presented in this work support a unified interpretation of the Neptunian ridge and the adjacent desert boundary within the HEM scenario, sharpening earlier HEM-based explanations. Matsakos & Königl (2016) showed that tidal-survival and secular-chaos curves can reproduce the basic desert geometry when combined with an empirical mass–radius relation. In a subsequent study, Owen & Lai (2018) showed that post-HEM circularisation typically occurs only slightly beyond the disruption boundary, implying a narrow range of final periods. Our key differences are that we confronted the HEM framework with the updated, bias-corrected desert–ridge boundaries and explicitly incorporated the observed dispersion in the mass–radius relation, which broadens the survival limit into a finite tidal survival band whose full extent naturally spans the Neptunian ridge and the hot-Jupiter pile-up. We note that recent models suggest that the HEM-driven circularisation of ridge Neptunes often causes runaway tidal inflation unless the planets are sufficiently metal-rich or tidal dissipation occurs mainly in their outer layers (Hallatt & Millholland 2026b). Atmospheric survival may therefore limit which HEM products remain in the ridge.

Independent observational trends further support a HEM origin for the Neptunian ridge overdensity. The excited eccentricities, polar orbits, and evidence of atmospheric erosion in several ridge planets are more consistent with a late, dynamically excited origin than with early, smooth disk-driven migration (e.g. Correia et al. 2020; Bourrier et al. 2023, 2025). Host-star properties point in the same direction: ridge and desert Neptunes preferentially orbit metal-rich stars relative to savanna planets (Vissapragada & Behmard 2025), consistent with a higher occurrence of massive companions capable of triggering HEM.

Finally, although HEM reproduces the desert–ridge geometry, it does not readily explain the savanna, which is mostly populated by nearly circular planets (Correia et al. 2020). In this context, Bourrier et al. (2025) propose that low-density Neptunes migrating early through the disk are strongly eroded at short periods, while those left farther out populate the savanna. By contrast, planets delivered later to the ridge through HEM would be less affected by early high-energy irradiation and could therefore survive at short periods (Bourrier et al. 2023; Attia et al. 2021).

5. Conclusions

We revisited the HEM scenario using the bias-corrected geometry of the Neptunian desert–ridge structure. The tidal disruption condition naturally reproduces the desert boundary, and the density dispersion broadens the disruption limit into a finite tidal survival band that traces the ridge. Since tidal dissipation increases steeply towards the disruption threshold, HEM survivors are expected to circularise just beyond this limit and thus cluster within the survival band, making the coupled desert boundary and ridge overdensity a natural outcome of HEM. The period–density plane provides an independent consistency check and reveals a non-Gaussian density distribution with a robust local maximum near ρ_p ≃ 1.7 g cm⁻³ that is primarily associated with ridge planets and whose origin warrants further investigation.

Acknowledgments

We thank the anonymous referee for a very careful, constructive, and insightful report, which significantly improved the clarity and overall quality of this work. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grant 51NF40_205606. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (project Spice Dune, grant agreement No 947634). ACMC acknowledges support from the FCT, Portugal, through the CFisUC project UID/04564/2025, with DOI identifier 10.54499/UID/04564/2025. This work has made use of the NASA Exoplanet Archive (Christiansen et al. 2025) and the Kepler DR25 catalogue (Thompson et al. 2018). This research made use of nep-des (available in https://github.com/castro-gzlz/nep-des).

References

Appendix A: Tidal survival and circularisation

A.1. Tidal disruption and minimum circularisation period

For M_p ≪ M_★, a planet of mass M_p and radius R_p orbiting a star of mass M_★ undergoes significant tidal stripping or disruption if its periastron approaches the characteristic tidal radius (Roche 1849; Chandrasekhar 1963),

$$\begin{array}{c}\hfill r_{\mathrm{tide}}=\eta R_{\mathrm{p}}{\left(\frac{M_{★}}{M_{\mathrm{p}}}\right)}^{1/3},\end{array}$$(A.1)

where the dimensionless parameter η encapsulates the effects of planetary internal structure and the hydrodynamics of the encounter (e.g. Guillochon et al. 2011). Following previous related studies, as a first approximation we assume a single value of η over the density range of gas-rich planets from sub-Neptunes to Jupiters (e.g. Matsakos & Königl 2016; Owen & Lai 2018).

In HEM scenarios, planets are driven to highly eccentric orbits (e → 1). If angular momentum is conserved, the final circularised semi-major axis is (Ford & Rasio 2006)

$$\begin{array}{c}\hfill a_{\mathrm{F}}\simeq 2r_{\mathrm{p}},\end{array}$$(A.2)

where r_p is the periastron distance. Requiring survival against tidal disruption implies

$$\begin{array}{c}\hfill a_{\mathrm{F}}\ge 2r_{\mathrm{tide}}.\end{array}$$(A.3)

Using Kepler’s third law, this condition translates into a minimum circularisation period,

$$\begin{array}{c}\hfill P_{\mathrm{F}}\ge 2\pi {\left(\frac{{(2r_{\mathrm{tide}})}^3}{G{M}_{★}}\right)}^{1/2}.\end{array}$$(A.4)

Substituting Eq. (A.1) into Eq. (A.4) yields

$$\begin{array}{c}\hfill P_{\mathrm{F}}\ge 2\pi {\left(\frac{{(2\eta R_{\mathrm{p}})}^3}{G{M}_{\mathrm{p}}}\right)}^{1/2}\propto {\left(\frac{R_{\mathrm{p}}^3}{M_{\mathrm{p}}}\right)}^{1/2}\propto {\rho }_{\mathrm{p}}^{-1/2},\end{array}$$(A.5)

showing that the tidal survival boundary is controlled primarily by the planet’s mean density and is independent of stellar mass. Once a mass–radius relation is specified, this condition defines a tidal survival boundary P_tide(R_p) in the period–radius plane.

A.2. Lidov-Kozai circularisation

Following Owen & Lai (2018), we estimate the maximum orbital separation from which tidal dissipation can circularise a highly eccentric orbit within a system age t_sys. This estimate is based on the equilibrium-tide formalism in the weak-friction approximation, appropriate for HEM driven by secular perturbations such as Lidov-Kozai cycles (e.g. Fabrycky & Tremaine 2007; Anderson et al. 2016).

Requiring that tidal dissipation circularises the orbit within t_sys yields an upper limit on the circularised semi-major axis,

$$\begin{array}{c}\hfill a_{\mathrm{F}}\le a_{\mathrm{circ}},\end{array}$$(A.6)

with

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill a_{\mathrm{circ}}& =0.05\mathrm{au}{\left(\frac{t_{\mathrm{sys}}}{1\mathrm{Gyr}}\right)}^{1/7}{\left(\frac{Q_p^{\prime }}{{10}^5}\right)}^{-1/7}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{2/7}\hfill \\ \hfill & \times {\left(\frac{a_0}{1\mathrm{au}}\right)}^{-1/7}{\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{J}}}\right)}^{-1/7}{\left(\frac{R_{\mathrm{p}}}{R_{\mathrm{J}}}\right)}^{5/7},\hfill \end{array}\end{array}$$(A.7)

where Q′_p is the effective planetary tidal quality factor and a₀ is the characteristic initial semi-major axis.

The corresponding maximum orbital period reachable by tidal circularisation is obtained via Kepler’s third law,

$$\begin{array}{c}\hfill P_{\mathrm{circ}}^{max}=2\pi {\left(\frac{a_{\mathrm{circ}}^3}{G{M}_{★}}\right)}^{1/2}.\end{array}$$(A.8)

Equation (A.8) provides a first-order estimate of the maximum final periods reachable on a timescale t_sys.

A.3. Secular-chaos circularisation

Following Matsakos & Königl (2016), we consider the secular-chaos circularisation scaling from Wu & Lithwick (2011). The characteristic post-circularisation semi-major axis is written as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill a_{\mathrm{c},\mathrm{SC}}& =0.03\mathrm{au}{\left(\frac{\alpha }{1/6}\right)}^{-3/5}{\left(\frac{M_{★}}{M_{\odot }}\right)}^{2/5}\hfill \\ \hfill & \times {\left(\frac{M_{\mathrm{pert}}}{M_{\mathrm{J}}}\right)}^{-1/5}{\left(\frac{M_{\mathrm{p}}}{M_{\mathrm{J}}}\right)}^{-2/5}\left(\frac{R_{\mathrm{p}}}{R_{\mathrm{J}}}\right),\hfill \end{array}\end{array}$$(A.9)

where α is a dimensionless normalisation parameter and M_pert is the mass of the perturbing companion. As in Matsakos & Königl (2016), we adopt α = 1/6, which sets the reference normalisation of the scaling. The corresponding orbital period is then obtained from Kepler’s third law,

$$\begin{array}{c}\hfill P_{\mathrm{SC}}=2\pi {\left(\frac{a_{\mathrm{c},\mathrm{SC}}^3}{G{M}_{★}}\right)}^{1/2}.\end{array}$$(A.10)

Appendix B: Density clustering diagnostics

We quantify the evidence for a persistent density feature near ρ_p ≃ 1.7 g cm⁻³ using complementary diagnostics applied consistently to the same parent sample. We select planets with measured masses and radii from the NASA Exoplanet Archive and apply the adopted radius range (4.0 ≤ R_p/R_⊕ ≤ 8.5) and precision cuts (σ_Mp/M_p ≤ 20%, σ_Rp/R_p ≤ 10%). Densities are computed as ρ_p = (M_p/R_p³) ρ_⊕, with uncertainties propagated as

$$\frac{{\sigma }_{\rho }}{{\rho }_{\mathrm{p}}}\simeq \sqrt{{\left(\frac{{\sigma }_{M_{\mathrm{p}}}}{M_{\mathrm{p}}}\right)}^2+{\left(3\frac{{\sigma }_{R_{\mathrm{p}}}}{R_{\mathrm{p}}}\right)}^2}.$$

Measurement uncertainties are propagated either analytically (for window occupancies) or via Monte Carlo resampling of ρ_p from per-planet Gaussian error models when assessing the stability of the diagnostics.

C.1. KDE peak stability near ρ_p ≃ 1.7 g cm⁻³

We analyse the density distribution using kernel-density estimation (KDE) in ρ_p with a Gaussian kernel. To avoid conclusions driven by a particular smoothing choice, we explore a grid of bandwidths spanning h = 0.06-0.22 g cm⁻³ (17 values), from under- to over-smoothed regimes. For each bandwidth, we evaluate whether the KDE profile exhibits a local maximum within a fixed narrow tolerance of the target density ρ_p = 1.7 g cm⁻³.

In the total sample, a local KDE maximum near 1.7 g cm⁻³ is recovered for 94% of the tested bandwidth choices. In the ridge subsample, the same peak is recovered for 100% of the bandwidth grid. When propagating measurement uncertainties through Monte Carlo resampling, the fraction of realisations retaining a local maximum at 1.7 g cm⁻³ is 0.58 for the total sample and 0.70 for the ridge sample. This bandwidth- and uncertainty-level robustness indicates that the ≃1.7 g cm⁻³ feature is unlikely to be an artefact of smoothing or measurement noise.

C.2. Gaussian-mixture evidence against a unimodal null

We fit one- and two-component Gaussian mixture models to the ρ_p distribution and compare them using the Bayesian information criterion (BIC), computing ΔBIC = BIC_K = 1 − BIC_K = 2. For the total sample, we obtain ΔBIC = 15.7, while for the ridge subsample ΔBIC = 3.35.

To assess whether these values can arise under a unimodal distribution, we perform a parametric bootstrap under a single-Gaussian null and recompute ΔBIC for each realisation. Under this unimodal null, the probability of obtaining ΔBIC at least as large as observed is p = 4.0 × 10⁻⁴ for the total sample and p = 3.6 × 10⁻³ for the ridge sample. Consistently, the unimodal null yields P(ΔBIC > 0)≃0.028 (total) and 0.009 (ridge).

When uncertainties are propagated via Monte Carlo resampling, the median ΔBIC remains positive for the total sample (median ≃12.5, 68% interval [7.2, 18.6]), while for the ridge subsample the distribution of ΔBIC spans zero (median ≃ − 0.1, 68% interval [-4.1, 4.4]), indicating weaker formal evidence for a fully separated bimodal decomposition in this subsample.

C.3. Local bump strength and window occupancy

We further quantify the prominence of the ≃1.7 g cm⁻³ feature using complementary local diagnostics. We define a ‘1.7-peak window’ of ρ_p ∈ [1.45,  2.05] g cm⁻³. A parametric local-excess test under the fitted unimodal null yields p = 0.236 (total) and p = 0.112 (ridge).

Accounting analytically for measurement uncertainties by summing the per-planet Gaussian probabilities of falling inside the window yields an uncertainty-aware expected number of planets in the window, $\mathbb{E}[N_{\mathrm{win}}]=19.97\pm 2.77$ (total) and 18.68 ± 2.60 (ridge). These values correspond to deviations of $z=(N_{\mathrm{obs}}-\mathbb{E}[N_{\mathrm{win}}])/\sigma \simeq 1.81$ (total) and z ≃ 2.04 (ridge).

We also define a KDE-based bump-strength statistic, S, at ρ_p = 1.7 g cm⁻³ as the logarithmic ratio between the KDE value at the target density and a local baseline estimated from symmetric sidebands. Using a non-parametric bootstrap, we obtain S = 0.364 with a 68% interval [0.116, 0.582] for the total sample and S = 0.477 with a 68% interval [0.242, 0.701] for the ridge subsample. By construction, S > 0 indicates a positive bump relative to the local baseline.

C.4. Unimodality diagnostics

Finally, we apply Hartigan’s dip test for unimodality. The resulting p-values are p = 0.936 (total) and p = 0.290 (ridge), indicating that unimodality cannot be rejected at high significance.

Overall, while formal unimodality tests do not require a fully separated bimodal structure, multiple independent diagnostics consistently reveal a statistically persistent local maximum near ρ_p ≃ 1.7 g cm⁻³, particularly within the ridge population.

All Figures

Fig. 1. Confrontation between tidal survival limits and the close-in planet population. (a) Period–radius plane. The solid black curve traces the Neptunian desert boundary, and the vertical dashed line marks the ridge limit (Castro-González et al. 2024a). Dashed red and blue curves show tidal survival boundaries Ptide(Rp; q) for the q = 0.84 and q = 0.16 mass–radius envelopes, adopting η = 4.5, and the shaded region indicates the survival band. Dotted curves show representative $P_{\mathrm{circ}}^{max}$ tracks with $Q{\prime }_{\mathrm{p}}={10}^5$ and M★ = 1 M⊙ (Owen & Lai 2018): orange and red correspond to tsys = 8 and 4 Gyr with a0 = 1.5 au along the q = 0.84 envelope, and blue to tsys = 8 Gyr with a0 = 1.0 au along the q = 0.16 envelope. Dash-dotted dark magenta and pink curves show illustrative secular-chaos circularisation curves (e.g. Matsakos & Königl 2016), computed using M(Rp; q = 0.5) and perturbing-companion masses of 0.7 and 0.4 MJ, respectively. (b) Mass–radius plane with the quantile envelopes M(Rp; q).

In the text

Fig. 2. Left: Period–density distribution of close-in Neptunian planets with precise masses and radii (4 ≤ Rp/R⊕ ≤ 8.5). The red curve shows the tidal survival limit adopting η = 4.0. Vertical dashed lines mark the desert–ridge and ridge–savanna boundaries (Castro-González et al. 2024a), and horizontal dotted blue and red lines indicate the low- and high-density quantiles and their intersections with the tidal limit. The grey-shaded region highlights the low-density depletion, which effectively limits typical HEM outcomes to the ridge interval. Right: Normalised density distribution, defining the savanna as Porb > 8 d.

In the text