Sub-Jovian desert of exoplanets at its boundaries. Parameter dependence along the main sequence

Context. The lack of sub-Jovian planets on orbits of $P_{\rm orb}<3$ days is a puzzling aspect of galaxy formation with regard to the distribution of exoplanets whose origins are currently unresolved. Aims. The possible explanations behind the formation of the sub-Jovian or Neptunian desert include several scenarios that can lead to different shapes for the boundary, predicting various dependencies between the position of the boundary and the stellar parameters. Methods. We explored the exoplanet distribution in various 2D and 3D projections, revealing the stellar-dependent substructures in the $P_{\rm orb}-M_{P}$ and the $P_{\rm orb}-R_{P}$ parameter plane. Results. We demonstrate that the upper boundary includes a range of planets, namely, inflated hot Jupiters and normal hot Jupiters, in the two parameter planes, respectively. We confirm the dependence of the boundary on several stellar parameters and, based on a fuzzy clustering analysis, we provide quantitative formulae for the dependencies in groups of smaller and larger planets. The overall period-radius distribution shows chemical substructures as well, with the boundary being dependent on volatiles and alpha-elements, alongside marginal (to none) dependence found for refractory elements. Conclusions. These findings confirm multiple plausible causes for the formation of the desert, particularly preferring those scenarios related to the irradiation-driven loss of the atmospheres of moderately massive planets as the predominant process in shaping planetary distributions.

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 et al. 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 volatilerich 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 P ∝ P − 1 3 , and at the lower edge, R P ∝ P 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.

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 1 ), 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 < 13M 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.16R J is labeled B. In the M P -P plane, we defined the C region between 0.032 − 0.305 M J and the D region between 0.001 − 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 C 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 lowtemperature 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 C are in the desert and B 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 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-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 ) and focused on observing as many planethost 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-meter Sloan Foundation Telescope (Gunn et al. 2006) and the 2.5-meter 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 15140 Å to 16940 Å. 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 multiepoch visits along with high-precision RV measurements are available (Majewski et al. 2017 hood 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 3 (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).

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 instrumentspecific selection distortion of some specific exoplanet discov- Residual log ρ P /ρ J Fig. 3. Scatterplots of the distribution of exoplanets in the R P -P and M P -P planes as in Fig. 2 but with synthetic densities subtracted from the observed ones (see text for details). ery 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).

Dependence of the boundary on stellar parameters
In this subsection, we analyze the R P -P and M P -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).

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 −5 and 10 −4 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.

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 −4 and 10 −5 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 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.

Stellar mass
The scatter plot of exoplanets with the mass of the host star used as the third parameter is shown in Fig. 6.  (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. Article number, page 5 of 11 A&A proofs: manuscript no. 44846corr

Stellar log g
Exploring the effects of log g on the distribution of exoplanets in the two parameter spaces, we divided the samples at log g = 4.46 in the A region, log g = 4.60 in B, log g = 4.48 in C, and log g = 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.

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 g, 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.

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 Earth-sized planets had a median radius of 0.21 R J (interquartile range of 0.14-0.27). When the clustering was done in the mass-metallicity 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]. Between the planet size and the effective temperature, we get: in the "large planet" cluster and log M P M J = 2.04(25) · log T eff 1 K − 9.09(91), log R P R J = 0.85(12) · log T eff 1 K − 3.83 (46), 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. 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 log R P R J = 0.17(7) · [M/H] − 0.70 (7), 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 −4 and 10 −16 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.

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. in the "large planet" cluster and log 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 availagble 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 4 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: 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: 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.

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, personal communication). 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 irradia-   (2007), additional important factors have been introduced to the formalism, such as atmospheric loss due to tidal forces and the effect of he 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.  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 d 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.

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.