OBSERVATIONAL EVIDENCE FOR TIDAL DESTRUCTION OF EXOPLANETS

The distribution of the orbits of close-in exoplanets provides important constraints on models of planet formation and evolution. These planets are unlikely to have formed in their current orbits (with semimajor axis a < 0.1 AU) because the protoplanetary gas disk from which they accreted was probably too warm so close to the host star (Lin et al. 1996). For example, the core-accretion model requires coagulation of solid material into a planetary core (Pollack et al. 1996). Within a few 0.01 AU of a star, however, temperatures are too high for condensation of solid material. Instead, the exoplanets observed today in close-in orbits likely formed several AU from their host stars, where temperatures were low enough for condensation of solid material, and later migrated inward (e.g., Lin & Papaloizou 1985). An important process driving the migration to ∼ 0.1 AU was probably gas drag in the protoplanetary gas nebula (Lin et al. 1996). Once planets are inside 0.1 AU, tides can play an important role (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Barnes et al. 2008; Jackson et al. 2008a, 2008b; Nagasawa et al. 2008). Here we consider the role of tides in shaping the observed orbital distribution.

Figure 1 shows the semimajor axes and estimated ages for many close-in exoplanets. (The data for Figure 1 are given in Table 1.) Most of these planetary systems are younger than 6 Gyr, and a-values are concentrated between 0.03 and 0.7 AU, centered on about 0.05 AU, which corresponds to an orbital period of about 3 days for a star of one solar mass. Older planets tend to be farther away from their host stars, although there is considerable uncertainty regarding stellar ages (e.g., Saffe et al. 2005; Soderblom 2009) and the sample size is limited. Also, transiting planets (the filled squares in Figure 1) tend to be younger than nontransiting planets: all but one transiting planet (XO-5 b) are younger than 6 Gyr, while there are many older nontransiting planets. Perhaps most significant, there are no planets with a < 0.017 AU, and this lower cutoff (approximated by the dashed line in Figure 1) increases with age. For example, among planets younger than 5 Gyr, a-values are as small as 0.017 AU, whereas among planets older than 5 Gyr, a-values are all greater than 0.036 AU. Given the uncertainty in stellar ages, we performed Monte Carlo simulations that incorporated uncertainties in age and a and found that this trend is persistent. As we will show, these trends are all consistent with the evolution and eventual destruction of planets by tides.

Zoom In Zoom Out Reset image size

Previous studies have addressed the apparent concentration of a-values centered at 0.05 AU. This concentration is often called the "three-day pile-up" (e.g., Ford & Rasio 2006; Cumming et al. 2008; Jackson et al. 2008a), a term that reflects an implicit assumption that inward migration deposited these planets at their current orbits, like rocks deposited in a moraine at the front edge of a glacier. Thus, studies of the inward migration have tried to account for a stopping mechanism at around 0.05 AU.

Models in which interactions between planets and a protoplanetary gas disk caused the planets to migrate often invoke a clearing of the gas close to the star to halt the migration (Lin et al. 1996; Trilling et al. 1998; Ward 1997; Papaloizou 2007). These studies have suggested a host star's magnetosphere can clear out the gas disk within a few 0.01 AU of the star, and once the migrating planets enter the cleared region, they stop migrating, resulting in a "pile-up."

Other studies have invoked a combination of gravitational scattering and tidal evolution to move planets inward to the "pile-up" (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Ford & Rasio 2006; Nagasawa et al. 2008). In such scenarios, interactions between planets in a dynamically unstable system scatter one of the planets into a highly eccentric orbit, with a pericenter distance close enough to the star that tides would affect the orbit. Ford & Rasio (2006) showed that, if scattering resulted in e ∼ 1 and pericenter distances close to the Roche limit of the star, tides raised on the planet by the host star would drive a-values to nearly twice the Roche limit as the eccentricity dropped to zero, at which point tidal evolution resulting only from tides raised on the planet would cease. This process would typically deposit a planet close to 0.05 AU.

However, such models have generally ignored the effect of the tide raised on the host star, which plays a significant role in the evolution of the orbits of close-in exoplanets (Trilling et al. 1998; Dobbs-Dixon et al. 2004; Adams & Laughlin 2006; Jackson et al. 2008a; Levrard et al. 2009; Barker & Ogilvie 2009). In this paper, we consider the effect of the tides raised on the star, which can dominate the orbital evolution for a < 0.05 AU. We find that the concentration of observed planets near 0.05 AU may result from inward migration and destruction of closer-in planets, rather than a pile-up of planets that migrated in from farther out. Previous theoretical studies have shown that tides may pull planets into stars in only a few Gyrs (Mardling & Lin 2004; Raymond et al. 2008; Levrard et al. 2009). In this study, we present observational evidence that this destruction has already occurred.

By modeling the tidal evolution of a hypothetical population of close-in planets, we show how the distribution of observed a-values with age provides evidence for tidal destruction. We find that many of the features of the observed distribution agree with our models if the ages of our model planetary systems are distributed uniformly. However, in the solar neighborhood, there are observed to be more younger than older stars, and the agreement between the observed and modeled distributions of a-values is significantly improved if we account for the prevalence of younger planetary systems. We also find that our model can reproduce the observed clustering of a-values if we assume the distribution of initial a-values was weighted to smaller a.

As we demonstrate below, the tendencies for older exoplanets to be farther from their host stars and for transiting planets to be young likely result from tidal evolution of orbits. We also show that the apparent clustering of orbital periods near three days does not require a mechanism to stop the inward migration of close-in exoplanets. Rather the observed cutoff in a-values is a natural result of tidal evolution. Our results indicate that many observed close-in exoplanets will probably be destroyed by tides in a few billion years.

Tides have been considered in many contexts in the study of extrasolar systems. For example, tidal heating may be sufficient to account for many of the discrepancies between the predicted and observed radii of transiting exoplanets (Bodenheimer et al. 2001, 2003; Burrows et al. 2007; Jackson et al. 2008b; Liu et al. 2008; Ibgui & Burrows 2009; N. Miller et al. 2009, in preparation). Tides also reduce orbital eccentricities e, probably leading to the observed dichotomy between the relatively narrow distribution of e-values for close-in exoplanets and the broad distribution of e-values for planets far from their host stars (Rasio et al. 1996; Jackson et al. 2008a).

Tides also change orbital semimajor axes a (Jackson et al. 2008a; Levrard et al. 2009), and thus play a key role in producing the observed distribution of a-values. As we show below, tides can reduce a-values so much that many exoplanets have probably already crossed the critical distance (e.g., the Roche limit) inside of which there would be tidal disruption and accretion of the planet by the host star. As shown below, this tidal destruction probably leads to the observational trends pointed out in Figure 1.

In order to examine how tides can shape the observed distribution (Figure 1), we invoke the classical equations of tidal evolution (cf. Goldreich 1963; Goldreich & Soter 1966; Jackson et al. 2008a; Ferraz-Mello et al. 2008):

$$\begin{eqnarray} \frac{1}{a}\frac{{{\rm d}a}}{{{\rm d}t}} &=& {-} \left({\frac{{63}}{2}(GM_* ^3)^{1/2} \frac{{R_p ^5 }}{{Q^\prime_p M_p }}e^2 + \frac{9}{2}(G/M_*)^{1/2} \frac{{R_* ^5 M_p }}{{Q^\prime_* }}}\right.\nonumber\\ &&\left.\times\left({1 + \frac{{57}}{4}e^2 } \right) \right)a^{ - 13/2} \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\frac{1}{e}\frac{{{\rm d}e}}{{{\rm d}t}} \!=\! {-} \left(\!{\frac{{63}}{4}(GM_* ^3)^{1/2}\! \frac{{R_p ^5 }}{{Q^\prime_p M_p }} + \frac{{225}}{{16}}(G/M_*)^{1/2} \frac{{R_* ^5 M_p }}{{Q^\prime_* }}} \!\right)\!a^{ - 13/2}\nonumber\\ \end{eqnarray} \tag{ 2 }$$

where G is the gravitational constant, R is a body's radius, M its mass, and Q' its modified tidal dissipation parameter (as defined in Goldreich & Soter 1966), and subscripts p and * refer to the planet and star, respectively. Note that Equations (1) and (2) include corrections to those presented by Jackson et al. (2008a). They include an additional term of order e² in Equation (1) due to tides raised on the star and a corrected numerical coefficient in the last term of Equation (2) (e.g., Ferraz-Mello et al. 2008). (Neither of these corrections is great enough to significantly affect the results in Jackson et al. 2008a, 2008b, 2008c, 2008d or Barnes et al. 2008.)

These equations involve several assumptions about tidal processes, many of which are detailed in Jackson et al. (2008a). One that is particularly relevant here is the assumption that the period of rotation for exoplanetary host stars is always significantly longer than the planet's orbital period, an assumption that is probably reasonable for most close-in exoplanets (see Section 4.2). This assumption yields the particular numerical coefficients in the terms containing Q'_* in Equations (1) and (2). Different assumptions about the periods would change the values and signs of these coefficients (Ferraz-Mello et al. 2008).

Several previous studies overestimate the time for a planetary orbit to decay into the star, the so-called "orbital decay timescale." These studies attempted to simplify the tidal evolution equations by calculating the current values of a/(da/dt) and using it to estimate the timescale, essentially assuming a damps exponentially with time. For example, Rasio et al. (1996) estimated that the orbital decay timescale for 51 Peg b was a trillion years, longer than the main sequence lifetime of its host star. However, the strong dependence of da/dt on a (Equation (1)) means the ratio a/(da/dt) is not constant, and therefore a does not damp exponentially. Instead, a numerical solution of the coupled Equations (1) and (2) shows that da/dt rapidly increases in magnitude as a drops. The orbit decays more and more rapidly, as was also pointed out by Jackson et al. (2008a) and Levrard et al. (2009). In fact, 51 Peg b's orbit could decay in only a few billion years, as discussed in more detail below.

Presumably once the planet's pericenter distance p = a(1 – e) is close enough to its host star, the tidal gravity of the star will disrupt the planet, and the planet will be destroyed. Figure 2 shows the coupled tidal evolution of pericenter p and eccentricity e from their current values forward in time for two close-in exoplanets for a range of stellar dissipation parameters Q'_*. The horizontal line near the bottom of each plot is the Roche limit a_R for that system, given by a_R = (R_p/0.462) (M_*/M_p)^1/3 (Faber et al. 2005). The classical Roche limit involves various assumptions about the planet's density, structure, and strength, so the actual distance of tidal break-up may be different depending on physical parameters. However, the exact value for the Roche limit is not important because once a tidally evolving planet gets within a few 0.01 AU of its host star, its orbit evolves so quickly that the timescale for the destruction of the planet is not sensitive to the exact break-up distance.

Zoom In Zoom Out Reset image size

In Figure 2(a), we show the tidal evolution of CoroT-7 b, which is a transiting exoplanet with the smallest radius found to date, about 2 Earth radii (exoplanet.eu), indicating it may be a rocky planet. Its current orbital eccentricity has been reported as zero, so p = a. As a result, only the tide raised on the star by the planet contributes to its tidal evolution (Equation (1)). As p begins to decrease, the rate of decrease accelerates rapidly. For the commonly adopted value Q'_* = 10⁵ (Mathieu 1994; Lin et al. 1996), CoroT-7 b's orbit crosses the Roche limit in a few tens of millions of years, but even for Q'_* = 10⁷, CoroT-7 b would be doomed in the next few billion years.

For 51 Peg b (Figure 2(b)), which was the first exoplanet discovered around a Sun-like star (Mayor & Queloz 1995), recent estimates give e = 0.013 (Butler et al. 2006). Again, as shown in Figure 2(b), tidal destruction is possible within the lifetime of the star. However, for this planet, the nonzero e-value means that the tide raised on the planet contributes slightly to the orbital evolution. As a result, before 51 Peg b plummets to the Roche limit, its eccentricity decreases, and p can rise a bit temporarily. While the pericenter distance does not decrease monotonically, it does eventually accelerate toward the star.

These examples show how dependence of da/dt on a can accelerate the removal of a planet. Consequently, it is the initial a-value that determines the time before a planet's orbit crosses the Roche limit, and not so much the initial eccentricity. The acceleration of tidal evolution with decreasing a means that a population of planets with small initial a-values (orbiting similar stars) would have their a-values spread apart more quickly than a group with large initial a-values. This process must have a strong effect on the density distribution of those planets as a function of a and age, and thus can help shape the observed distribution shown in Figure 1.

3.1. Effects of Decreasing a

We can illustrate this spreading process and its effect on the orbital distribution by considering a hypothetical population of Jupiter-like planets (with M_p = 1 Jupiter mass and R_p = 1 Jupiter radius) in initial orbits around Sun-like stars (with M_* = 1 solar mass, R_* = 1 solar radius and Q'_* = 10⁶) with a range of a-values uniformly distributed between 0.01 AU (near the Roche limit) and 0.2 AU (beyond which tides have little effect, e.g., Jackson et al. 2008a). We set initial e-values to zero, because (as discussed above) they have little effect on the lifetime before the orbits cross the Roche limit. Here we are considering the evolution of many individual planets, each in its own planetary system, rather than many planets in a single system. For each planet, we modeled the tidal evolution forward in time according to Equation (1). (Equation (2) is not relevant in this case of circular orbits.) Figure 3 shows the evolution of the population over 15 Gyr. The closest planets plunge rapidly in toward their stars, while planets farther out hardly evolve. The planets closest to their stars become less densely distributed in a-values.

Zoom In Zoom Out Reset image size

The decrease in the density of planets for smaller a-values is quantified in the histogram in Figure 4, which shows the evolution of number density (number of planets/increment in a) as a function of semimajor axis, and how this function changes with age. The results shown in Figure 4 are based on calculations similar to Figure 3, but with many more planets: 10,000 planets spread initially from a = 0.01 to 0.2 AU. For Figure 4, planets are binned by a-value in bins 0.005 AU wide, and then the number of planets in that bin is normalized by dividing by the total number of planets in the original population, giving the normalized number density. Initially (t = 0), all bins are filled evenly, reflecting the assumed uniform spread in a-values. However, as time passes, the innermost bins are rapidly cleared out, the outermost bins are unaffected, and planets from intermediate bins move inward to smaller a-values. As suggested by Figure 2, the rate at which bins clear out depends on the value of Q'_*, but for any value of Q'_*, the innermost bins clear first.

Zoom In Zoom Out Reset image size

If the initial a-values for real planets were uniformly distributed (and by "initial" we mean when tides began to dominate the evolution), we would expect the distribution of planets to be similar to Figure 4, with relatively few in a-bins with low densities and more in bins with high densities. For example, according to Figure 4, among planets that are 1 Gyr old, we might expect to observe about twice as many with a-values near 0.03 AU as we find with a-values near 0.025 AU.

In order to compare these models of the evolution of number density to the observed distribution, we recast Figure 4 as a plot of semimajor axis versus stellar age, analogous to Figure 1. From Figure 4, we can extract the locus of points with any fixed number density on a plot of semimajor axis versus age, yielding the solid curves in Figure 5. For example, consider a normalized number density of 0.015 in Figure 4. We see that after 1 Gyr, the histogram bin corresponding to a density of 0.015 is centered on 0.035 AU. After 2 Gyr, the bin with that same density is centered on 0.04 AU. In order to map the 0.015 density contour in Figure 5, we draw a curve connecting the point at 1 Gyr and 0.035 AU to the point at 2 Gyr and 0.04 AU (and also to all the intermediate and exterior points). This increase with age in the a-value corresponding to a fixed number density is reflected by the positive slope of the contours in Figure 5.

Zoom In Zoom Out Reset image size

If the population represented by these contours were searched observationally, depending on the total number of planets and the limits of completeness of the search, some bins might not be populated in the observational data. That is, we would not expect to have observed planets in the parts of the distribution below a particular limiting density. Thus the contour line for that density should represent a boundary, below which planets would not appear in the observed population. Note, however, that without knowing the actual number of planets or the completeness of the search to date, we cannot say which specific contour should give that boundary.

However, we can test the plausibility of our model (that the low-a cutoff is due to tidal migration into the star) by determining whether any contour of the model distribution follows the cutoff in the observed distribution. Figure 5 includes the observed distribution and the apparent cutoff (from Figure 1) for comparison with the model contours. We see that, for all three values of Q'_*, the model contours do predict a paucity of planets with small a, and the contours do have slopes toward larger a with greater age. But none of the contours matches the relatively steep slope of the observed cutoff (dashed line in Figures 1 and 5).

3.2. Effects of the Age Distribution

The model above does not include an important selection effect: there are inherently fewer older stars than there are younger stars in the solar neighborhood (Takeda et al. 2007). Figure 6 shows a histogram of stellar ages for a large population of nearby F, G, and K stars, as reported by Takeda et al. (2007). The number of stars in a given age bin declines roughly linearly with age (although stellar ages do remain very uncertain, e.g., Saffe et al. 2005; Soderblom 2009). We can also see this trend in Figure 1, especially at large a-values (where tides have little effect). For this version of our model, we represent this reduction in the number of older stars with a straight-line fit (as shown in Figure 6).We use that function as the relative probability for a star to be a certain age.

Zoom In Zoom Out Reset image size

We apply this reduction in the number of older stars to the population shown in Figure 4, reducing the number of planets (uniformly in all a-bins) in accord with their increasing age. Figure 7 shows the results of modifying simulated population in this way.

Zoom In Zoom Out Reset image size

The density contours corresponding to Figure 7 are shown on a plot of semimajor axis versus age in Figure 8. Compared with Figure 5, the shapes of the contours at low a are steeper and, for great enough age, the contours become vertical, reflecting the reduction in all a-bins of the number of planets with age.

Zoom In Zoom Out Reset image size

Comparing the slopes of the model contours to the low-a cutoff in the observed population, we find good fits for all values of Q'_*. For Q'_* = 10⁵, 10⁶, and 10⁷, respectively, there is a reasonable fit for the density contours with values 0.00015, 0.0015, and 0.0045. In addition, the contours fit the reduction in numbers of old planets for larger a-values. Thus, incorporating the reduction in the number of older stars shows that observations are consistent with destruction of planets through tidal decay of their orbits.

3.3. Effects of the Initial a-distribution

The analysis above assumes that planets begin with a uniform distribution in a. We also investigate the effects of choosing an alternative distribution of initial a-values. For example, consider an initial population of planets whose number density decreases with increasing a, as in the linear distribution at t = 0 in Figure 9. Such a distribution might reflect the effects of gas disk migration, for example. Figure 9 shows the evolution of such a population due to tidal variation of a and the loss with age (according to Figure 6). As time passes, the number density for the innermost bins drops rapidly, and the density in the bins farther out remains mostly unchanged, just as in previous calculations. Other assumptions about the initial distribution give qualitatively similar results.

Zoom In Zoom Out Reset image size

Density contours for semimajor axis versus age are shown in Figure 10 for this case, analogous to Figures 5 and 8. Here contour lines show that the density drops for increasing a where a > 0.1 AU. At the lower-a cutoff, the contours lines follow the slope of the observed distribution (dashed line). In general, the model contours fit the observed distribution reasonably well. The agreement is especially good for Q'_* = 10⁷, where the model yields a concentration of a-values around 0.05 AU, similar to that observed. These results suggest that we can account for the observed clustering of a-values and the low-a cutoff by invoking an initial distribution that favors closer planets, followed by tidal orbital decay. The concentration at 0.05 AU may not be a pile-up, but rather it could be the result of carving away part of the initial population. The orbits of many planets have decayed into their parent stars, and more will do so in the future.

Zoom In Zoom Out Reset image size

4.1. Transiting Planets

4.2. Stellar Rotation Rates

Our results suggest that many planets have already been accreted into their host star. One potentially observable effect would be an increase in stellar spin rate, as a planet's orbital angular momentum is transferred to the star. The change in rotation rate for a star is given by

$$\begin{eqnarray} &&\Delta \Omega _* = \frac{{M_p \sqrt {GM_* a(1 - e^2)} }}{{C_* M_* R_*^2 }} \end{eqnarray} \tag{ 3 }$$

where C_* is the stellar moment of inertia coefficient. For example, if a Jupiter-like planet in a circular orbit (e = 0) with initial a = 0.05 AU plunges into a Sun-like star, with a typical 30 day period and C_* ∼ 0.1 (Massarotti 2008), the rotation period is reduced to only ∼ 5 days. Because main-sequence stars would otherwise tend to spin down as they age (Skumanich 1972; Barnes 2007b), an anomalously high spin rate could be evidence for accretion of a tidally destroyed planet. Some anomalously rapid rotators have been found among red giants, and their rapid rotation rates have been taken as evidence for accretion of planets by the stars (Massarotti 2008). Perhaps candidates could also eventually be identified among main sequence stars.

4.3. Stellar Compositions

4.4. Stripped Planetary Cores

4.5. Uncertainties and Assumptions

Tidal decay of close-in exoplanetary orbits can, and probably has, led to the destruction of many close-in exoplanets. Evidence for this destruction comes from the orbital distribution of observed close-in exoplanets: older planets tend to be farther from their host stars than younger planets; transiting planets tend to be younger than nontransiting planets; the observed orbits show a cutoff in small a-values, which increases to larger a-values with increasing age. These trends can be explained by tidal evolution and destruction of close-in exoplanets. Tidal decay rates are quantitatively consistent with the observed distribution of planetary semimajor axes and ages for a wide range of tidal dissipation parameters.

While it is widely assumed that some mechanism must stop the inward migration of close-in exoplanets in order to explain the observed orbital distribution, our results show that it is not necessary. Our model also explains the short life expectancies of some of the closest-in planets (e.g., WASP-12 b, CoroT-7 b, and OGLE-TR-56 b). They are simply the next in line to be tidally destroyed after the many that have already spiraled down into their stars. Transiting planets especially, because they usually pass close to their host stars to be observed in transit, will quickly be destroyed by tides, which may explain why transiting planets tend to be relatively young.

Future observations may provide further evidence for tidal destruction of close-in exoplanets. Stars that have recently consumed a close-in planet may have unusual compositions and may rotate more quickly than expected. Tidal stripping of a gaseous planet's mass once the planet nears the Roche limit may strand the rocky cores of gaseous planets with orbits very close to their host stars, creating a distinctive class of planets that may be seen with improving surveys.

As stellar age estimates become more accurate, more close-in exoplanets are discovered, and our understanding of tidal effects improve, these results will need to be revisited. Nevertheless, our results suggest tidal destruction of close-in exoplanets is common, and many of the exoplanets we see today may soon be gone.

The authors acknowledge useful conversations with Josh Eisner, Chris Laws, Michael Meyer, Dylan Morgan, Andrew Skemer, and Christa Van Laerhoven. The paper also benefited from comments by an anonymous referee. This work was supported by a grant from NASA's Planetary Geology and Geophysics program and a fellowship to B.J. from the NASA Earth and Space Science Fellowship Program. R.B. acknowledges funding from NASA Astrobiology Institute's Virtual Planetary Laboratory lead team, supported by NASA under Cooperative Agreement No. NNH05ZDA001C.