ON THE ORIGIN OF THE SUB-JOVIAN DESERT IN THE ORBITAL-PERIOD–PLANETARY-MASS PLANE

As has been known for some time now, radial velocity surveys of exoplanets exhibit an abrupt drop in the number of hot Jupiters (HJs) for orbital periods P_orb ≲ 3 days (e.g., Cumming et al. 2008), corresponding to a pileup of such planets near P_orb = 3 days (e.g., Gaudi et al. 2005) and their paucity at shorter periods (e.g., Zucker & Mazeh 2002). (We define HJs as planets that have masses in the range M_p ∼ 0.3–3 M_J, where M_J is Jupiter's mass, and periods P_orb ≲ 10 days.) Planet candidates identified by transit measurements exhibit a similar sharp drop in the number count below P_orb ∼ 3–4 days for objects with radii ${R}_{{\rm{p}}}\gtrsim 4\;{R}_{\oplus }$ (e.g., Howard et al. 2012; Fressin et al. 2013). One interpretation of these results invokes planet–planet scattering events that place one of the interacting planets on a highly eccentric orbit: that planet thereby attains a small pericenter distance, where its orbit is circularized through tidal interaction with the host star (e.g., Ford & Rasio 2006). In this picture, the shortest possible pericenter distance is given by the Roche limit a_R, where the planet starts to be tidally disrupted, and the planet's resulting circular radius must therefore exceed ∼2 a_R (where the factor of 2 follows from conservation of orbital angular momentum). A planet can also be placed on a highly eccentric orbit by a more gradual process such as Kozai (e.g., Wu & Murray 2003; Fabrycky & Tremaine 2007) or secular (e.g., Wu & Lithwick 2011; Petrovich 2015a) migration. An alternative interpretation attributes the scarcity of high-M_p planets on short-period orbits to tidal exchange of angular momentum with the star, which causes close-in giant planets to spiral inward and get ingested by their host (e.g., Jackson et al. 2009; Teitler & Königl 2014).

As data continued to accumulate, more details emerged about the structure of the planet distribution at low values of P_orb. In particular, Szabó & Kiss (2011) considered the (P_orb, M_p) plane and identified a region (defined by P_orb < 2.5 days and M_p between 0.02 and 0.8 M_J) with a pronounced dearth of planets, which they termed the "sub-Jupiter desert."¹ They noted that this region has both an upper boundary, consisting of HJs, and a lower boundary, consisting of close-in super-Earths (SEs), but is devoid of hot Neptunes. Similar identifications were made by Beaugé & Nesvorný (2013) and by Mazeh et al. (2016; with the latter authors approximating the boundary by two intersecting straight lines in the ($\mathrm{log}{P}_{\mathrm{orb}}$, $\mathrm{log}{M}_{{\rm{p}}}$) plane), and it is generally accepted that this feature is not the result of observational bias.² However, the origin of the desert is still being debated (e.g., Szabó & Kiss 2011; Benítez-Llambay et al. 2011; Beaugé & Nesvorný 2013; Kurokawa & Nakamoto 2014; Batygin et al. 2016; Mazeh et al. 2016). In one proposed scenario (Kurokawa & Nakamoto 2014), the desert is attributed to the evaporation of close-in giant planets by the stellar radiation field, followed by Roche-lobe overflow. However, this picture also implies the complete loss of the gaseous envelopes of the less massive SEs, which is inconsistent with the fact that the inferred radii of the planets near the lower boundary of the desert are generally well in excess of $1.6\;{R}_{\oplus }$ and are thus unlikely to be purely rocky (see Rogers 2015). In an alternative scenario (Valsecchi et al. 2014, 2015), HJs arrive at the RL and undergo rapid mass loss through Roche-lobe overflow, with photoevaporation further contributing to the removal of their envelopes. It was proposed that this mechanism could convert HJs into SEs, but it is not obvious that the resulting low-mass planet distribution would be consistent with the observed shape of the desert's lower boundary. Yet another possibility is that the desert arises from the in situ formation of HJs by gas accretion onto SE cores (Boley et al. 2016; Batygin et al. 2016). However, this picture also faces various challenges (e.g., Inamdar & Schlichting 2015).

In this Letter we propose that the boundary of the desert reflects the locus of the innermost circularized orbits of planets that arrived in the vicinity of the host star through planet–planet scattering, Kozai migration (arising from an interaction with either a stellar or a planetary companion), or a secular process long after the protoplanetary disk had dispersed. As was already pointed out by Ford & Rasio (2006), the Roche limit a_R (which is taken to be the distance of closest approach for a planet arriving via either scattering or Kozai migration) is a function of both R_p and M_p (for a given host mass ${M}_{*}$), and it can therefore be expressed (using the appropriate mass–radius relation) as a function of M_p. A similar expression can be written down for the expected distance of closest approach in the secular chaos scenario. Using an empirical R_p(M_p) relation, we demonstrate that this picture can account for the basic shape of the desert's boundary, with its oppositely sloped upper and lower segments. We further show that the detailed structure of the boundary can be qualitatively reproduced when one also takes into account the orbital evolution induced in close-in planets by their tidal interaction with the host star. These results on the behavior of late-arriving close-in planets complement recent inferences about the evolution of an earlier generation of HJs that had migrated though the natal disk (Matsakos & Königl 2015).

Since most of the planets in the vicinity of the desert's boundary are inferred to have low eccentricities (e < 0.1; see Figure 1 below), we focus on scenarios in which planets are placed on high-eccentricity orbits that are circularized by internal tidal dissipation after the planets arrive in the vicinity of the host star. For any given system, we label by t_arr the time at which the planet under consideration arrives at the stellar vicinity with orbital eccentricity e₀. After arrival, its orbit is circularized to a radius a_c (where its orbital period is P_orb) on a timescale ${\tau }_{{\rm{c}}}\propto {Q}_{{\rm{p}}}^{\prime }{P}_{\mathrm{orb}}^{13/3}{M}_{{\rm{p}}}{R}_{{\rm{p}}}^{-5}$, where Q'_p is the planet's tidal quality factor (e.g., Matsumura et al. 2010a). Planets along the lower boundary of the desert have low values of R_p, which has the effect of increasing τ_c. However, the smallest such planets would be mostly rocky and would thus also have low values of Q'_p (e.g., Goldreich & Soter 1966), so their orbits might still be circularized over the planets' typical ages (t_age ≳ 1 Gyr). In this work we relate R_p and M_p through the empirical relation obtained by Weiss et al. (2013). They divided their sample into small and large planets (superscripts S and L, respectively), separated at M_p ≈ 150 M_⊕, and wrote down power-law relationships that we further approximate here by adopting the median incident flux of the data set (F = 8.6 × 10⁸ erg s⁻¹ cm⁻²) and rounding to the first decimal figure in each exponent:

$$\begin{eqnarray}&&\displaystyle \frac{{R}_{{\rm{p}}}^{{\rm{S}}}}{{R}_{{\rm{J}}}}\sim 1.6\;{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{J}}}}\right)}^{1/2},\quad \quad \displaystyle \frac{{R}_{{\rm{p}}}^{{\rm{L}}}}{{R}_{{\rm{J}}}}\sim 1.5\;{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{J}}}}\right)}^{0}.\end{eqnarray} \tag{ 1 }$$

Zoom In Zoom Out Reset image size

The tidal interaction with the star can also lead to orbital decay and eventual ingestion by the host. This occurs on a nominal timescale ${\tau }_{{\rm{d}}}\sim [{({Q}_{*}^{\prime }/{Q}_{{\rm{p}}}^{\prime })({M}_{*}/{M}_{{\rm{p}}})}^{2}{({R}_{{\rm{p}}}/{R}_{*})}^{5}]\;{\tau }_{{\rm{c}}}$ (where an asterisk again denotes the star), which is typically ≫τ_c (e.g., Matsumura et al. 2010a). Given that ${\tau }_{{\rm{d}}}\propto 1/{M}_{{\rm{p}}}$, the planets on the lower boundary of the desert are hardly affected by this interaction; however, the initial period distribution of the circularized orbits of large planets is measurably modified.

2.1. Initial Locus of Innermost Planets with Circularized Orbits

The post-circularization locus in the (P_orb, M_p) plane of the lowest-P_orb planets can be obtained for both the planet–planet scattering and Kozai migration cases from the expression for the Roche limit,

$$\begin{eqnarray}&&{a}_{{\rm{R}}}=q\;{({M}_{*}/{M}_{{\rm{p}}})}^{1/3}\;{R}_{{\rm{p}}},\end{eqnarray} \tag{ 2 }$$

where the precise value of the coefficient q depends on the planet's structural and orbital characteristics. Two of the most frequently cited estimates are q = 2.16 (Paczyński 1971) and q = 2.7 (Guillochon et al. 2011), but there are indications that its value in real systems could be higher (up to ∼3.6–3.8; Valsecchi & Rasio 2014; Petrovich 2015b). As it turns out, the interpretation of the sub-Jovian desert in terms of the Roche limit also favors a comparatively large value for q (≃3.5; see Section 3). Under the assumption that e₀ ≈ 1, the circularization radius is ≃2 a_R and thus can be written as

$$\begin{eqnarray}{a}_{{\rm{c}},\mathrm{RL}} & \approx & 0.03\;\left(\displaystyle \frac{q}{3}\right){\left(\displaystyle \frac{{M}_{*}}{{M}_{\odot }}\right)}^{1/3}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{J}}}}\right)}^{-1/3}\times \left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{{\rm{J}}}}\right)\;\mathrm{au}.\end{eqnarray} \tag{ 3 }$$

In the secular chaos (SC) model, the circularization radius was estimated by Wu & Lithwick (2011) as twice the pericenter distance obtained by equating the precession rate of a planet's longitude of pericenter—due to a secular interaction with another planet located farther out—with the orbit-averaged precession rate due to the tidal quadrupole induced on the planet by the star. This yields

$$\begin{eqnarray}{a}_{{\rm{c}},\mathrm{SC}} & \approx & 0.03{\left(\displaystyle \frac{\alpha }{1/6}\right)}^{-3/5}{\left(\displaystyle \frac{{M}_{*}}{{M}_{\odot }}\right)}^{2/5}{\left(\displaystyle \frac{{M}_{\mathrm{pert}}}{{M}_{{\rm{J}}}}\right)}^{-1/5}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{{\rm{J}}}}\right)}^{-1/5}\left(\displaystyle \frac{{R}_{{\rm{p}}}}{{R}_{{\rm{J}}}}\right)\;\mathrm{au},\end{eqnarray} \tag{ 4 }$$

where M_pert is the mass of the perturbing planet and α is the ratio of the semimajor axes of the two planets. For the sake of simplicity, we henceforth fix α at its normalization value and treat M_pert as the relevant model parameter.³

By combining Equations (1), (3), and (4), and using ${P}_{\mathrm{orb}}\propto {a}_{{\rm{c}}}^{-3/2}$, we obtain the initial (post-circularization) locus of the innermost planets for these two scenarios,

$$\begin{eqnarray}&&\mathrm{RL}:\left\{\begin{array}{l}{P}_{\mathrm{orb},\mathrm{RL}}^{{\rm{L}}}\propto {M}_{{\rm{p}}}^{-1/2}\\ {P}_{\mathrm{orb},\mathrm{RL}}^{{\rm{S}}}\propto {M}_{{\rm{p}}}^{1/4}\end{array}\right.,\quad \mathrm{SC}:\left\{\begin{array}{l}{P}_{\mathrm{orb},\mathrm{SC}}^{{\rm{L}}}\propto {M}_{{\rm{p}}}^{-3/5}\\ {P}_{\mathrm{orb},\mathrm{SC}}^{{\rm{S}}}\propto {M}_{{\rm{p}}}^{3/20}\end{array}\right..\end{eqnarray} \tag{ 5 }$$

It is seen that in both cases the inferred boundary has a negatively sloped upper branch and a positively sloped lower branch. This basic property of the observed sub-Jovian desert is a generic feature of our proposed interpretation, arising from the different slopes in the empirical R_p(M_p) relation for large and small planets.

2.2. Effect of Tidally Induced Orbital Evolution

We calculate this effect following the procedure outlined in Matsakos & Königl (2015). We perform Monte Carlo simulations for a population of 30,000 planetary systems, assuming solar-type hosts (${M}_{*}={M}_{\odot }$, ${R}_{*}={R}_{\odot }$) with initial rotation periods distributed uniformly in the range 5–10 days, and drawing ages in the range 1–8 Gyr from the empirical distribution of Walkowicz & Basri (2013). We also adopt Q'_* = 10⁶ (see Section 3). The values of P_orb and R_p are chosen from the observationally inferred distributions presented in Youdin (2011), with M_p deduced from the ${R}_{{\rm{p}}}({M}_{{\rm{p}}})$ compilation of Weiss et al. (2013; using their inferred power-law fit for small planets but accounting for the observed scatter in both radius and mass for those with ${M}_{{\rm{p}}}\gtrsim 150\;{M}_{\oplus }$).⁴ We further assume a random distribution for the initial angle between the stellar spin and the orbital plane, and that the distribution of planet arrival times is uniform in $\mathrm{log}{t}_{\mathrm{arr}}$ (see Section 4). These choices are summarized in Table 1.

Table 1.  Model Parameters

Parameter	Sampled Value	Distribution
Number of planets	5	⋯
${P}_{\mathrm{orb}}\;(\mathrm{days})$	0.5–50 (1)	$f(\mathrm{ln}P)\propto {P}^{0.47}$
${R}_{{\rm{p}}}\;({R}_{\oplus })$ for ${R}_{{\rm{p}}}\lt 12\;{R}_{\oplus }$	3–12 (1)	$f(\mathrm{ln}R)\propto {R}^{-0.66}$
${M}_{{\rm{p}}}\;({M}_{\oplus })\;\mathrm{for}\;{R}_{{\rm{p}}}\lt 12\;{R}_{\oplus }$	9–144 (2)	${({R}_{{\rm{p}}}/{R}_{\oplus })}^{2}\;{M}_{\oplus }$
${R}_{{\rm{p}}}\;({M}_{\oplus })$ for ${R}_{{\rm{p}}}\gt 12\;{R}_{\oplus }$	9–20 (2)	Fit to data
${M}_{{\rm{p}}}\;({M}_{\oplus })\;\mathrm{for}\;{R}_{{\rm{p}}}\gt 12\;{R}_{\oplus }$	56–5620 (2)	Fit to data
${t}_{\mathrm{arr}}$ (Gyr)	0.01–10	Uniform in log time

References. (1) Youdin (2011), (2) Weiss et al. (2013).

Download table as:  ASCIITypeset image

For each system, we integrate the evolution equations for the stellar and orbital angular momenta, taking into account the effects of equilibrium tides and of magnetic braking. We assume that all orbits in a multi-planet system remain coplanar, but we neglect planet–planet interactions. We carry out the calculations for each of the two distributions obtained in Section 2.1, removing from consideration any planet drawn with an initial semimajor axis that is less than ${a}_{{\rm{c}},\mathrm{RL}}$ (Equation (3)) or ${a}_{{\rm{c}},\mathrm{SC}}$ (Equation (4)). Given that we do not include systems with ${t}_{\mathrm{age}}\leqslant 1\;{\rm{Gyr}}$ in our population counts, our predicted distributions do not reveal the possible presence of close-in giant planets that arrive at the stellar vicinity by migration through the natal disk. The contribution of such early-arriving planets to the observed number count is, however, small due to their relatively rapid tidal ingestion by the star. Nevertheless, as was noted by Matsakos & Königl (2015), these HJs could have a strong influence on the observed distribution of the angle between the stellar spin and the orbital plane. Based on that work, we take account of this effect by including, at time t = 0, a "stranded HJ" (SHJ) characterized by M_SHJ = 0.6 M_J, R_SHJ = R_J, and P_SHJ = 2 days in (randomly selected) 50% of the modeled systems. A giant planet of this type would undergo tidal ingestion by a G-type star on a timescale of ∼0.7 Gyr.

Figure 1 presents the predicted distributions in the (P_orb, $\mathrm{log}{M}_{{\rm{p}}}$) plane for the RL and SC models. For the RL model we also plot the relation obtained from combining Equations (1) and (3) for three values of the parameter q,⁵ whereas for the SC model we show the corresponding relation obtained from Equations (1) and (4) for three values of M_pert.⁶ It is seen that each of these scenarios can reproduce the observed shape of the sub-Jovian desert's boundary quite well. The initial circularization locus accounts for the basic "bird's beak" configuration of the boundary, but only the shape of its low-mass segment (along which orbital evolution effects are negligible) and the location of the vertex (where the upper and lower segments intersect) are preserved in the final distribution. The mass coordinate of the vertex is fixed by the transition point (M_p ≈ 150 M_⊕) of the empirical mass–radius relation and is independent of the details of the underlying planet migration model; the fact that it matches the observations so accurately is a strong indication that this aspect of the small/large planet dichotomy plays a key role in determining the desert's shape. On the other hand, the plots in the left panels of Figure 1 demonstrate that the period coordinate of the vertex is sensitive to the value of the relevant model parameter. In each case, however, the inferred best-fit value has a plausible magnitude. In the RL picture, the comparatively large indicated value of q (≃3.5) is consistent with other recent determinations (e.g., Valsecchi & Rasio 2014, who derived the initial locations of HJs detected inside the Roche limit; and Petrovich 2015b, who investigated the formation of HJs by Kozai migration). In the case of secular chaos, the inference that ${M}_{\mathrm{pert}}\approx {M}_{{\rm{p}}}$ is also eminently reasonable. For Jupiter-mass planets, the excitation of secular interactions would be much weaker if the companion planet's mass were measurably less than that of the innermost planet (e.g., Wu & Lithwick 2011), whereas for SE planets, comparable-mass objects are the most frequent companions and are thus the most likely to act as perturbers. Note that the value of the model parameter (q or M_pert) also affects the predicted shape of the lower branch of the desert's boundary, which further constrains the choice of the best match.

The shape of the upper boundary is largely independent of the details of the initial distribution of high-mass planets, and is determined mainly by the ensuing orbital evolution. It is well reproduced—in the context of the equilibrium tidal interaction model that we employ—with a stellar tidal quality factor ${Q}_{*}^{\prime }\sim {10}^{6}$, which is consistent with previous inferences from modeling the P_orb distribution of HJs (e.g., Teitler & Königl 2014; Essick & Weinberg 2016). It is noteworthy that the two alternative models originally proposed to explain the spatial distribution of HJs—circularization of highly eccentric orbits and tidal exchange of angular momentum with the star—are both found to be relevant to the interpretation of the detailed shape of the desert's boundary in our picture.

Our proposed interpretation of the sub-Jovian desert's shape—that it is a natural consequence of the orbital circularization of planets that arrive at the stellar vicinity on high-eccentricity orbits and of their subsequent tidal angular-momentum exchange with the star—is consistent with the finding that the eccentricity distribution of giant planets broadens with increasing period (e.g., Winn & Fabrycky 2015) and with the evidence for ongoing orbital evolution of the innermost HJs (e.g., Valsecchi & Rasio 2014; Teitler & Königl 2014; Essick & Weinberg 2016). The latter result, in turn, supports the view that the observed HJs are mostly late-arriving planets rather than the product of migration in the protoplanetary disk. The high-eccentricity migration scenario, in which the ingoing planet is placed on a high-e₀ orbit through gravitational interaction with one or more massive bodies (other planets or a binary star), is consistent with the indicated high occurrence rate of planetary and stellar companions in systems that harbor HJs (e.g., Knutson et al. 2014; Ngo et al. 2015). We considered three possible pathways to such an outcome—planet–planet scattering, Kozai migration, and secular drift—and inferred that in principle they could all play a role. However, a more detailed examination is required to determine the actual contribution of each of these processes. One pertinent question is whether the influence of the process extends to late times; for example, secular chaos is inherently a long-lasting interaction, for which the expected distribution of arrival times is approximately uniform in $\mathrm{log}{t}_{\mathrm{arr}}$ (Y. Lithwick 2015, personal communication), whereas scattering likely only plays a role at early times (≲10⁸ years; e.g., Chatterjee et al. 2008; Jurić & Tremaine 2008). Another relevant question concerns the rate of eccentricity growth in the high-e migration model; for example, Wu & Lithwick (2011) argued that the observed ∼3 day pileup of HJs is more likely to have been caused by a slow process such as Kozai or secular migration than by sudden scattering events (see also Nagasawa et al. 2008). Further studies are needed to fully address this issue.

One attractive feature of the proposed interpretation is that it can simultaneously reproduce both the upper and the lower segments of the desert's boundary. An alternative possibility for the low-mass branch is that it arises from the in situ formation of close-in SEs (Mazeh et al. 2016), which has been recently explored in the literature (e.g., Lee & Chiang 2016). However, arguments have also been given in favor of formation at larger distances (up to a few au) and subsequent inward migration (e.g., Schlichting 2014; Inamdar & Schlichting 2015). Our model for the sub-Jovian desert is consistent with the latter picture.

The scenario considered in this work, which involves late-arriving planets, complements the SHJ model presented in Matsakos & Königl (2015), which is concerned with early-arriving HJs. That model explains the good alignment exhibited by a significant fraction of HJs around cool (G-type) stars—as well as the good alignment inferred for more distant planets around such stars—versus the broad range of obliquities exhibited by HJs around hot (F-type) stars. As we noted in Section 2, the earlier population of HJs would not affect the observed properties of the sub-Jovian desert because of the relatively short SHJ ingestion time. However, consistency with the SHJ scenario requires the orbital plane of any late-arriving planet to roughly coincide with that of the natal disk. This additional constraint on the high-eccentricity migration model need not, however, be too difficult to fulfill. For example, Matsumura et al. (2010b) found that only ∼15% of planets in a three-planet system that emerges out of a gas disk have orbital inclinations >10°, and Lithwick & Wu (2014) discovered that, even if larger initial inclinations are allowed, 60% of the HJs formed through secular migration involving a three-planet system have projected obliquities <10°. Furthermore, Petrovich (2015a) demonstrated the feasibility of producing HJs through coplanar high-e migration (a secular process involving two planets).

We are grateful to Tsevi Mazeh for alerting us to the "bird's beak" shape of the sub-Jovian desert's boundary. We thank him as well as Dan Fabrycky, Yoram Lithwick, Leslie Rogers, and the referee for helpful input. This work was supported in part by NASA ATP grant NNX13AH56G.