© The Authors 2024

1 Introduction

One of the most intriguing features in the various findings in the realm of exoplanetary science is the scarcity of exoplanets that spans from sub-Jupiters to super-Earths at small orbital distances (period of ≲4 days). This is known as the hot sub-Jovian/ Neptune desert (e.g. Mazeh et al. 2016; Szabó & Kiss 2011). We adopt this nomenclature here. Mazeh et al. (2005) first pointed out that a negative linear relation exists between the planetary mass and the orbital period. As more planets were discovered in the following years, several other authors examined this tendency (e.g. Szabó & Kiss 2011; Beaugé & Nesvorný 2013) and broadened this correlation to the planetary radius (e.g. Beaugé & Nesvorný 2013; Mazeh et al. 2016). In particular, Szabó & Kiss (2011) found a sparsely populated region in the P− M_p diagram with P ≤ 2.5 days, spanning 0.02–0.8 M_J, which they termed the sub-Jovian desert. Similar findings were also obtained by Beaugé & Nesvorný (2013), who extended the analysis to the P – R_p diagram: this diagram shows that planets with 3 R_⊕ ≤ R_p ≤ 11R_⊕ that orbit their host star in ≲4 days are very rare. These works pointed out that this region is constrained by two boundaries: an upper boundary delimited by hot Jupiters, and a lower boundary consisting of super-Earths and sub-Neptunes. Mazeh et al. (2016) obtained an analytic description of the two boundaries delimiting the desert in the log P − log M_p and log P − log R_p planes. Although the debate about the origin of the desert is still ongoing, the distribution of close-in exoplan-ets is the result of a non-trivial combination of atmospheric and dynamical processes. The most accredited scenario today explains the lower boundary by means of the photoevaporation of the H/He atmospheres of an initially low-mass planet (Davis & Wheatley 2009; Lopez & Fortney 2014; Owen & Lai 2018), while the upper boundary is consistent with numerical simulations based on the high-eccentricity migration scenario and the subsequent tidal angular-momentum exchange of the planet with the star (Benítez-Llambay et al. 2011; Matsakos & Königl 2016; Owen & Lai 2018; Owen 2019). Photoevaporation is the physical process through which exoplanets shed their H/He envelopes mainly during the first few hundred million years, primarily due to the influence of ultraviolet (UV) and X-ray radiation from their host stars (Owen & Wu 2013). This reduces the planet size. In this scenario, Neptune-sized planets that either formed within the desert or migrated into it undergo a rapid transformation and become super-Earths. It is worth noting that the lower boundary is influenced by the cumulative UV/X-ray radiation exposure over the planetary lifetime, along with the spectral characteristics of the host star (McDonald et al. 2019). Observational constraints upon the photoevaporative processes that sculpt the desert, or its lower boundary, might come from using the helium as a tracer of escape in UV spectroscopy (Oklopčić & Hirata 2018; Allart et al. 2018; Nortmann et al. 2018; Salz et al. 2018; Allart et al. 2019). Guilluy et al. (2023) conducted a systematic He I survey of nine exoplanets that belong to the edges of the sub-Jovian and Neptune desert, but found no signature of He absorption at 3σ. Although the recent findings of Thorngren et al. 2023 would suggest that photoevaporation still partially causes the formation of the upper boundary, it cannot fully explain the upper edge of the sub-Jovian and Neptune desert as more massive planets (>0.5 M_J) withstand the effects of photoevaporation. Numerical simulations by Ionov et al. (2018) and Owen & Lai (2018) have shown that if photoevaporation were the sole mechanism that sculpts the upper limit, the P − R_p plane would be populated with sub-Jovian planets that are located at extremely short orbital periods, which contradicts the observed distribution. Observational evidence that supports these simulations was recently provided by Vissapragada et al. (2022). In their exploration of helium absorption signals from seven gas-giant planets near the upper boundary, they found that the atmosphere lifetimes of these planets significantly exceed ten billion years. An alternative explanation for the upper limit involves a high-eccentricity migration scenario, in which gravitational interactions can compel a planet into an extremely elliptical orbit that subsequently circularises as the planet draws closer to its parent star. When the orbit is circularised, planets with masses greater than 1 M_J can migrate even closer to their parent star, which is facilitated by the tidal decay caused by the stellar gravitational influence. Observational constraints upon the migration scenario are likely to come from spin-orbit angle measurements of planets lying in proximity of the sub-Jovian and Neptune desert, which is the ultimate aim of the Desert-Rim Exoplanets Atmosphere and Migration (DREAM; Bourrier et al. 2023) project. In particular, Bourrier et al. (2023) conducted a homogeneous analysis of 14 close-in exoplanets using the Rossiter-McLaughlin revolution technique and found 9 planets on misaligned orbits. In a follow-up paper, Attia et al. (2023) investigated how the spin-orbit misalignments are influenced by the star-planet interactions.

It is generally accepted that the sub-Jovian and Neptune desert cannot be the result of observational biases. Most of the Kepler candidates lying in the desert were found to be a false-positive detection (Sanchis-Ojeda et al. 2014). However, an increasing number of TESS¹ hot-Neptune planet candidates is spectroscopically confirmed as bona fide planets (e.g. Armstrong et al. 2020; Burt et al. 2020; Murgas et al. 2021; Persson et al. 2022; Naponiello et al. 2023), as is also shown through the contribution of ground-based transit surveys (e.g. Eigmüller et al. 2017; West et al. 2019; Jenkins et al. 2020; Smith et al. 2021). The desert is thus not completely empty, and some authors referred to it as an oasis (Murgas et al. 2021) or savanna (Bourrier et al. 2023; Szabó et al. 2023; Kálmán et al. 2023).

Castro-González et al. (2024) introduced the term Neptunian ridge to refer to the region of the period-radius space that separates the desert (no planets at the 3σ level) from the savanna with 3.2 days ≲ P ≲ 5.7 days. In addition, Melton et al. (2024a) applied a statistical procedure to a sample of 0.9 million light curves and identified, among other things, hundreds of new TESS candidates that lie within the sub-Jovian and Neptune desert (Melton et al. 2024b, 2023). If a large fraction of this sample were confirmed by future follow-up efforts, then it would indicate that the desert is not as barren as it appears, but might be caused by an unknown bias in the Kepler survey. Nonetheless, the recent vetting efforts of Magliano et al. (2023a,b) based on TESS candidates showed that the false-positive occurrence rate is ≈75% for candidates in the sub-Jovian and Neptune desert and ≈34% outside of it. This demonstrates that hot-Neptune candidate planets are twice as likely to be false-positive detections.

In this tangle of conflicting results, a different perspective in which to contextualise this issue might deliver new and interesting cues. Lundkvist et al. (2016) were the first to study the hot-super-Earth desert issue and considered the incident flux on the planetary surface rather than the orbital period. The orbital period is not a completely informative variable for defining a planet as ‘hot’ because it does not take the characteristics of the host star into account. Szabó & Kálmán (2019) and Szabó et al. (2023) showed that the boundary of the desert depends on the stellar effective temperature, metallicity, and log 𝑔, in order of significance. They suggested that the dependence on the fundamental stellar parameters may be an indication of the predominance of photoevaporative processes in determining the boundary, while the non-dependence on tidal forces and log 𝑔_p of the planet is consistent with dynamical processes that partly cause the desert. Thus, the term ‘hot’ ought to also be correlated to the properties of the star in order to capture the full picture of these systems. Pioneering works that first addressed the sub-Jovian and Neptune desert conundrum (e.g. Mazeh et al. 2016) employed the orbital period as the key parameter since it is the most affordable observable. Through the astronomical efforts of different missions (e.g. Gaia; Gaia Collaboration 2016), a preliminary characterisation of the parent star is almost always feasible today. For this reason, we study the issue of the sub-Jovian and Neptune desert by replacing the orbital period with the incident flux on the planetary surface.

The outline of the paper is as follows. In Sect. 2 we describe our dataset and its features in the F − R_p and F − M_p diagrams. In Sect. 3 we then retrieve the analytic expressions of the two boundaries in the two diagrams within a Bayesian framework. Finally, we summarise our conclusions in Sect. 4.

2 Methods

The insulation flux F corresponds to the amount of stellar energy received per unit area on the planetary surface and is defined as $$F=\frac{L_{*}}{4\pi d^2},$$

where L_* is the luminosity of the star, and d is the distance of the planet from the star (semi-major axis of the orbit). For context, Earth’s insulation flux F_⊕ is approximately 1361 W/m², and it provides a useful reference point for comparing the energy received by other planets. However, the insulation flux F that a planet receives from its host star (modelled as a black-body source) in units of the Earth value F_⊕ in the circular orbit approximation can also be expressed as $$\frac{F}{F_{\oplus }}={\left(\frac{{\rho }_{*}}{{\rho }_{\odot }}\right)}^{-2/3}{\left(\frac{P}{1\text{\hspace{0.17em}yr}}\right)}^{-4/3}{\left(\frac{T_{*}}{T_{\odot }}\right)}^4,$$(1)

where ρ_*. is the stellar mean density, ρ_⊙ = 1.408 g/cm³ is the mean density of the Sun, P is the orbital period of the planet, T_* is the stellar effective temperature, and T_⊙ = 5778 K is the effective temperature of the Sun. Thus, the insulation flux simultaneously takes the orbital period of the planet and the host star temperature and mean density into account. At a fixed orbital period, the radiation received by the planet grows when it orbits a hotter and less dense star.

We initially used a dataset that consisted of 3335 confirmed exoplanets² with a known orbital period and planetary radius that orbit a star with an estimated mass and radius. The corresponding stellar parameters were retrieved from the NASA Exoplanet Archive. The requirement of a known planetary radius implies that approximately 99% of the planets in the sample were discovered using the transit method. This observational bias likely results in an overrepresentation of larger, short-period exoplanets, particularly those orbiting bright stars. While this bias does not significantly affect the primary findings of this work, such as the observed paucity of highly irradiated Neptune-and sub-Jovian-sized planets, it may influence the accuracy of the boundary calculations at lower flux values, where smaller or more distant planets are potentially underrepresented. Future exoplanetary missions that aim to expand the demographic reach of exoplanet detection will be crucial in strengthening and refining the accuracy of these results.

We removed 1808 (≈54% of the total) exoplanets from this initial sample for one of the following reasons: (i) Their radius or period estimates had an uncertainty larger than 10%, or (ii) they orbit stars whose radius or mass estimates have an uncertainty larger than 10%. After this filtering, our final sample (S₁ hereafter) was composed by 1527 exoplanets. We used a threshold value of 10% for the maximum uncertainty in order to minimise the propagated error over the insulation flux and to have the cleanest and most complete dataset possible. Only 0.6% of the planets in S₁ have ∆F/F > 40%, and none exceeds 50%.

2.1 Orbital period against incident flux

In order to investigate whether F is a more appropriate physical quantity to describe the sub-Jovian and Neptune desert, we preliminary divided the S₁ sample into four groups of exoplanets: (i) the first group was composed of 586 planets that received an insulation flux F < 50 F_⊕, (ii) the second group consisted of 364 planets whose incident flux was between 50 F_⊕ and 200 F_⊕, (iii) 238 planets in the third group had 200 F_⊕ < F ≤ 550 F_⊕, and iv) the fourth group comprised the remaining 339 planets with F > 550 F_⊕. We emphasise that the flux binning we used in this step was completely arbitrary and was simply provided as a guideline to distinguish between extremely and moderately irradiated objects. The selection of four bins aimed to balance the need for granularity in the analysis with the necessity of maintaining sufficient statistical power within each bin. We ensured that each bin contained a sufficient number of data points that allowed us to apply robust statistical methods. In Fig. 1 we depict the distribution of the four groups in the P – R_p plane. This chart would suggest that the less irradiated exoplanets are more uniformly distributed than the highly irradiated ones. At F < 50 F_⊕, there is no evidence of a sub-Jovian and Neptune desert. However, as F increases, a gap opens in the distribution, and the Neptune/sub-Jupiter region in R_p is voided of planets. The lower left and right panels (200 F_⊕ < F ≤ 550 F_⊕ and F > 550 F_⊕, respectively) show a clear dichotomy of highly irradiated exoplanets: hot Jupiters or hot super-Earths/sub-Neptunes. It therefore appears that the insulation flux plays the main role in sculpting the sub-Jovian and Neptune desert, while the orbital period is an incomplete physical parameter, as we show in Sect. 3. We also point out the well-known lack of inflated Jupiters (R_p ≳ 15 R_J) that receive less than 550 F_⊕ (Laughlin et al. 2011; Miller & Fortney 2011; Demory & Seager 2011; Laughlin & Lissauer 2015; Thorngren & Fortney 2018). These planets, similar to Jupiter in composition, appear to be puffed up and less dense. The cause of this inflation is not fully understood, but may involve factors such as stellar irradiation, internal heat, and tidal forces. In order to statistically demonstrate the differences among the four aforementioned planetary populations, we employed a two-sample Kolmogorov-Smirnov test (e.g. Feigelson & Babu 2012) on the planetary radius distributions of these groups. Specifically, for each pairing of the four datasets, we derived a p value lower than 10⁻¹⁵, which conclusively establishes the statistical disparity in the planetary radius distributions in all four groups.

We concluded this preliminary analysis of our dataset by replicating the steps performed by Szabó & Kálmán (2019) on our S₁ sample, which is nearly three times larger than the sample they used, despite the stricter constraints we adopted on the uncertainties over each parameter. Szabó & Kálmán (2019) employed the two-sample Kolmogorov-Smirnov test to assess differences in the planetary distributions in the P – R_p diagram according to a third variable (e.g. stellar mass, radius, or metal-licity). By comparing the characteristics of exoplanets inside and outside the desert region, the study shed light on the factors that contribute to its formation. We confirmed the dependence of the desert edges in the P – R_p diagram on the stellar effective temperature, mass, and metallicity found by Szabó & Kálmán (2019). In particular, we recovered that planets in the desert are likely to have a shorter period around cooler, less massive stars and more metal-rich stars, as obtained by Szabó & Kálmán (2019). However, the role of the metallicity in sculpting the desert was statistically downsized by the latest work of Szabó et al. (2023) based on a new dataset in which the stellar parameters were uniformly retrieved in the Apache Point Observatory Galactic Evolution Experiment (APOGEE) DR17 sample (Majewski et al. 2017; Zasowski et al. 2017).

2.2 The sub-Jovian/Neptunian desert in the flux-radius plane

In the left panel of Fig. 2, we display the distribution of the 1527 planets in S₁ in the F_⊕ R_p plane. This parameter space shows a scarcity of Neptune- and sub-Jupiter-like planets that receive at least ≈200 F_⊕ and form a vaguely triangular region that is delimited by two boundaries, as delineated in grey in Fig. 2 based on visual inspection. The lower edge comprises highly irradiated super-Earths and sub-Neptunes, and the upper edge consists of hot Jupiters. Furthermore, our sample contains 221 planets that lie between the boundaries of the desert obtained by Mazeh et al. (2016) in the P – R_p plane. These objects would be referred to as hot Neptunes if we only considered their distribution in the P – R_p plane, but most of them are found in crowded regions in the F – R_p plane, thus leaving the desert (see Sect. 3.1).

As shown in the right panel of Fig. 2, planets with 3 R_⊕ ≤ R_p ≤ 11 R_⊕ are not so rare when the flux received from their host star is <200 F_⊕. In particular, Neptunes and sub-Jupiters that experience an incident flux F ≤ 200 F_⊕ are three times more common than those with F > 200 F_⊕. In order to describe the sub-Jovian and Neptune desert in this new parameter space, we studied the radius distribution of the planets in S₁ in terms of the received insulation flux. The left panel of Fig. 2 suggests that there are two distinct populations of exoplanets: (i) sub-Neptunes and super-Earths with an incident flux spanning almost five orders of magnitudes, and (ii) Jupiters and inflated hot Jupiters that are found to withstand more stellar radiation on average. We applied a Gaussian mixture model (GMM; Ivezić et al. 2020) to distinguish between the two classes of planets. A GMM is a probabilistic model that is used for clustering and a density estimation. It assumes that the data points are generated from a mixture of several Gaussian distributions, each characterised by its mean and covariance. For the multivariate case, a GMM with K components is described by a density probability p(x) given by $$p(x)={\displaystyle \sum _{i=1}^K{\varphi }_i}\mathcal{N}\left(x\mid {\mu }_i,{\text{Σ}}_i\right),$$(2)

where µ_i and Σ_i correspond to the mean and the covariance of the ith component, respectively, 𝒩 refers to a Gaussian density probability, and ϕ_i is the mixture component weight of the ith component, such that ${\displaystyle {\sum }_{i=1}^K{\varphi }_i}=1$. We calculated the best parameters for each component by using the expectation-maximisation algorithm (Meng & Van Dyk 2002). The expectation stage calculates the probability of each data point belonging to each Gaussian component, and the maximisation step fine-tunes the model parameters based on these probabilities. We performed a two-component (i.e. K = 2) bivariate GMM in the log F -log R_p space in order to cluster Si into two groups, as shown in Fig. 3. We opted to work with the log values since the planetary fluxes are spread over a wide range of magnitudes. Figure 3 clearly shows two well-separated class of exoplanets. The first component of the model, group A_R hereafter, consists of a variety of exoplanet types: super-Earths/sub-Neptunes, hot Neptunes/Saturns, and cold giants. We stress that there is still a less dense (but not empty) populated region of highly irradiated sub-Jovians and Neptunes that would represent the savanna, similar to what we found in the P − R_p and P – M_p diagrams. We found the mean of group A_R to be located at ${\mathit{\mu }}_{A_R}\simeq \left(50F_{\oplus },2R_{\oplus }\right)$. On the other hand, the second component corresponds to hot Jupiters and inflated hot Jupiters, group B_R hereafter. The mean of group B_R corresponds to ${\mathit{\mu }}_{B_R}\simeq \left(523F_{\oplus },13R_{\oplus }\right)$. The GMM allowed us to label each point of our diagram as belonging to one of the two groups. Since the focus of this work is on highly irradiated exoplanets, we only considered the bodies that received at least 200 F_⊕. Thus, we took the 75th percentile of the incident flux distribution of planets with 3 R_⊕ ≤ R_p ≤ 11 R_⊕ as the threshold value above which the sub-Jovian and Neptune desert is defined.

The next step was to determine the borders in the log F − log R_p diagram that define the lower and upper edges of the desert. In particular, we defined the edges of the two groups with a quantitative criterion, and the results are shown in Fig. 4. We sliced the diagram in terms of the incident flux in order to study the planetary radius distribution of groups AR and BR in each slice. For both groups, we defined a custom grid of F values using a binning that ensured at least 20 elements per group. This approach allowed us to calculate statistically reliable uncertainties for the representative value of each group. In particular, the main purpose of the grid-based analysis was to calculate the points that were used below to derive the borders (see Sect. 3.1). Group AR is slightly larger than group BR, and its grid had ten slices, while the grid of group BR had seven slices. The two grids span [200 F_⊕, 4000 F_⊕], which provides enough points for a statistical analysis. Moreover, as shown in Fig. 3, the Rp distribution of group AR is skewed towards larger radii by the cold and warm sub-Jovians/Jupiters. We also applied a Shapiro-Wilk test (Shapiro & Wilk 1965) to verify that it is not normally distributed. We rejected the null hypothesis with p < 10⁻¹³. In particular, the scarcity of warm and hot sub-Jovians belonging to group A_R, within the savanna and nearby the exterior ellipsoids of group B_R, would still contaminate the lower border calculation. In order to avoid contamination from hot sub-Jovians of group AR inside the desert, we removed them by performing a 2σ clipping for the R_p distribution of group A_R in each slice.

We defined the lower boundary point in the ith slice of the flux grid related to group A_R as $$R_{\text{low},i}={\mu }_{A_R,i}+2{\sigma }_{A_R,i}i=1,\mathrm{..},10,$$(3)

where ${\mu }_{A_R,i}$ and ${\sigma }_{A_R,i}$ represent the mean and the standard deviation of the the Rp distributions of group AR in the ith slice, respectively. The incident flux F_low,i, associated with R_low,i, was taken as the average flux within the ith slice. The uncertainties ∆R_low,i and ∆F_low,i· on the mean planet radius and incident flux, respectively, were estimated as the error on the standard deviation of the Rp distribution for group A_R and the standard deviation of the F distribution for group A_R, both evaluated for the ith slice. We note that if we had not employed the sigma-clipping process, ${\sigma }_{A_R,i}$ in each slice would have been affected by the presence of a few hot sub-Jovians and would thus have yielded a point (F_low,i, R_low,i) that would have fallen inside the desert instead of at its lower edge. This methodological difficulty is caused by the Neptunian ridge (Castro-González et al. 2024) in the flux-radius diagram that separates the desert from the underpopulated savanna.

We applied an analogous treatment for the planets falling within group B_R in order to build up the upper boundary of the desert. We did not employ any sigma clipping since no planet of group B_R contaminated the upper edge calculation. In this case, the upper boundary in planet radius is defined as $$R_{\text{up},j}={\mu }_{B_R,j}-2{\sigma }_{B_R,j}j=1,\dots 7.$$(4)

The decision to define the boundary points as the points with a z-score equal to 2 in each distribution is arbitrary, but was made in light of an a posteriori visual analysis. As shown by the confidence ellipsoids in Fig. 3, if we used z = 3 or z = 1, the points would be significantly farther out from the edges on the outside and inside, respectively. In particular, in each slice, more than 90% of exoplanets of groups A_R and B_R are found to be above (below) the upper (lower) boundary. Moreover, although at least 20 objects in each bin might be an arbitrary choice, this was the best trade-off between having solid statistical results and a sufficient number of boundary points to later fit with high-Dlmensional And multi-MOdal NesteD Sampling (DIAMONDS, Corsaro & De Ridder 2014). When we decreased the number of objects per bin to 10 in order to almost duplicate the number of boundary points, those describing the lower edge were mostly unreliable. This is mainly due to contaminants in the savanna whose weights in the boundary point computation are not attenuated by the poor statistics.

We observed that the R_p distribution of the super-Earths and sub-Neptunes (group A_R) slowly narrows and shifts towards a smaller mean planetary radius as F increases. On the other hand, the R_pdistribution of the giant planets (group B_R) broadens and quickly tends to a larger mean planetary radius when F increases. These two effects at high fluxes lead to the observed dichotomy of highly irradiated exoplanets and to the peculiar triangular shape of the sub-Jovian and Neptune desert. In Fig. 4 we show the points that delimit the lower and upper edges in (F, Rp) space within the range [200 F_⊕, 4000 F_⊕].

Fig. 1 Distribution of the confirmed exoplanets of the S1 sample in the log P – log Rp diagram for different ranges of the insulation flux. The upper left panel shows the 586 planets that receive ≤50 F⊕, the upper right panel shows the 364 planets that receive between 50 F⊕ and 200 F⊕, the lower left panel depicts the 238 planets that receive between 200 F⊕ and 550 F⊕, and the lower right panel illustrates the 339 planets with F > 550 F⊕.

Fig. 2 Analysis of the planetary distribution in the F – Rp plane and the flux cumulative distribution function. Left panel: distribution of the 1527 planets of the S1 set in the F – Rp plane. The yellow triangles represent planets that would be referred to as hot Neptunes according to the limits obtained by Mazeh et al. (2016). The grey triangle is a qualitative representation of the desert. Right panel: cumulative distribution function of incident flux values F of planets with 3 R⊕ ≤ Rp ≤ 11 R⊕. The vertical dashed line in both panels represents the threshold value of 200 F⊕ below which lie 75% of the planets with 3 R⊕ ≤ Rp ≤ 11 R⊕.

Fig. 3 Application of a two-component bivariate GMM over S1 sample in the log F – log Rp space. We found two separate populations: group AR (purple dots) with a mean value of (1.70,0.34), and Group BR (yellow dots) with a mean value of (2.72, 1.12).

Fig. 4 Calculated boundaries in (F, Rp) space: (Flow,i, Rlow,i) pairs (purple points) sample the lower boundary, and (Fup,i, Rup,j) pairs (red points) mark the upper boundary.

2.3 The sub-Jovian/Neptunian desert in the flux-mass plane

Within our initial sample, 609 planets out of 1527 also have a planetary mass estimate with ∆M_p/M_p < 30%. We refer to this set as sample S 2. We stress that the higher threshold value for ∆M_p /M_p with respect to ∆R_p /R_p was set to guarantee enough statistics in the grid-analysis step. If we had required a ∆M_p/M_p < 10%, it would have halved the S₂ sample. In Fig. 5 we show the distribution of the objects belonging to the S 2 sample in the F − M_p plane. In this parameter space, we still observe a paucity of Neptunes and sub-Jupiters that receive at least ≈550 F_⊕. In this case as well, only a few Hot Neptunes, defined as the objects within the boundaries of Mazeh et al. (2016), lie within the underpopulated region (see Sect. 3.2). The right panel of Fig. 5 shows that planets with 10 M_⊕ ≤ M_p ≤ 200 M_⊕ are common when the flux received from their host star is ≲550 F_⊕. Neptunes and sub-Jupiters that experience an incident flux F ≤ 550F_⊕ are three times more frequent than those with F > 550 F_⊕. Therefore, in contrast to the flux-radius diagram, the vertex of the triangle (namely the 75th percentile of the incident flux distribution of planets with mass 10 M_⊕ ≤ M_p ≤ 200 M_⊕) here appears to be shifted towards higher incident flux (see Sect. 2.4). As done in Sect. 2.1, we confirmed all the results found by Szabó & Kálmán (2019) and Szabó et al. (2023) in the P − M_p plane, with planets in the desert likely to have a shorter period around cooler and more metal-rich stars. We also recovered their finding that the stellar mass does not play a critical role in sculpting the desert in the P − M_p diagram.

We repeated the same steps as in Sect. 2.2 considering the planetary mass in order to find the analytic expression of the two boundaries in the (F, M_p) space. We report the result of the two-component bivariate GMM in log F − log M_p space in Fig. 6. Group AM consists of a wide variety of exoplanet types spanning several orders of magnitude in incident flux (from ≈0.1 F_⊕ to ≈10⁴ F_⊕) as well as a broad spectrum of masses (from ≈0.1 M_⊕ to ≈8000 M_⊕). We also point out that the most massive objects in group AM tend to experience low levels of incident flux. We found the mean of group A_M to be located at ${\mu }_{A_M}\simeq \left(78F_{\oplus },27M_{\oplus }\right)$. On the other hand, group B_M corresponds to planets that are more massive and more irradiated than those in group A_M on average. We found the mean of group B_M to be located at ${\mu }_{B_M}\simeq \left(939F_{\oplus },347M_{\oplus }\right)$. It also covers a narrower flux and mass ranges, as shown by the confidence bands of the ellipsoids. However, group B_M still spans two orders of magnitudes in planetary mass beyond the threshold flux of 550 F_⊕, and it is not Gaussian (Shapiro-Wilk test, p < 10⁻⁶) but a skewed distribution towards more massive planets, as shown by the marginal histogram in Fig. 6.

Thus, we needed to remove extremely massive exoplanets that otherwise would lead to large spreads in the giant populations and would distort the upper boundary computation. For this reason, we discarded all the planets with a planetary mass M_p > 600 M_⊕, thus resulting in a final sample composed of 512 exoplanets with a planetary mass estimate (S₃ sample hereafter). The choice of using a cutoff value of 600 M_⊕ was arbitrary, but is capable of guaranteeing enough statistics in each slice.

This operation does not affect the robustness of our results since we are interested in describing the low-mass edge of the group B_M distribution. However, the S₃ sample contains fewer than half of the objects that were contained in the initial sample we used in Sect. 2.2, so that the statistical analysis yields larger errors over the upper and lower boundary points. Similarly, the fitted parameters through DIAMONDS are more uncertain than those retrieved in Sect. 2.2. A summary of the three samples defined in this work is given in Table 1

For the boundary calculations, we only considered exoplanets that received at lest 550 F_⊕. For both groups in our sample, we created a custom grid of F values spaced such that each slice contained at least 20 elements per group. As anticipated, since the sample is much less numerous in this space of parameter, we needed to choose larger slices than we defined in Sect. 2.2. For this reason, the grid used to sample group A_M consisted of only three slices and that used to slice group B_M contained just four slices. The two grids span the range [550 F_⊕, 4000 F_⊕], and the upper limit was set to ensure enough points for a statistical analysis. In this case as well, the contamination of highly irradiated sub-Jovian planets within the savanna forced us to perform a 2σ clipping for the M distribution of group A_M in each slice.

In contrast with the boundary definitions given in Eqs. (3)– (4), since the planetary masses were spread over several orders of magnitudes, we defined the boundary points as points with a z-score equal to one, as suggested by the confidence ellipsoids in Fig. 6. In this case as well, the choice of the proper ɀ score was based on an a posteriori visual inspection of the points delimiting the two boundaries. However, due to the poorer statistics in each slice, more than 75% of the exoplanets in groups A_M and B_M were found to be above (or below) the upper (or lower) boundary.

In particular, the lower boundary point in the ith slice of the flux grid related to group A_M was calculated as $$M_{\text{low },i}={\mu }_{A_M,i}+{\sigma }_{A_M,i}i=1,\mathrm{..},3.$$(5)

The upper boundary point in the jth slice of the flux grid related to group B_M is given by $$M_{\text{up},j}={\mu }_{B_M,j}-{\sigma }_{B_M,j}j=1,\dots ,4.$$(6)

In Fig. 7 we show the points that delimit the lower and upper boundary in the (F, M_p) space within the range [550 F_⊕, 4000 F_⊕].

Fig. 5 Analysis of the planetary distribution in the F – Mp plane and the flux cumulative distribution function. Left panel: distribution of the 609 planets of the S2 set in the F – Mp plane. The yellow triangles represent planets that would be referred to as hot Neptunes according to the limits obtained by Mazeh et al. (2016). The grey triangle is a qualitative representation of the sub-Jovian and Neptune desert. Right panel: cumulative distribution function of incident flux values F of planets with 10 M⊕ ≤ Mp ≤ 200 M⊕ . The vertical dashed line represents the threshold value of 550 F⊕ in both panels below which 75% of the planets with 10M⊕ ≤ Mp ≤ 200 M⊕ lie.

Fig. 6 Application of a two-component bivariate GMM on the sample S2 in the log F – log Mp space. We identified two separate populations: group AM (purple dots) with a mean value of (1.8,1.4), and group BM (yellow dots) with a mean value of (2.9, 2.5).

Fig. 7 Calculated boundaries in the (F, Mp) space: (Flow,i, Mlow,i) pairs (purple points) sample the lower boundary, and (Fup,i, Mup, j) pairs (red points) mark the upper boundary of the desert.

2.4 About the threshold incident flux

In the F – R_p and F – M_p planes (see Sects. 2.2 and 2.3), we defined the threshold incident flux from which to start our analysis as the 75th percentile of the distribution of exoplanets with 3 R⊕ ≤ R_p ≤ 11 R_⊕ and 10 M_⊕ ≤ M_p ≤ 200 M_⊕, respectively. We refer to this value as F_75th. We found F_75th ≈ 200 F_⊕ in the F – R_p plane and F_75th ≈ 550 F_⊕ in the F – M_p plane. This raises the intriguing question whether the variation in F_75th between F – R_p and F – M_p reflects an underlying physical cause, or if is it merely a consequence of the sparser data in the mass plane. In order to address the second option, we employed two different strategies: (i) We reduced the size of the S₁ sample to match that of S₂ and calculated F_75th for planets with 3R_⊕ ≤ R_p ≤ 11 R_⊕, and (ii) we enlarged the size of the S 2 sample to match that of S₁ and calculated F_75th for planets with 10 M_⊕ ≤ M_p ≤ 200 M_⊕.

2.4.1 Downsizing the S₁ sample

We performed 10 000 Monte Carlo simulations, in each of which the S₁ sample (comprising 1527 planets) was reduced to the size of the S₂ sample (609 planets) by randomly removing 918 objects. In each simulation, we calculated the F_75th of exoplanets with 3 R_⊕ ≤ R_p ≤ 11 R_⊕ , whose distribution is shown in Fig. 8. We obtained a non-Gaussian F_75th distribution with a mean value of ≈163 F_⊕ and a standard deviation of ≈25 F_⊕. The distribution shows a slight right-handed tail that is consistent with a skewness of ≈0.4. The non-Gaussianity was also investigated using the Anderson-Darling (AD) test (Anderson & Darling 1952), whose statistics of 65 exceeds the highest critical value by far (1.092) at the 1% significance level. The Gaussianity null hypothesis is therefore rejected at all of the provided significance levels. We opted for the AD test instead of the more commonly used Shapiro-Wilk test since the latter is known to perform better for small sample sizes (Shapiro & Wilk 1965). We found that 99.84% of our simulations yielded F_75th values lower than 250 F⊕ . Moreover, only one simulations out of 10 000 returned a value of F_75th ≳ 300 F_⊕ , which is still much lower than the 550 F⊕ we obtained for the F – M_p diagram.

Fig. 8 Histogram of the variable F75th for planets with 3 R⊕ ≤ Rp ≤ 11 R⊕, generated from 10 000 Monte Carlo simulations. Each simulation reduced the S1 sample to the size of the S2 and calculated the 75th percentile of exoplanets with 3 R⊕ ≤ Rp ≤ 11 R⊕ . The histogram is normalised to show the probability density.

2.4.2 Augmenting the S₂ sample

In order to understand the effect of the small sample size on the definition of F_75th and ultimately on the borders of the desert, we considered an augmented version of S₂ by assigning an approximate mass estimate, derived using forecaster (Chen & Kipping 2016), to all objects in the S₁ sample without a measured planetary masses. We also included in the new enlarged sample planets with ∆M_p /M_p > 30%. This procedure resulted in a S₂ sample of 1527 confirmed planets, each with an estimated mass, regardless of the magnitude of the associated uncertainty. Nevertheless, even in this scenario, the F_75th for planets with 10M_⊕ ≤ M_p ≤ 200 M_⊕ remains slightly above 400F_⊕.

Thus, our statistical analysis suggests that the threshold value of the incident flux F_75th is not affected by the poorer statistics in the flux-mass plane, but might be the result of a non-trivial physical effect. A similar result was obtained by Ma & Ghosh (2021): by employing a Bayesian mixture model, they found that the incident flux has a non-negligible effect on the M_p – R_p relation. Future multivariate approaches that simultaneously take into account (F, R_p, M_p) might shed light upon this issue.

2.5 The mass-radius relation

We also checked whether the GMM clustering processes performed in the two different spaces of parameters yielded any inconsistent results. In particular, one might wonder whether a planet of group A_R (small radius) in the F– R_p diagram still belongs to group A_M (low mass) in the F – M_p space or if it falls within the domain of group B_M (high mass). If this were not the case, it would mean that planets exist whose classification depends on the physical parameter used to describe them. We performed this analysis using the S₃ set of 512 exoplanets with a well-constrained planetary mass and planetary radius measurements. Our results are listed below.

1. We found that 271 planets (≈53% of S₃ sample) are classified as belonging to group A according to both radius and mass (i.e. small radius A_R and low mass A_M). We term this class A_RA_M hereafter.

2. We found that 29 planets are classified as belonging to group B with respect to the radius and to group A with respect to the mass (i.e. large radius B_R and low mass A_M). We term this class B_RA_M hereafter.

3. We found that one planet (K2-60 b; Eigmüller et al. 2017) was classified as belonging to group A with respect to the radius and to group B with respect to the mass (i.e. small radius A_R and high mass B_M). We term this class A_RB_M hereafter.

4. We found that 211 planets (≈41% of S₃ sample) were classified as belonging to group B for the radius and mass (i.e. large radius B_R and high mass B_M). We term this class B_RB_M hereafter.

In Fig. 9 we show the distribution of 512 objects within the log R_p – log M_p plane, each coloured according to its classification with respect to both radius and mass. It is interesting to note that the 29 planets classified as B_RA_M belong to the transition of the M(R) curve. The average mean planetary density of the B_RA_M group is ≈0.68 g/cm³ (approximately half of the mean Jupiter density) with a standard deviation of ≈0.61 g/cm³ . Thus, these objects might represent the well-known class of inflated hot Jupiters (Sestovic et al. 2018). We remark that some of the hybrid objects (i.e. belonging to groups A_RB_M and B_RA_M) might simply be the result of the precision limit of the GMM clustering. K2-60 b (Eigmüller et al. 2017) has a mean planetary density of ≈1.6g/cm³. Future spectroscopic follow-up investigations of these objects could shed light upon their proper position within the log R_p – log M_p plane.

2.6 Bayesian inference and model comparison with DIAMONDS

In order to fit the boundaries in the F – R_p plane that delimit the region of the desert, we exploited a series of different analytical models whose statistical significance was subsequently tested by means of a Bayesian model comparison process. All the fits were performed using the public Bayesian inference tool DIAMONDS³ (Corsaro & De Ridder 2014). DIAMONDS allowed us to simultaneously perform a Bayesian parameter estimation and model comparison using the nested-sampling Monte Carlo algorithm (Skilling 2004), without limitations regarding the size of the dataset or the number of free parameters involved in the model. DIAMONDS returned the Bayesian evidence (E) for each tested model. The Bayesian evidence represents the probability of the data in the light of a given model, that is, it averages the likelihood distribution over the parameter space set by the priors, namely by our choice of values of the free model parameters. In particular, E is defined as $$ℰ\equiv {\displaystyle {\int }_{{\text{Σ}}_{ℳ}}ℒ}(\theta )\pi (\theta \mid ℳ)\text{d}\theta ,$$(7)

where θ is the parameter vector, ℒ is the likelihood, and π is the prior on the free parameters as conditional to the choice of the model ℳ, the latter defining the parameter space Σ_ℳ over which the integral is performed. When we compare the two models M1 and ℳ₂, the corresponding Bayes factor B_1,2 of ℳ₁ relative to ℳ₂ is defined as the ratio of the model evidence, or equivalently, as ln B_1,2 ≡ ln ℰ₁ – ln ℰ₂ when we adopt a logarithmic form. For ln B_1,2 ≥ 5 (Jeffreys 1961), the evidence condition in favour of model ℳ₁ is strong, meaning that this model can be considered with a probability higher than 99% to be the favoured choice.

The first model we took into account was a power-law model of the form $$\frac{R_p}{R_{\oplus }}=\beta {\left(\frac{F}{F_{\oplus }}\right)}^{\alpha },$$(8)

with the free parameters α and β. In order to incorporate the uncertainties on the independent variable F, we linearised this model by adopting a logarithmic scale, $$\ln \left(\frac{R_p}{R_{\oplus }}\right)=\alpha \ln \left(\frac{F}{F_{\oplus }}\right)+\ln \beta .$$(9)

This model, hereafter model ${ℳ}_1^R$, has the free parameters (α, ln β).

As a second option, we considered the linear model $$\frac{R_p}{R_{\oplus }}=a\left(\frac{F}{F_{\oplus }}\right)+b,$$(10)

hereafter termed model ${ℳ}_2^R$, with (a, b) as free parameters.

We performed a Bayesian inference on the two models presented above by adopting a Gaussian likelihood that incorporated the uncertainties on the planet radii and fluxes. The corresponding log-likelihood reads $$\text{Λ}(\theta )={\text{Λ}}_0(\theta )-\frac{1}{2}{\displaystyle \sum _{i=1}^N{\left[\frac{{\text{Δ}}_i(\theta )}{{\tilde{\sigma }}_i(\theta )}\right]}^2},$$(11)

where θ is the parameters vector (e.g. (α, ln β) for model ℳ₁), N is the total number of points considered in each boundary, ${\tilde{\sigma }}_i(\mathit{\theta })$ is the relative uncertainty over θ, and Λ₀(θ) is a term depending on the relative uncertainties, given by $${\text{Λ}}_0(\theta )=-{\displaystyle \sum _{i=1}^N\text{\hspace{0.17em}ln\hspace{0.17em}}}\sqrt{2\pi }{\tilde{\sigma }}_i(\theta ).$$(12)

Lastly, we considered a third model, termed ${ℳ}_3^R$, which extended the linearity to a quadratic regime, but this time, without the uncertainties in the covariates F. The model reads $$\frac{R_p}{R_{\oplus }}=\tilde{a}\left(\frac{F}{F_{\oplus }}\right)+\tilde{b}{\left(\frac{F}{F_{\oplus }}\right)}^2+\tilde{c}$$(13)

and the free parameters are $(\tilde{a},\tilde{b},\tilde{c})$. The outputs of DIAMONDS in the flux-radius diagram are discussed in Sect. 3.1.

In addition to the investigation of the F – R_p plane, we performed a similar analysis in the F – M_p space. In particular, we analysed the analogous of models ${ℳ}_1^R$ and ${ℳ}_2^R$ by replacing the planet radius R_p with the planet mass M_p, from which we obtained models ${ℳ}_1^M$ and ${ℳ}_2^M$. We did not consider the quadratic model ${ℳ}_3^R$ because too few points are available to constrain the larger number of free parameters. The outputs of DIAMONDS in the flux-mass diagram are discussed in Sect. 3.2.

Fig. 9 Distribution of 512 objects of sample S3 within the log Rp – log Mp plane, coloured according to their classification with respect to the radius and mass. 271 blue dots (ARAM) represent the planets classified as not giants for the radius and mass, 211 red dots (BRBM) represent the planets labelled as giants for the radius and mass, 29 green points (BRAM) are giant planets according to their radius and not giants according to their mass, and the one yellow point (ARBM) corresponds to a not-giant in terms of the radius and to a giant with respect to the planetary mass.

3 Results

The analysis presented in Sect. 2 shows that a paucity of highly irradiated Neptunes and sub-Jovians exists in the flux-radius and flux-mass diagrams. This feature suggests a dichotomy of planets at a high level of incident flux (≳200 F_⊕): sub- Neptunes/super-Earths, and Jupiters/super-Jupiters. By using the GMM algorithm, we determined the upper and lower edges in the flux-radius and flux-mass planes as described in Sects. 2.2 and 2.3. We then fitted each edge in the two planes using DIAMONDS, as shown in Sect. 2.6. Thus, the depletion region delimited by the two boundaries here takes the orbital period of the planet and the fundamental properties of the host star into account.

3.1 Flux-radius plane

The Bayesian evidence E computed through DIAMONDS is listed in Table 2. The model that is clearly favoured for the lower and upper boundary is the power-law model ${ℳ}_1^R$, for which the Bayes factors largely exceeds the strong-evidence condition (ln B > 5.0). We noted that the model involving a quadratic dependence upon F is subject to significant departures from the data regime for values outside the range imposed by the data. As appears evident from the reported Bayesian evidence, the larger number of free parameters of this latter model significantly penalises its statistical weight as compared to the other models that we took into account in our study (Jeffreys 1961).

In Table 3 we report the best-fitting values for α and ln β of the power-law model ${ℳ}_1^R$ as retrieved by DIAMONDS. We name the area comprised between the boundaries obtained with the ${ℳ}_1^R$ model irradiation desert in the F – R_p plane.

In Fig. 10 we show the distribution of S₁ sample in the P–R_p diagram alongside the Mazeh et al. (2016) boundaries and within the F – R_p diagram alongside the lower and upper boundaries obtained in this work. Out of 221 hot Neptunes, identified in the P – R_p diagram, 194 fall outside the irradiation desert in the F – R_p plane or lie along its edges. Future refinements of the orbital and physical parameters of the objects lying on the edges could shed light upon their membership of the desert. Only 27 out of 221 planets occupy the statistical scarcely populated region within the F – R_p diagram. Thus, according to this new perspective, the 194 planets are not so rare, even though they were earlier classified as hot Neptunes. Furthermore, we also checked whether any planet of our sample that resides outside the desert as defined by Mazeh et al. (2016) (in the P – R_p diagram) would fall inside the desert derived in this work. Only three planets satisfy this condition, and they are Kepler-1270 b ($R_p={3.32}_{-0.65}^{+0.30}{R}_{\oplus }$ and P = 6.03356196 ± 0.00006551 days; Morton et al. 2016), Kepler-815 b ($R_p={4.11}_{-0.33}^{+0.19}{R}_{\oplus }$ and P = 8.57503552 ± 0.0000591 days) and KOI-94 c (R_p = 4.32 ± 0.41 R_⊕ and P = 10.423707 ± 0.000026 days; Weiss et al. 2013; Masuda et al. 2013). As shown in Fig. 10, while these planets lie nearby but outside of the lower boundary as defined by Mazeh et al. (2016), they fall well inside the desert found in this work.

3.2 Flux-Mass plane

We also fitted the two edges of the sub-Jovian and Neptune desert within the F – M_p plane using DIAMONDS. As in the previous case, the model that is clearly favoured for the lower and upper boundary is the power-law model ${ℳ}_1^R$, and its Bayes factors largely exceeds the strong-evidence condition (ln B > 5.0), as shown in Table 2.

In Table 3 we report the best-fitting values for α and ln β of the power-law model ${ℳ}_1^M$ as retrieved by DIAMONDS. We name the area between the boundaries obtained with the ${ℳ}_1^M$ model the irradiation desert in the F – M_p plane. With respect to the results obtained for the flux-radius plane (see Sect. 3.1), the errors on the fitted parameters are relatively larger because the statistics in the flux-mass diagram are larger, as already discussed in Sect. 2.3. Within the S₃ sample, 110 planets are found to populate the desert obtained by Mazeh et al. (2016) in the P – M_p diagram. Figure 11 depicts the distribution of these planets in the F – M_p diagram, together with the lower and upper boundaries obtained in this work. As discussed in Sect. 2.3, because the statistics in the flux-mass plane are lower, the confidence bands of the edges overlap when we take the uncertainties on the fit in the flux-mass plane into account, as shown in Fig. 11. However, 63 planets out of 110 still receive an incident flux F ≤ 550 F_⊕ and therefore do not populate the irradiation desert in the F – M_p plane. This might be the result of using a different F_75th between the F – R_p and F – M_p planes because this effect does not clearly arise when we considered the orbital period in place of the incident flux.

In Sect. 2.4.2 we constructed an augmented version of sample S₂ by using forecaster to obtain an M_p estimate for planets without mass measurements and keeping objects regardless of the magnitude of the uncertainty ∆M_p /M_p. In principle, we could have used this sample to compute the analytical expression of the lower and upper boundary in the F – M_p plane. However, even though this approach would yield a larger data set to be fitted with DIAMONDS, we definitely lose much in accuracy. Typically, the uncertainties in M_p estimated by forecaster range from approximately 30 to 50% for well-characterised exoplanets, but this can be significantly higher for planets with less precise data. In particular, for our data set, we found the uncertainties to be normally distributed around ≈60% with a standard deviation of about 5%. Consequently, even though employing an augmented sample might help us to mitigate the issues associated with limited statistics in the F – M_p plane, the results would ultimately rely on a highly uncertain sample, which would lead to unreliable boundary points.

Fig. 10 Comparison of the S–1 sample distribution in the P – Rp and F – Rp planes to highlight the irradiation desert. Left: distribution of the S–1 sample in the P – Rp diagram along the edges obtained by Mazeh et al. (2016). The yellow triangles represent the 221 hot Neptunes within the S1 sample according to the definition of Mazeh et al. (2016) in the P – Rp plane. The black dots represent all the planets that lie outside the desert. Right: distribution of the S1 sample in the F – Rp plane. The orange and blue lines represent the lower and upper boundary obtained in Sect. 3.1 alongside their confidence bands (1σ). The yellow triangles show the 194 out of 221 objects that do not lie within the irradiation desert defined in this work. The yellow triangles used to designate the hot Neptunes in the P – Rp diagram still refer to the same planets in the right plot. The red stars indicate planets that fall within our defined boundaries but not within those outlined by Mazeh et al. (2016).

3.3 Planetary mean density

In Sects. 3.1 and 3.2 we obtained the analytic expressions of the boundaries of the sub-Jovian and Neptune desert in the F – R_p and F – M_p planes. For each edge, by combining R_p (F) and M_p (F), we can obtain a mass-radius relation in order to understand what types of planet populate them. In particular, the mass-radius relation for the upper boundary is given by $$\frac{R_p}{R_{\oplus }}=(3.1\pm 2.5){\left(\frac{M_p}{M_{\oplus }}\right)}^{0.23\pm 0.15},$$(14)

while the mass-radius relation describing the lower boundary is $$\frac{R_p}{R_{\oplus }}=(0.7\pm 0.7){\left(\frac{M_p}{M_{\oplus }}\right)}^{0.39\pm 0.07}.$$(15)

We show the mass-radius relations for the two boundaries in Fig. 12. Planets populating the upper edge of the sub-Jovian and Neptune desert in the F – R_p and F – M_p planes are compatible with very low planetary density values (≈0.1 g/cm³), thus suggesting a gaseous nature. The resulting exponent (0.23 ± 0.15) agrees well with the findings of Mazeh et al. (2016) within the numerical uncertainties. On the other hand, planets lying on the lower edge systematically have a higher density. The spread of the curve, which is mostly due to the uncertainty on the mass estimates, leads to a certain level of ambiguity: These planets could either be covered by a gaseous envelope that drastically reduces the mean planetary density, or they might be silicate worlds. The former case might be compatible with a low-mass planet covered by a H/He envelope that undergoes photoevaporation due to the stellar irradiation of the host. The latter might represent the final stage of the previous scenario where the planet has lost part of its envelope.

Fig. 11 Comparison of the S3 sample distribution in the P – Mp and F – Mp planes to highlight the irradiation desert. Left: distribution of the S3 sample in the P– Mp diagram along the edges obtained by Mazeh et al. (2016). The yellow dots represent the 110 hot Neptunes within S3 according to the definition of Mazeh et al. (2016) in the P – Mp plane. The black dots represent planets outside the desert. Right: distribution of our sample in the F – Mp plane. The orange and blue lines represent the lower and upper boundary obtained in Sect. 3.2 alongside their confidence bands (1σ). The red dots show that 63 out of 110 objects experience an incident flux lower than 550 F⊕. The yellow triangles used to designate the hot Neptunes in the P – Mp diagram still refer to the same planets in the right plot.

Fig. 12 Mass-radius relations for the upper (orange line) and lower boundary (blue line) given by Eqs. (14) and (15), respectively. We present different curves at fixed density (left). While the upper boundary is consistent with a very low mean planetary density (≈0.1 g/cm3), the lower boundary comprises different values spanning from ≈1 g/cm3 to ≈2g/cm3. We also present different mass-radius curves for different planetary compositions (right).

4 Conclusions

We have expanded upon the seminal work of Mazeh et al. (2016) by studying the sub-Jovian and Neptune desert in more detail. This is the region of the P – R_p and P – M_p space parameters from which Neptunian exoplanets with short orbital periods are notably absent. The orbital period is not a complete physical variable for understanding the nature of this phenomenon because it does not account for stellar radiation. Our analysis therefore incorporated the incident flux from the host star in order to provide a more comprehensive picture of the factors influencing the formation and evolution of these intriguing exoplanets.

Our analysis revealed that the incident flux plays an important role in shaping the boundaries of the sub-Jovian and Neptune desert. We analysed the dependences of the upper and lower edges on the incident flux in the flux-radius and flux-mass space parameters. We employed DIAMONDS, a novel Bayesian nested-sampling tool, to infer the best models that fit the two boundaries in both planes. Even when we accounted for the incident flux, we still found a paucity of known Neptune- and Saturn-like exoplanets experiencing high-intensity levels of stellar radiation. We called the region in which they lie the irradiation desert. Moreover, ≈87% of planets in the desert defined by Mazeh et al. (2016) in the period-radius plane do not lie within the irradiation desert we defined in the flux-radius diagram. A similar result was obtained in the flux-mass plane. Thus, the sub- Jovian and Neptune desert is likely much more desolate than we would expect. We stress that our method yielded boundaries that delimit a region that admits contamination by a few moderately irradiated Neptunes and sub-Jovians.

We also retrieved the mass-radius relation for the planets lying close to the two boundaries. While the upper boundary is populated by exoplanets with a very low mean planetary density, which is compatible with gas giants that are stable against the photoevaporation scenario, the lower boundary comprises a broader range of planetary types: from planets with gaseous envelops to objects with a higher mean planetary density, which might be compatible with worlds dominated by silicates.

Future observations with advanced telescopes such as the James Webb Space Telescope (JWST, Gardner et al. 2023) and exoplanetary missions, as well as an ongoing analysis of the thousands of planet candidates that have yet to be confirmed, will help us to refine the calculations of the two boundaries and will thus allow us to potentially reveal new insights into the nature and origins of Neptune-sized exoplanets in close orbits. In particular, the PLAnetary Transits and Oscillations of stars (PLATO, Rauer et al. 2014, 2024) mission is expected to discover several hundreds to thousands of Neptune-sized exoplanets orbiting their star in fewer than 5 days, depending on the assumed planet occurrence rates (Matuszewski et al. 2023). PLATO will not only enhance our understanding of planetary architectures by expanding the exoplanetary populations, but its high-photometric precision, combined with a host star characterisation via asteroseismology, will allow us to refine the planetary and stellar parameters. This will reduce the uncertainty on the data used in this work.

As the volume of these discoveries increases, there will be a critical need for efficient vetting processes, both through human expertise and automated machine-learning techniques (e.g. McCauliff et al. 2015; Coughlin et al. 2016; Mislis et al. 2016; Shallue & Vanderburg 2018; Yu et al. 2019; Tey et al. 2023; Fiscale et al. 2023), to ensure the accuracy and reliability of confirmed planets. However, the complexity of machinelearning models, which are often perceived as black-box systems, can limit transparency and scientific understanding. Therefore, integrating explainable machine-learning methods (Maratea & Ferone 2019) into future vetting pipelines will be essential. These techniques will allow scientists to maintain control over and insight into the decision-making process, enabling greater confidence in the discoveries and fostering collaboration between computational tools and human expertise in the validation of new exoplanet candidates.

Acknowledgements

This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. We thank the anonymous referee whose careful reading helped to improve the quality of our manuscript. E.C. acknowledges support from PLATO ASI INAF agreement no. 2022-28-HH.0 “PLATO Fase D”. G.C. thanks the University of Napoli Federico II (project: FRA-CosmoHab, CUP E65F22000050001). Software: DIAMONDS (Corsaro & De Ridder 2014), forecaster (Chen & Kipping 2016).

References

All Tables

All Figures

Fig. 1 Distribution of the confirmed exoplanets of the S1 sample in the log P – log Rp diagram for different ranges of the insulation flux. The upper left panel shows the 586 planets that receive ≤50 F⊕, the upper right panel shows the 364 planets that receive between 50 F⊕ and 200 F⊕, the lower left panel depicts the 238 planets that receive between 200 F⊕ and 550 F⊕, and the lower right panel illustrates the 339 planets with F > 550 F⊕.

In the text

Fig. 2 Analysis of the planetary distribution in the F – Rp plane and the flux cumulative distribution function. Left panel: distribution of the 1527 planets of the S1 set in the F – Rp plane. The yellow triangles represent planets that would be referred to as hot Neptunes according to the limits obtained by Mazeh et al. (2016). The grey triangle is a qualitative representation of the desert. Right panel: cumulative distribution function of incident flux values F of planets with 3 R⊕ ≤ Rp ≤ 11 R⊕. The vertical dashed line in both panels represents the threshold value of 200 F⊕ below which lie 75% of the planets with 3 R⊕ ≤ Rp ≤ 11 R⊕.

In the text

	Fig. 3 Application of a two-component bivariate GMM over S1 sample in the log F – log Rp space. We found two separate populations: group AR (purple dots) with a mean value of (1.70,0.34), and Group BR (yellow dots) with a mean value of (2.72, 1.12).
In the text

	Fig. 4 Calculated boundaries in (F, Rp) space: (Flow,i, Rlow,i) pairs (purple points) sample the lower boundary, and (Fup,i, Rup,j) pairs (red points) mark the upper boundary.
In the text

Fig. 5 Analysis of the planetary distribution in the F – Mp plane and the flux cumulative distribution function. Left panel: distribution of the 609 planets of the S2 set in the F – Mp plane. The yellow triangles represent planets that would be referred to as hot Neptunes according to the limits obtained by Mazeh et al. (2016). The grey triangle is a qualitative representation of the sub-Jovian and Neptune desert. Right panel: cumulative distribution function of incident flux values F of planets with 10 M⊕ ≤ Mp ≤ 200 M⊕ . The vertical dashed line represents the threshold value of 550 F⊕ in both panels below which 75% of the planets with 10M⊕ ≤ Mp ≤ 200 M⊕ lie.

In the text

	Fig. 6 Application of a two-component bivariate GMM on the sample S2 in the log F – log Mp space. We identified two separate populations: group AM (purple dots) with a mean value of (1.8,1.4), and group BM (yellow dots) with a mean value of (2.9, 2.5).
In the text

	Fig. 7 Calculated boundaries in the (F, Mp) space: (Flow,i, Mlow,i) pairs (purple points) sample the lower boundary, and (Fup,i, Mup, j) pairs (red points) mark the upper boundary of the desert.
In the text

	Fig. 8 Histogram of the variable F75th for planets with 3 R⊕ ≤ Rp ≤ 11 R⊕, generated from 10 000 Monte Carlo simulations. Each simulation reduced the S1 sample to the size of the S2 and calculated the 75th percentile of exoplanets with 3 R⊕ ≤ Rp ≤ 11 R⊕ . The histogram is normalised to show the probability density.
In the text

Fig. 9 Distribution of 512 objects of sample S3 within the log Rp – log Mp plane, coloured according to their classification with respect to the radius and mass. 271 blue dots (ARAM) represent the planets classified as not giants for the radius and mass, 211 red dots (BRBM) represent the planets labelled as giants for the radius and mass, 29 green points (BRAM) are giant planets according to their radius and not giants according to their mass, and the one yellow point (ARBM) corresponds to a not-giant in terms of the radius and to a giant with respect to the planetary mass.

In the text

Fig. 10 Comparison of the S–1 sample distribution in the P – Rp and F – Rp planes to highlight the irradiation desert. Left: distribution of the S–1 sample in the P – Rp diagram along the edges obtained by Mazeh et al. (2016). The yellow triangles represent the 221 hot Neptunes within the S1 sample according to the definition of Mazeh et al. (2016) in the P – Rp plane. The black dots represent all the planets that lie outside the desert. Right: distribution of the S1 sample in the F – Rp plane. The orange and blue lines represent the lower and upper boundary obtained in Sect. 3.1 alongside their confidence bands (1σ). The yellow triangles show the 194 out of 221 objects that do not lie within the irradiation desert defined in this work. The yellow triangles used to designate the hot Neptunes in the P – Rp diagram still refer to the same planets in the right plot. The red stars indicate planets that fall within our defined boundaries but not within those outlined by Mazeh et al. (2016).

In the text

Fig. 11 Comparison of the S3 sample distribution in the P – Mp and F – Mp planes to highlight the irradiation desert. Left: distribution of the S3 sample in the P– Mp diagram along the edges obtained by Mazeh et al. (2016). The yellow dots represent the 110 hot Neptunes within S3 according to the definition of Mazeh et al. (2016) in the P – Mp plane. The black dots represent planets outside the desert. Right: distribution of our sample in the F – Mp plane. The orange and blue lines represent the lower and upper boundary obtained in Sect. 3.2 alongside their confidence bands (1σ). The red dots show that 63 out of 110 objects experience an incident flux lower than 550 F⊕. The yellow triangles used to designate the hot Neptunes in the P – Mp diagram still refer to the same planets in the right plot.

In the text

Fig. 12 Mass-radius relations for the upper (orange line) and lower boundary (blue line) given by Eqs. (14) and (15), respectively. We present different curves at fixed density (left). While the upper boundary is consistent with a very low mean planetary density (≈0.1 g/cm3), the lower boundary comprises different values spanning from ≈1 g/cm3 to ≈2g/cm3. We also present different mass-radius curves for different planetary compositions (right).

In the text