© The Authors 2023

1 Introduction

The observed absence of short-period planets (≲3 days) below the hot Jupiter clump and above the hot super-Earths is known as the “sub-Jovian desert” (Szabó & Kiss 2011; Benítez-Llambay et al. 2011; Sanchis-Ojeda et al. 2014; Colón et al. 2015; Matsakos & Königl 2016; Eigmüller et al. 2017; Owen & Lai 2018; Szabó & Kálmán 2019) and “Neptunian desert” (Mazeh et al. 2016; Demangeon et al. 2018; Ionov et al. 2018; Mori et al. 2022). In recent years, several planets have been observed in these desert regions (Demangeon et al. 2018; Dragomir et al. 2020; Jenkins et al. 2020; Armstrong et al. 2020; Murgas et al. 2021; Mori et al. 2022), with some deserts possibly undergoing a conversion into a “savanna.” In this paper, we describe the desert–savanna boшidaries in the parameter space.

There is a known planet population both above and below its boundaries, therefore, the sub-Jovian or Neptunian desert or savanna itself is a puzzling structure of planetary distribution. The location of the desert is a meeting point for high-energy physics and atmospheric non-local thermodynamic equilibrium (NLTE) processes (representing the stellar activity and irradiation), planetary thermodynamics, and magneto-hydrodynamical interactions between the planet and the star, end-points of planetary migration, accretion processes in general, and planet formation and evolution in extreme environments as well. Therefore, a test for planet formation theories is to check how they can explain the presence of the desert and its boundaries.

In Fig.1, we summarize the possible paths of how a young planet, once formed, can leave the desert. The indicative direction of “leaving the desert” is denoted by numbers in the figure. The processes suggested for the different scenarios thus far include: (1) the hyperinflation of low-mass gas giants at the boundary (Mordasini et al. 2015); (2) planet-growth-based processes, namely, it is has been well established that hot Jupiters have higher-than-expected (i.e., inflated) radii (Thorngren & Fortney 2018), which can be attributed to irradiation from the host (e.g., Sarkis et al. 2021). As the desert is observed in both the R_P–P and M_P–P planes (and at sub-Jovian planet sizes), these processes would necessarily have to include accretion as well as inflation. The theoretical calculations done so far (Emsenhuber et al. 2021a,b) were insufficient to explore this scenario in detail; (3) high-eccentricity migration (Giacalone et al. 2017; Owen & Lai 2018) suggests that the orbits of close-in planets with initially high eccentricities circularize with higher orbital periods; (4) decreasing planet size can be the result of photoevaporation (Lopez & Fortney 2014; Lundkvist et al. 2016; Owen & Lai 2018), Roche-lobe-overflow (Kurokawa & Nakamoto 2014), “boil-off” (Owen & Wu 2016), and impact-based mass-loss (Schlichting etal. 2015; Schlichting & Mukhopadhyay 2018). Finally, (5) the tidal migration of hot Neptunes (Lin & Papaloizou 1986; Rozner et al. 2022), tidally trapped outward migration (Masset et al. 2006), and disk migration (Mazeh et al. 2016; Ataiee & Kley 2021) are other alternative scenarios.

Besides scenarios involving migration and evaporation, Bailey & Batygin (2018) found that an accretion parameter correlated with planet mass and disk inner edge can explain the upper right side of the desert and explains its suggested slope of –2/7. Ataiee & Kley (2021) suggested a further possibility, where orbital resonances change the orbital periods of the inner planets. The curious case of TOI-2196 b (Persson et al. 2022) would fit in well with migration-based scenarios, being a small and volatile-rich planet present in the savanna.

Mazeh et al. (2016) derived relationships for the lower and upper boundaries of the desert in the R_P–P plane finding that at the upper edge, $R_{\text{P}}\propto P^{-\frac{1}{3}}$, and at the lower edge, $R_{\text{P}}\propto P^{\frac{2}{3}}$, thus suggesting different origins in the two size-regimes. Matsakos & Königl (2016) explained the two boundaries as the tidal disruption barrier for gas giants following their high-eccentricity migrations and because of different mass-radius laws of low-mass and high-mass exoplanets, the two edges are formed differently. Owen & Lai (2018) suggested that the photoevaporation and the high-eccentricity migration to form the upper and lower boundaries, respectively. Owen (2019) summarizes the processes that could have formed the desert boundaries from an evaporation approach, with a very detailed introduction to the astrophysical processes. A somewhat similar structure in the period-radius distribution is related to the observed increasing sub-Saturn (≈4–8 R_⊕) occurrence rate up to ~300 days orbital period. Hallatt & Lee (2022) argued for a radiation-induced process (resulting in atmoshperic mass-loss), which is efficient at transforming sub-Saturns into sub-Neptunes (≲4 R_⊕), with short orbital periods. The short-period end of this process is similar to the “process 3” in our Fig. 1 and could explain the lower boundary of the desert as well.

In our earlier paper Szabó & Kálmán (2019), we found that the boundaries of the desert depend more on fundamental stellar parameters and less on planetary mass, density, and size. In that work, the sample of exoplanets were smaller, which limited the further inference. In this paper, we re-apply the same analyses to the larger exoplanet data set at hand and we include a few new tests as well. We show that the boundaries of the desert in the period-mass (P–M_P) and the period-radius (P–R_P) planes are marked by a different planet population. The desert in the two main projections can be considered as a sign of multiple scenarios of planetary evolution acting differently in the various parameter spaces. Therefore, we revisited the P–M_P and P–R_P boundaries separately, which allowed for the findings of Szabó & Kálmán (2019) to be quantified. The dependence of the P–M_P and P–R_P boundaries on the stellar parameters are to be considered as another piece of evidence to support this assertion. Finally, we give expressions for both the planetary radius and mass in terms of stellar parameters and via fuzzy clustering. We conclude that in all projections, planets form two main groups that behave differently, in a possible connection with the unresolved processes behind the formation of the sub-Jovian or Neptunian desert.

Fig. 1 Sketch of possible paths for already formed planets leaving the sub-Jovian desert (red arrows) and planets not allowed to reach the desert (blue arrows) overplotted to the distribution of known exoplanets. Theories and references behind the scenarios are summarized in the text.

2 Methods and data selection

2.1 Data on confirmed planets from the NASA exoplanet archive

The input to this analysis consists of two data sets. We followed Szabó & Kálmán (2019) in building up the sample from the NASA Exoplanet Archive (5009 planets in total¹), following the same filtering as in Szabó & Kálmán (2019). This filtering has left out planets with less constrained parameters and kept the “golden sample” only. In building up the test sample, we defined the following selection criteria: (i) planetary mass M_P < 13 M_J and density ρ_P < 25ρ_J and (ii) relative uncertainty for the planetary radius and mass of <20% and <60%, respectively, leaving us with a sample of 650 exoplanets. This is an increment with regard to the 406 planets in our previous analysis (Szabó & Kálmán 2019) and our current sample also represents a better quality because it also contains the filtering for the mass determination errors. We accepted the stellar parameters provided in the original NASA Exoplanet Archive records.

In Fig. 1, we show the distribution of our filtered sample (of 650 exoplanets) in the R_P–P plane. Figures 2–8 show the same sample in both the radius-period and mass-period parameter spaces. This design enables the visualization of the exact values of the continous variable, but has the limitation of being merely qualitative. To be able to quantitatively see the structures in the sample, we applied two-sampled Kolmogorov-Smirnov (KS) tests (see e.g., Feigelson & Babu 2012) to the projected parameters in different regions (univariate bands) of the data distribution.

Still following Szabó & Kálmán (2019), we defined two statistical regions in the R_P–P plane as: the area between 0.28 and 0.63 R_J is labeled A, and the one between 0.06 and 0.16 R_J is labeled B. In the M_P–P plane, we defined the С region between 0.032 and 0.305 M_J and the D region between 0.001 and 0.010 M_J.

To observe whether a stellar or planetary parameter is selective at the position of the border (in other words, the position of the border depends on the examined parameter), we compared the period distribution of two subsamples in A to those in B; and two subsamples in С to those in D. The subsamples were separated at the median of the examined parameter. For example, we compared the period distribution of high-temperature and low-temperature stars in A and B, etc. We compared the resulting period distributions in a two-sampled KS test in all four regions and computed the p values. If the resulting p is small, we can conclude that the two compared samples are drawn from different distributions; or, in other words, the third parameter can affect the distribution of exoplanets. Because regions A and С are in the desert and В and D are outside, this comparison helps identify those parameters that affect the formation and the boundaries of the desert.

Following the technique of colored scatter plots, we sometimes identified several groups of exoplanets following a specific pattern. These features can be proven with clustering methods. Here, we followed the fuzzy clustering algorithm of Bezdek (1981), as implemented in the FKM function in the R software package².

Fig. 2 Sub-Jovian/Neptune desert of exoplanets in the RP–P and MP–P parameter spaces (top and bottom panel), colored by the planetary density. In the top panel, the upper boundary is marked by “normal” Jupiters with higher density (>1ρJ), while in the bottom panel, the population of inflated planets (red points) appear below the normal Jupiters and, thus, the inflated Jupiters mark the actual boundary.

Fig. 3 Scatterplots of the distribution of exoplanets in the RP–P and MP–P planes as in Fig. 2 but with synthetic densities subtracted from the observed ones (see text for details).

Fig. 4 Scatterplot showing the distribution of exoplanets in the RP–P (top left) and MP–P (top right) parameter spaces, colored by the temperature of the host star (see the color bar for the values). The middle and bottom panels show the cumulative distribution of the two samples divided at the median of Teff in the A, B, C, and D regions. The p-values of the KS-test are also shown.

Fig. 5 Scatterplots showing the distribution of exoplanets in the RP–P and MP–P parameter spaces with Rs used as the third parameter (top row) and the respective cumulative distributions (middle and bottom rows; same as Fig. 4).

Fig. 6 Scatterplots showing the distribution of exoplanets in the RP–P and MP–P parameter spaces with MS used as the third parameter (top row) and the respective cumulative distributions (middle and bottom rows; same as Fig. 4).

Fig. 7 Same as Fig. 4 but with log ɡ used as the third parameter instead of Teff.

2.2 APOGEE catalog of planet-host stars

The second database we used for the analyses is based on Apache Point Observatory Galactic Evolution Experiment (APOGEE, Majewski et al. 2017), which provides homogeneously derived atmospheric stellar parameters and eliminates uncertainties originating from the usage of different model atmospheres, spectral synthesis codes, line lists, and so on. We built an APOGEE exoplanet host stars catalog to test the findings on this homogeneous data set and we plot the data in Figs. 10 and 11.

As part of the fourth iteration of Sloan Digital Sky Survey (SDSS-IV, Blanton et al. 2017), in its final data release (DR17, Abdurro’uf et al. 2022), the APOGEE-2 survey derived and reported main atmospheric parameters and individual element abundances for 733 901 stars across the entire sky (Zasowski et al. 2017) and focused on observing as many planet-host stars as possible, so it is an excellent choice to cross-match our sample with that of DR17. Observations were made with the identical SDSS spectrographs (Wilson et al. 2019) mounted on the 2.5-m Sloan Foundation Telescope (Gunn et al. 2006) and the 2.5-m Irénée du Pont Telescope of Las Campanas Observatory, Chile. The resolving power is R ~ 22 500 and the spectral coverage of the near-infrared H-band – from 15 140 Å to 16 940 Å. This makes it possible to select planet-host stars near the Milky Way disk and dust obscured regions of the Milky Way. Besides observing stellar targets, APOGEE is sensitive to substellar mass companions (exoplanet candidates) too, as multi-epoch visits along with high-precision RV measurements are available (Majewski et al. 2017).

The APOGEE DR17 collaboration published both raw and calibrated values of effective temperature, surface gravity, and applied offset corrections to [M/H] based on solar neighborhood stars. Here, we take the calibrated effective temperatures and overall metallicities into account as T_eff is in a better agreement with the infrared flux method (IRFM) scale than the raw spectroscopic values (Abdurro’uf et al. 2022). For this analysis, we carefully select the most accurate and precise APOGEE data. We globally filter out APOGEE stars labeled by the quality control flag STAR_BAD³ (Abdurro’uf et al. 2022, but retained STAR_WARN), then we cross-matched the exoplanet catalog with the APOGEE DR17 data set in TOPCAT version 4.8−2 (Tool for Operations on Catalogues And Tables, Taylor 2005) based on the equatorial coordinates of the host stars (∆ < 1 arcsec). The next step involved the implementation of several quality criteria: relative error lower than 20% for R_P, VSCATTER < 1 km s⁻¹ for RV variations, VERR < 1 km s⁻¹ for RV error to filter out potential variable stars, and M_H_ERR < 0.1 dex for the metallicity uncertainty.

We note that the selected elements fall into the APOGEE data reduction categories of “most reliable” and “reliable” with regards to the overall quailty of their derivation (see the APOGEE DR17 website for more details).

Fig. 8 Same as Fig. 4 but with [M/H] used as the third parameter instead of Teff.

Fig. 9 Correlation between planet mass (top row), planetary radius (bottom row), and the main stellar parameters: Teff (left column), MS (middle column), and [M/H] (right column). The samples are divided into two groups via a fuzzy clustering algorithm applied to the distributions in the right panels. Blue points: Cluster of “larger planets.” Red points: Cluster of ‘smaller planets.” The fitted linear trends of Eqs. (2)–(11) are plotted with solid blue and red lines. Pearson’s r values, corresponding to the giant and smaller exoplanets are also displayed with blue and red, respectively.

Fig. 10 Correlations between RP–Teff (left panel) and RP–[M/H] (same as the lower left and lower right panels of Fig. 9) showing the fuzzy clustering of APOGEE data and with the linear regressions of Eqs. (14)–(16) overplotted.

Fig. 11 Distribution of exoplanets from the APOGEE planet host sample with [M/H] as the coloring parameter (left panel) and the residuals after a bilinear fit (right panel).

3 Results

3.1 Different planets on the borders in projections

In Fig. 2, we show the filtered NASA Exoplanet Archive planets in the R_P–P and M_P–P planes with the planet’s mean density as a coloring variable. The comparison of the two scatter plots reveals an interesting feature. Planets with ≳ρ_J density, plotted as blue dots, follow a similar distribution in both planes, with the most important difference of the compressed size range of hot Jupiters (upper panel) that is related to well-understood planetary physics. However, the group of low-density planets (plotted in red) behave very differently between the two plots. These planets have ≲M_J mass and they can be found at the middle of the mass distribution (lower panel), but they are among the most inflated planets, and can be observed on the top of hot Jupiters in the radius distribution (upper panel).

This distinct clump of planets suffers from a very significant dislocation compared to the non-inflated or moderately inflated planets, which requires explanation. Also, the low-density planets themselves form the upper boundary of the sub-Jovian (or Neptunian) desert in the M_P–P projection, while they tend to be farther from the boundary than the hot Jupiters in the R_P–P plane. Indeed, in the R_P–P plane, the hot Jupiters themselves form the upper boundary of the desert.

The anonymous referee has pointed out that similar structures could be simple byproducts of the well-known increasing bulk density above a 1.0 M_J planetary mass, due to the increasing self-compression of the gas that builds up the planetary body. To test the robustness of this finding, we used an empirical model to describe the mass-density relations in general and we removed this contribution from the scatter plot.

We calculated synthetic densities for each planet, following Traub (2011); Mordasini et al. (2012) and based on mean radii (Eq. (24) of Mordasini et al. 2012). Since the applied model equation is designed to reproduce theoretical radii for planets, with semi-major axis a < 0.1 AU and 0.1 AU <a < 1 AU, we further restricted our sample to include planets with a < 1 AU. Subtracting the logarithm of the resulting densities yielded the distribution shown in Fig. 3. We find that the upper edge of the desert in the M_P–P plane (bottom panel of Fig. 3) is still generally dominated by low residual-density planets, which are then then observed to be further away from the sub-Jovian desert in the R_P–P plane (top panel of Fig. 3).

This suggests that the general density-mass relation does not explain the observed phenomenon. Indeed, the scatter plots colored with the residual log densities (Fig. 3) show a very similar distribution to what we see in Fig. 2 and we can have the subjective visual impression that Fig. 3 is even clearer. In any case, our conclusion is that the found double population of close-in giants is not a consequence of the density-mass relations, but it is rather superimposed onto the general pattern of density distributions. Taking the position of the planets with ρ/ρ_J < 0 into account, it seems to be quite plausible that their presence is primordially influenced by the desert boundaries. We also point out that the exoplanets with longer orbital periods (> 10 days) which have lower densities (red points on Fig. 2) remain roughly in the same position relative to the desert in both planes of Fig. 2.

The analysis of the bulk densities (Fig. 2) thus suggests the upper boundary of the desert has a multiform nature: it is marked by ordinary hot Jupiters in one projection and inflated hot Jupiters in the other projection. Since the planetary physics behind these two distinct classes of hot Jupiters are different, this also brings up the possibility that the formation of the desert is a result of multiple processes (at least at the upper boundary): one process acting on ordinary jot Jupiters and one another acting on inflated hot Jupiters.

The possible biases behind the observed distribution can be excluded by simple quantitative considerations. In particular, 137 exoplanets falling into the ρ_P < 0.35 ρ_J category were discovered primarily with 18 different photometric facilities including both ground-based and space telescopes, so their distribution in the parameter-space cannot be devoted to some instrument-specific selection distortion of some specific exoplanet discovery surveys. Also, from the entire sample of 650 exoplanets, only 23 of them were not discovered with the transit method (mostly radial velocity planets that had later been observed in transit as well) and none of them belong to the inflated group of planets (according to the definition given at the beginning of this paragraph). Therefore, the presence of radial velocity discoveries should not be considered as sources of biased data. Out of the 137 exoplanets with ρ_P < 0.35 ρ_J, an overwhelming majority (116) are members of planetary systems with a single known planetary mass object, with 19 systems containing 2 planets and 2 containing 2 and 3 planets. The brightness distribution of host stars in the two samples is also very similar, with V = 10.56 ± 1.94 mag for the larger, 650-member sample and 10.89 ± 1.59 mag for the low-density sample. We therefore see no evidence for observational biases. We can consider two possible explanations for the multiple planet populations at the upper boundary: (i) either the inflated planets are in the early stages of planetary evolution and they are still in the process of initial contraction or losing their atmospheres; or (ii) they are more mature planets, but have a different inner structure than the higher density hot Jupiters. In either case, the processes that form the desert must be selective for these formation scenarios, leading to the two distinct planet populations at the two different projections of the upper boundary.

To test these scenarios, a wider sample with more precisely determined exoplanet parameters is necessary. Our understanding of these questions will be greatly improved by the two upcoming ESA missions: Ariel (Tinetti et al. 2021) and PLATO (Rauer et al. 2014).

3.2 Dependence of the boundary on stellar parameters

In this subsection, we analyze the RP-P and MP-P scatter plots with regard to the distribution of stellar parameters (effective temperature, stellar radius, stellar mass, log g, and metallicity), complementing and extending the results described in Szabó & Kálmán (2019).

3.2.1 Effective temperature of the host

Following the recipe described in Sect. 2, we split the sample of exoplanets at the median of the effective temperature of the host stars in different regions separately. The split values were 5370 K, 4677 K, 5370 K, and 3388 K for the A, B, C, and D regions, respectively (see Fig. 4). We found that in both parameter spaces, planets in A and C regions tend to have shorter orbital periods around cooler hosts. The p-values of the KS tests were: 10⁻⁵ and 10⁻⁴ in A and C, respectively. This discrepancy in the cumulative distribution has not been observed in the B or D regions (p-values: 0.45 and 0.07), proving that the observed dependence is specifically located at the right boundary of the desert. Therefore, we confirm the findings of Szabó & Kálmán (2019) that the effective temperature of the host star plays a key role in the formation of the desert and it does so in both parameter spaces.

3.2.2 Stellar radius

The boundary of the desert also shows dependence on R_S in both the R_P–P and M_P–P planes, as visible in the scatter plots of Fig. 5. The stars below the median radii (0.94 R_⊙ and 0.90 R_⊙ in the A and C regions, respectively) are more likely to host planets with shorter orbital periods. We note that the KS test yields p-values of 2 × 10⁻⁴ and 10⁻⁵ in A and C, respectively. There is no significant difference in the cumulative distribution planets in the B and D regions. In B, the p-value is 0.43; while in D, the p = 0.013 still corresponds to 2.48σ, proving the statistical insignificance of the (otherwise apparent) differences seen in the bottom right panel of Fig. 5. As we find that stars with lower radii tend to host planets with shorter orbital periods and that this phenomenon is not present in the control regions, we conclude that the stellar radius also plays a key role in the formation of the desert.

3.2.3 Stellar mass

The scatter plot of exoplanets with the mass of the host star used as the third parameter is shown in Fig. 6. The split values of M_S are 0.92 M_⊙, 0.76 M_⊙, 0.92 M_⊙, and 0.34 M_⊙ in A, B, C, and D, respectively. Dividing the samples in the four regions at these values leads to an ambiguous detection of M_S dependence of the desert boundary (i.e., the distribution of planets at the size ranges of the desert). This is expressed in terms of the p-values of 6 × 10⁻⁵ opposed to 0.45 in A and B regions, and then 0.055 opposed to 0.013 in C and D regions. The dependence is present in the R_P–P plane and is insignificant in the M_P–P plane, consistently with the results of Szabó & Kálmán (2019). Therefore, we find that stellar mass has some influence most importantly on the radius distribution, but does not play an essential role in the formation of the desert on its own.

3.2.4 Stellar log ɡ

Exploring the effects of log ɡ on the distribution of exoplanets in the two parameter spaces, we divided the samples at log ɡ = 4.46 in the A region, log ɡ = 4.60 in B, log ɡ = 4.48 in C, and log ɡ = 4.87 in D. We find that there is no statistical difference in the cumulative distributions in either the A or B regions of the R_P–P plane. This is expressed by the p-values of the KS-test: 0.040 in A (corresponding to ~2.1σ) and 0.22 in B. In the size ranges of the desert of the M_P–P parameter space, there is a hint of a discrepancy between the cumulative distributions of the sample divided at the median at p = 0.0055 (corresponding to 2.78σ). There is no such difference in the sample of the control region (p = 0.34).

We therefore find no firm dependence of the desert on the surface gravity of the host star, contradicting to the earlier results of Szabó & Kálmán (2019; based on fewer planets and unconstrained for a reliable mass determination). The structure of Fig. 7 the most inflated hot Jupiters are hosted by stars with the lowest log g values (<3.2) from the sample.

3.2.5 Stellar metallicity

In the filtered sample of 650 used in the exploration of the previous parameters, there are 10 stars with unknown metallicities and we omitted those from the current analysis. The samples in A, B, C, and D are split at [M/H] = 0.16, 0.01, 0.11, and −0.01, respectively. There is no significant difference in the cumulative distribution of planets in the two examined areas of the R_P–P plane (p-values are 0.088 and 0.27 in A and B). In the case of the M_P–P plane, there is also only a hint that planets from the C region are more likely to have shorter orbital periods around younger, more metal-rich stars, expressed via p = 0.0075 (corresponding to 2.67σ) in C compared to 0.30 in D.

As in the case of log ɡ, we find no firm dependence of the sub-Jovian desert of exoplanets on the metallicity of the host stars. This is somewhat contradictory to the results of Dong et al. (2018); Petigura et al. (2018) and Szabó & Kálmán (2019).

3.3 Multimodal dependence of the planetary mass and radius on stellar parameters

The experienced bimodality of the planet groups at the border of the desert (Fig. 4) and the experienced differences behind the R_P–P and M_P–P relations reflect the presence of two different type of planets at the border: hot Jupiters at the upper boundary and super-Earths at the lower boundary. As the boundary of the desert in the mass and radius projections depend on stellar parameters, we expect to see various groups of planets with different parameter dependencies between the radius, mass, and stellar parameters.

To explore distinct correlation laws between the planet and stellar parameters, we divided our sample into two subgroups in both the mass-metallicity and the radius-metallicity planes via a k-means fuzzy clustering algorithm (Ferraro et al. 2019). In the radius-metallicity plane, the resulting cluster consisting hot Jupiters had a median radius of 1.15 R_J (interquartile range of 1.01–1.33) and a cluster consisting Neptune- and Earthsized planets had a median radius of 0.21 R_J (interquartile range of 0.14–0.27). When the clustering was done in the massmetallicity distribution, the cluster of large planets had a median mass of 0.93 M_J (interquartile range of 0.57–1.80) and the subgroup of smaller planets has a mean mass of 0.026 M_J (interquartile range of 0.015–0.048).

We fit linear regression models describing the dependence of planet mass and radius on T_eff, M_S, and [M/H] (Fig. 9). Between the planet size and the effective temperature, we get: $$\begin{array}{cc}\log \left(\frac{M_{\text{p}}}{M_{\text{J}}}\right)=1.91\left(34\right)\cdot \log \left(\frac{T_{\text{eff}}}{1\text{\hspace{0.17em}\hspace{0.17em}K}}\right)-7.15\left(1.27\right),& r=0.48,\end{array}$$(1) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=2.05\left(9\right)\cdot \log \left(\frac{T_{\text{eff}}}{1\text{\hspace{0.17em}\hspace{0.17em}K}}\right)-3.88\left(33\right),& r=0.50,\end{array}$$(2)

in the “large planet” cluster and $$\begin{array}{cc}\log \left(\frac{M_{\text{P}}}{M_{\text{J}}}\right)=2.04\left(25\right)\cdot \log \left(\frac{T_{\text{eff}}}{1\text{\hspace{0.17em}\hspace{0.17em}K}}\right)-9.09\left(91\right),& r=0.27,\end{array}$$(3) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=0.85\left(12\right)\cdot \log \left(\frac{T_{\text{eff}}}{1\text{\hspace{0.17em}\hspace{0.17em}K}}\right)-3.83\left(46\right),& r=0.43,\end{array}$$(4)

in the “small planets” cluster. In the expressions, the ambiguity (standard deviation) of the last digits are shown in parentheses after the value of each coefficients.

The expressions involving the stellar mass are: $$\begin{array}{cc}\log \left(\frac{M_{\text{P}}}{M_{\text{J}}}\right)=0.78\left(17\right)\cdot \log \left(\frac{M_{\text{S}}}{M_{\odot }}\right)-0.02\left(2\right),& r=0.55,\end{array}$$(5) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=0.48\left(4\right)\cdot \log \left(\frac{M_{\text{S}}}{M_{\odot }}\right)-0.04\left(1\right),& r=0.46,\end{array}$$(6)

in the “large planet” cluster and $$\begin{array}{cc}\log \left(\frac{M_{\text{P}}}{M_{\text{J}}}\right)=0.94\left(9\right)\cdot \log \left(\frac{M_{\text{S}}}{M_{\odot }}\right)-1.46\left(3\right),& r=0.22,\end{array}$$(7) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=0.37\left(5\right)\cdot \log \left(\frac{M_{\text{S}}}{M_{\odot }}\right)-0.64\left(1\right),& r=0.47,\end{array}$$(8)

in the “small planet” cluster. Such a correlation between planet size and stellar mass (but without clustering) was also noted by Fortney et al. (2007); Wu (2019), and Lozovsky et al. (2021). Wu (2019) found the linear coefficient between log M_P and log M_S to be unity, which is in good agreement with Eq. (7). These findings are also consistent with the qualitative observations from Fig. 6 that the most massive stars from our sample (M_S > 1.5 M_⊙) host the largest and least dense hot Jupiters, those with radii of >1.4 R_J.

For the stellar metallicity dependencies, we get $$\begin{array}{cc}\log \left(\frac{M_{\text{P}}}{M_{\text{J}}}\right)=-0.07\left(11\right)\cdot \left[\text{M}/\text{H}\right]+0.01\left(2\right),& r=-0.03,\end{array}$$(9) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=-0.12\left(3\right)\cdot \left[\text{M}/\text{H}\right]+0.07\left(1\right),& r=-0.18,\end{array}$$(10)

in the “large planet” cluster and $$\begin{array}{cc}\log \left(\frac{M_{\text{P}}}{M_{\text{J}}}\right)=0.94\left(9\right)\cdot \left[\text{M}/\text{H}\right]-1.45\left(3\right),& r=0.29,\end{array}$$(11) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=0.17\left(7\right)\cdot \left[\text{M}/\text{H}\right]-0.70\left(7\right),& r=0.17,\end{array}$$(12)

in the “small planet” cluster. Wu (2019) suggested that there is no correlation between stellar metallicity and planetary mass, which is compatible to the low values of correlation coefficients found in our analysis.

The p value of these multiple correlations are in the range of 10⁻⁴ and 10⁻¹⁶ for all equations – except the ones involving [M/H] (Eqs. (9)–(12)), which have p values between 0.01–0.5, showing marginally significant or even, insignificant correlations. The bimodality of the correlations illustrate that there is no one-to-one relationship between planet sizes and stellar parameters. It can be said, however, that in almost every case (with the exception of R_P–M_S), a stronger correlation is found in the case of large exoplanets. In these linear models, no underlying astrophysical connections are considered, they are based solely on statistics; however, in future works, they can be used to further consider the reasons and processes which shape these dependencies.

3.4 Elemental abundances from APOGEE

The APOGEE stellar parameters were taken from a homogeneous, state-of-the-art quality catalog of stellar data, while the sources of stellar data in the NASA Exoplanet Archive are inhomogeneous, although the latter one includes many more planets, which can increase the significance of the statistical analysis based upon it.

From the APOGEE planet sample (Sect. 2.2), we repeated the fits described in Eqs. (2), (6), (3), and (7) leading to the following regressions: $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=1.21\left(43\right)\cdot \log \left(\frac{T_{\text{eff}}}{1\text{\hspace{0.17em}\hspace{0.17em}K}}\right)-4.65\left(1.60\right),& r=0.25\end{array}$$(13) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=0.03\left(5\right)\cdot \left[\text{M}/\text{H}\right]-0.11\left(3\right),& r=-0.05\end{array}$$(14)

in the “large planet” cluster and $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=0.61\left(14\right)\cdot \log \left(\frac{T_{eff}}{1\text{\hspace{0.17em}\hspace{0.17em}K}}\right)-3.05\left(53\right),& r=0.18\end{array}$$(15) $$\begin{array}{cc}\log \left(\frac{R_{\text{P}}}{R_{\text{J}}}\right)=0.02\left(5\right)\cdot \left[\text{M}/\text{H}\right]-0.79\left(1\right),& r=0.02\end{array}$$(16)

in the “small planet” cluster (these clusters and the linear models are shown in Fig. 10). Here, we can see that the APOGEE data reproduced the previously determined coefficients of T_eff in both clusters, with larger ambiguities due to the significantly fewer planets in the currently available APOGEE sample. Also, the r regression coefficients are smaller. In the case of the R–T_eff correlations, we see no significant dependence in the APOGEE data (the ambiguity of the coefficient of T_eff is larger than the value, hence, 0 is within the range) and the r regression coefficients are also close to zero. Here, the APOGEE analysis reproduced the small correlation coefficients between R_P and [M/H]. We note that as the uncertainties in the APOGEE data were also derived in a homogeneous way, we used these as weights for our regression.

The APOGEE data also contains the derived abundances⁴ for a number of elements. We also repeated the linear modeling described above for individual elements inculding C, C I (unionized C), N, O, Na, Mg, Al, Si, K, Ca, Mn, Co, Ni, Ce, and [C/O] to check for correlations between the relative abundance ratio of these elements and the planetary radius (Table 1). In both clusters of “larger” and “smaller” planets, the relationship between the planetary radius and any given X abundance is characterized by the slope of the line (k) and the intercept with the ordinate (n), such that: $$\log \frac{R_{\text{P}}}{R_{\text{J}}}=k\log \frac{X}{X_{\odot }}+n.$$(17)

To show an example of the fitted distributions, the left panel of Fig. 11 shows the exoplanets in the R_P–P plane with the metallicity as the coloring parameter taken from the APOGEE planet host sample. The right panel shows the residuals after removing the bilinear trend according to Eq. (18), which is very similar to the uncorrected one. The comparison of the two panels show that the significant features which can be partly explained by the metallicity includes a negative metallicity gradient along the period (at the boundary, stars with a metallicity higher than 0 are overabundant), mostly affecting the group of smaller planets, and a gradient along the radius. The large end of “smaller planets” and hot Jupiters tend to have increased metallicity. The appropriate KS-test here did not reveal differences of the metallicity distributions near and far from the edge of the desert. The plots for individual elements reproduced a very similar pattern, which we discuss below.

Table 1 shows those elements which we have found to be significantly correlated with the planetary radius within either the “larger planet” or the “smaller planet” sample. As a conclusion, O, Na, and K can influence the size distribution within the smaller planet group, while marginal significance was found in the case of N and Si in the “smaller planet” group as well. We found that only Mg influences the radius of planets in the “large” groups.

Correlations extending to the group of all planets, including both the “larger planet” and the “smaller planet” samples are fitted with a bilinear regression in the period–radius plane. This equation is expressed as follows: $$\log \frac{X}{X_{\odot }}=a\cdot \log P+b\cdot \log \frac{R_{\text{P}}}{R_{\text{J}}}+c,$$(18)

where a, b, and c are the coefficients of the linear regression. Table 2 shows those elements where a significant (p < 0.05) correlation was found. The most significant chemical parameter in the period-radius plane is the [M/H] itself, being very significant in both coordinates and we consider the overal metallicity as the most important forming parameter of the desert as well. The individual abundances, on the other hand, give us more insights to how the overall metallicity determines the chemical dependence of the desert boundaries.

The elements with significant effect of the period-radius distribution, besides what is explained by the [Fe/H], include C, N, O, Mg, Al, Si, K, and Mn. In addition, we did not find significant dependence involving Co and Ce (although these elements were determined only for 14 exoplanet-host stars), as well as [Ca/Fe], [Ni/Fe], and [C/O].

A comparison of Tables 1 and 2 shows that volatiles from the CNO process (especially N and O), moderately refractory alpha-process elements (Mg and Si) and K influence both the global structure of the planetary distribution – surrounding the desert as well – and they have positive correlations with the radius of smaller planets. Stars with increased N, O, Na, and Si, or decreased K abundances tend to form larger planets is the “smaller planet” groups. In the “large planet” group, only Mg is inversely correlated to the planet sizes. This revealed a set of elements that have a significant internal influence to the sizes of smaller and larger planets are disjointed, which serves as strong evidence that the formation and evolution of smaller and larger planets follow a chemical bimodality.

Table 2 shows that the difference between the smaller/larger groups of planets and the global period–radius maps is related to mostly those elements that also appeared in Table 1. Additional elements appear in Table 2 as C, Al, and Mn. This also corroborates the interpretation that elements forming in CNO cycle or alpha process in stellar nucleosynthesis have an influence on the period-radius distribution of the planetary system. Interestingly, all refractory elements showed no or insignificant correlations in both tests, which elements (having the highest condensation temperatures, T_cond > 1500 K) are often considered as the initially condensed material in the protoplanetary disks (Scott & Sanders 2009), and the condensation center for less refractory and volatile materials in further steps. The distribution of planets in the period-radius plane and around the desert seems to be related to more or less volatile elements in the atmosphere, rather than the refractory elements condensed in the cores. Therefore, the desert itself appears to be a predominantly atmospheric feature, rather than a “core” feature, confirming the significance of atmospheric evaporation in its formation.

4 Discussion

Based on the data presented in Sects. 3.2.1–3.2.5, we could generally confirm the results of Szabó & Kálmán (2019). Most importantly, the dependence of the desert boundary on the stellar effective temperature, which had been interpreted as an evidence for the photoevaporation, have been confirmed in the present study as well. An important new detail to this early interpretation is the detection of a multiple and distinct population of hot Jupiters, having inflated and normal hot Jupiters at the boundary in the R_P–P and the M_P–P planes, respectively. Inflated hot Jupiters on the border can be a naturally seen as a new piece of evidence for the scenarios invoking photoevaporation (e.g., Owen & Lai 2018), while the not inflated Jupiters at the boundary in the R_P–P plane suggests a multiple track leading to the formation of the desert.

We also found that the previously claimed metallicity dependence of the border are in general less significant than expected previously. The earlier interpretation of the claimed detection was also interpreted in connection to the evaporation scenarios, as the more effective cooling of a metal-rich atmosphere was claimed to form more irradiation-endurant exoplanets. Our new results have demonstrated that planetary occurrence indeed reveals a dependence between the period, radius, and metallicity (in an agreement with the results of e.g. Demangeon et al. 2018 and Petigura et al. 2018), even though it has not been confirmed as a specific characteristic of the planets near the desert boundaries.

The crowd of the exoplanets with spectroscopically detected atmospheric escape (dos Santos et al. 2021) can also be considered as an evidence of irradiation effects forming the desert, because most of these planets are seen at the boundaries, or at least very close to it (Lecavelier des Etangs, priv. comm.). It can also be pointed out from Fig. 4 that the hottest stars in our sample (with T_eff ~ 8000–10000 K) tend to host the largest and most massive hot Jupiters. A general trend that larger planets are more likely to be hosted by hotter stars can be observed from the median T_eff values from A, B, C, and D regions. Similar statements can be made upon examination of Figs. 5 and 6; namely that stars with larger radius or mass are more likely to host larger planets. Lozovsky et al. (2021) suggested that larger radii of planets surrounding larger stars are a result of different compositions: larger stars tend to host planets with larger H-He mass fractions. Stars with higher metallicities are likely to host larger planets, extending to periods as long as ≈20 days at least (Owen & Murray-Clay 2018). The explanation of the occurrence rate includes the cores of planets around more metal-rich stars to be more massive, hence, they can collect more initial atmospheric mass; and the photoevaporation of a more metal-rich atmosphere is more resistant to stellar irradiation due to the increased cooling related to the photoevaporation of such an atmosphere.

In Lecavelier Des Etangs (2007), additional important factors have been introduced to the formalism, such as atmospheric loss due to tidal forces and the effect of the inclination. These two factors together have an influence on the lifetime of an evaporating atmosphere, leading to the conclusion that in the case of higher inclination values, increased density is required to retain the atmosphere during a specified lifetime. This result suggests that stars with higher temperatures (which tend to host more planets on high inclinations, Winn et al. 2010; Albrecht et al. 2012; Winn & Fabrycky 2015; Zhou et al. 2019) can host only gaseous planets with a higher average atmospheric density for a longer time. This conclusion would at least partly explain the dependence of the desert boundary on the stellar effective temperature, and the general dependence of planetary occurrence on metallicity.

The most important result of our analysis is the detection of multiple relations between M_P and R_P on one side and T_eff and M_S on the other side. These detections reflect the structural difference between large, atmosphere-dominated planets and the smaller Neptunes and super-Earths. These multiple relations complicate the picture at the desert boundary because it is precisely these two kinds of planets that meet there. There are different kind of planets at the upper boundary and the lower boundary, which can naturally lead to different mass and radius dependencies on the stellar parameters. However, in the case of the metallicity dependencies, the sign of the slope is also different and we see a negative correlation between R_P and [M/H] for large planets, along with a positive correlation for small planets (Eqs. (6) and (8)). Therefore, we consider such multiple relations as an evidence for different processes forming the upper and lower boundaries of the desert (while there are at least two different processes at the upper boundary itself, as well; confirming the results from Ionov et al. 2018).

A widely claimed scenario here is the high-eccentricity migration (Owen & Lai 2018), as well as the tidal disruption of giant planets (Matsakos & Königl 2016) or an XUV photoevaporation affecting sub-Saturns (Hallatt & Lee 2022). Deciding which one is the dominant at the lower boundary of the desert requires detailed studies of individual planets. A possible detection of the “sculpted Saturns” scenario would corroborate the prediction of the possible surviving super low-density sub-Saturns to the present day if they are born with even larger atmospheres than they currently harbor – in particular, Kepler 223 d has been claimed to be an example for such a planet on a 14.79 days period orbit (Hallatt & Lee 2022).

There are several other questions that we were not able to answer directly from the present data – the most important being how the evaporation processes and the claimed high-eccentricity migration processes actually form the sub-Jovian/Neptunian desert. Also, because we found that the desert boundaries depend on stellar parameters, the simple linear and log-linear laws that mark the boundary of the desert (Mazeh et al. 2016) are in some way transited to a multilinear dependence. Here, we find a lack of predictions about the desert boundary in the R–P and M–P planes as a function of stellar temperature and metallicity, for instance, and a comparison is not readily possible. Instead, our aim in the present analysis is to show that the formation of the desert is indeed a complex astrophysical process itself. Thus, presenting planetary evolution processes with respect to fundamental stellar parameters should indeed be the next step toward improving our understanding of this process.

Our current conclusions were based on the currently available NASA Exoplanet catalog, combined with the SDSS-IV APOGEE-2 catalog of exoplanet host stars with a limited number of entries (<700 planets in case of all elements). In the future, an extended stellar catalog from the SDSS-V Milky Way Mapper (MWM, Kollmeier et al. 2017) program will cover a much larger sample of planet-host stars, and a similar analysis involving these data will reveal deeper details of these correlations.

5 Summary

In this paper, we analyzed exoplanet occurrence with respect to stellar and planetary parameters. Our main results are as follows:

• We demonstrated the multifaced nature of the upper boundary of the sub-Jovian or Neptunian desert: in the M_P–P plane, the boundary is marked by inflated hot Jupiters; while in the R_P–P plane, there are normal hot Jupiters at the boundary.

• We confirmed the dependence of the period boundary on stellar parameters, such as effective temperature, stellar mass, and stellar radius.

• The multiple populations and the parameter dependence of the boundary suggests multiple formation mechanisms.

• With fuzzy clustering, we also investigated double parameter relations between planet mass and radius on one side, and effective temperature, stellar mass, and metallicity on the other.

• The demarcation of these planet groups coincides with the position of the desert, suggesting that all the relationships connecting the planetary size and the fundamental stellar parameters conjoin at the level of the size ranges of the desert and play a role in its formation.

• In light of these results, we considered photoevaporation as the main process shaping the desert boundaries. We have discussed other possibilities such as high-eccentricity migration, abrasion by irradiation, tidal loss of atmosphere, and internal structural differences.

Acknowledgements

We acknowledge the support of the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, a PRODEX Experiment Agreement No. 4000137122, the Lendület LP2018-7/2022 grant of the Hungarian Academy of Science and the support of the city of Szombathely. LBo acknowledges the funding support from Italian Space Agency (ASI) regulated by “Accordo ASI-INAF n. 2013-016-R.0 del 9 luglio 2013 e integrazione del 9 luglio 2015”. Prepared with the professional support of the Doctoral Student Scholarship Program of the Co-operative Doctoral Program of the Ministry of Innovation and Technology financed from the National Research, Development and Innovation Fund.

References

All Tables

All Figures

	Fig. 1 Sketch of possible paths for already formed planets leaving the sub-Jovian desert (red arrows) and planets not allowed to reach the desert (blue arrows) overplotted to the distribution of known exoplanets. Theories and references behind the scenarios are summarized in the text.
In the text

Fig. 2 Sub-Jovian/Neptune desert of exoplanets in the RP–P and MP–P parameter spaces (top and bottom panel), colored by the planetary density. In the top panel, the upper boundary is marked by “normal” Jupiters with higher density (>1ρJ), while in the bottom panel, the population of inflated planets (red points) appear below the normal Jupiters and, thus, the inflated Jupiters mark the actual boundary.

In the text

	Fig. 3 Scatterplots of the distribution of exoplanets in the RP–P and MP–P planes as in Fig. 2 but with synthetic densities subtracted from the observed ones (see text for details).
In the text

	Fig. 4 Scatterplot showing the distribution of exoplanets in the RP–P (top left) and MP–P (top right) parameter spaces, colored by the temperature of the host star (see the color bar for the values). The middle and bottom panels show the cumulative distribution of the two samples divided at the median of Teff in the A, B, C, and D regions. The p-values of the KS-test are also shown.
In the text

	Fig. 5 Scatterplots showing the distribution of exoplanets in the RP–P and MP–P parameter spaces with Rs used as the third parameter (top row) and the respective cumulative distributions (middle and bottom rows; same as Fig. 4).
In the text

	Fig. 6 Scatterplots showing the distribution of exoplanets in the RP–P and MP–P parameter spaces with MS used as the third parameter (top row) and the respective cumulative distributions (middle and bottom rows; same as Fig. 4).
In the text

	Fig. 7 Same as Fig. 4 but with log ɡ used as the third parameter instead of Teff.
In the text

	Fig. 8 Same as Fig. 4 but with [M/H] used as the third parameter instead of Teff.
In the text

Fig. 9 Correlation between planet mass (top row), planetary radius (bottom row), and the main stellar parameters: Teff (left column), MS (middle column), and [M/H] (right column). The samples are divided into two groups via a fuzzy clustering algorithm applied to the distributions in the right panels. Blue points: Cluster of “larger planets.” Red points: Cluster of ‘smaller planets.” The fitted linear trends of Eqs. (2)–(11) are plotted with solid blue and red lines. Pearson’s r values, corresponding to the giant and smaller exoplanets are also displayed with blue and red, respectively.

In the text

	Fig. 10 Correlations between RP–Teff (left panel) and RP–[M/H] (same as the lower left and lower right panels of Fig. 9) showing the fuzzy clustering of APOGEE data and with the linear regressions of Eqs. (14)–(16) overplotted.
In the text

	Fig. 11 Distribution of exoplanets from the APOGEE planet host sample with [M/H] as the coloring parameter (left panel) and the residuals after a bilinear fit (right panel).
In the text