Introduction
Planetary habitability is a complex function of orbits, composition, atmospheric evolution and geophysical processes. Most searches for habitable environments begin with the comparison of a planet's orbit relative to the host star's habitable zone (HZ), the region around a star for which an Earth-like planet can support water on its surface (Dole Reference Dole1964; Kasting et al. Reference Kasting, Whitmire and Reynolds1993; Selsis et al. Reference Selsis2007; Kopparapu et al. Reference Kopparapu2013). A critical feature of this definition is the presence of a solid surface, i.e., the planet must be rocky. Hence, this determination is also crucial in the identification of potentially habitable environments.
Unfortunately a robust and universal definition of the boundary between rocky and gaseous worlds has remained elusive. Initially, research was theoretical and identified the mass of a solid body, a ‘core,’ that was large enough to permit the capture of protoplanetary gas (e.g., Pollack et al. Reference Pollack1996; Ikoma et al. Reference Ikoma, Emori and Nakazawa2001; Guillot Reference Guillot2005; Hubickyj et al. Reference Hubickyj, Bodenheimer and Lissauer2005; Militzer et al. Reference Militzer, Hubbard, Vorberger, Tamblyn and Bonev2008; Lissauer et al. Reference Lissauer, Hubickyj, D'Angelo and Bodenheimer2009; Movshovitz et al. Reference Movshovitz, Bodenheimer, Podolak and Lissauer2010). These results showed that a wide range of ‘critical masses’ is possible, from ≲1–20 MEarth. Undoubtedly, the actual critical mass is a function of the local protoplanetary disc's properties (e.g., temperature and viscosity), and one should expect the critical mass to vary from system to system. However, there should exist a maximum mass (or radius) below which planets are terrestrial-like because small mass planets do not possess enough gravitational force to hold on to hydrogen and helium. In this study, I examine the possibility that this boundary can be identified with transit data coupled with expectations from tidal theory.
Transiting exoplanets offer the opportunity to measure both the planetary radius and mass. The planetary radius can be constrained if the stellar radius is known, usually determined from spectral information (Torres et al. Reference Torres, Andersen and Giménez2010; Everett et al. Reference Everett, Howell, Silva and Szkody2013), but sometimes directly through interferometric observations (von Braun et al. Reference von Braun2011; Boyajian et al. Reference Boyajian2012). The mass can then be measured through radial velocity measurements, which no longer suffer from the mass-inclination degeneracy as the viewing geometry is known (e.g., Batalha et al. Reference Batalha2011). In multiple planet systems, masses can also be measured by transit timing variations (Agol et al. Reference Agol, Steffen, Sari and Clarkson2005; Holman & Murray Reference Holman and Murray2005), as for Kepler-9 (Holman et al. Reference Holman2010) and Kepler-11 (Lissauer et al. Reference Lissauer2011). However, for many transiting planets, these methods may not be available, as terrestrial planets tend to produce small radial velocity and timing variation signals. Thus, direct measurements of the masses of small planets orbiting FGK stars (0.7–1.4 MSun) will be challenging.
The transit detection method is biased towards the discovery of planets on close-in orbits, ≲0.1 AU for FGK stars. NASA's Kepler spacecraft has detected over 1000 planet candidates in this range, opening up the possibility that statistical analyses of these (uninhabitable) planets could reveal the critical radius between terrestrial and gaseous planets, R crit. As theoretical models of planet formation have found that the critical mass and radius lies near 10 MEarth and 2 REarth, planets that are smaller could be rocky. Although very few data points exist, it does appear that R crit<2 REarth, e.g Kepler-10 b with mass M p=4.56 MEarth and radius R p=1.42 REarth (Batalha et al. Reference Batalha2011) and CoRoT-7 b with M p<8 MEarth and R p=1.58 REarth (Léger et al. Reference Léger2009; Queloz et al. Reference Queloz2009; Ferraz-Mello et al. Reference Ferraz-Mello, Tadeu Dos Santos, Beaugé, Michtchenko and Rodríguez2011). Here I will call rocky planets larger than the Earth ‘super-Earths’, and gaseous planets less than 10 MEarth ‘mini-Neptunes.’
Planets amenable to both transit and radial velocity measurements tend to lie close to their stars, which is a different environment than any of the planets in our Solar System. These planets are subjected to more radiation, stellar outbursts (Ribas et al. Reference Ribas, Guinan, Güdel and Audard2005), and tidal effects (Rasio et al. Reference Rasio, Tout, Lubow and Livio1996; Jackson et al. Reference Jackson, Greenberg and Barnes2008b). While these effects can act in tandem (Jackson et al. Reference Jackson2010), here I only consider the tidal effects, as in isolation they can be used to calculate R crit.
The key feature is that the expected rates of tidal dissipation in terrestrial planets are orders of magnitude larger than for gaseous worlds. In the classical equilibrium tidal theory (e.g. Darwin Reference Darwin1880; Goldreich & Soter Reference Goldreich and Soter1966; Hut Reference Hut1981; Jackson et al. Reference Jackson, Greenberg and Barnes2008b), tidal dissipation is inversely proportional to the ‘tidal quality factor’ Q. For rocky planets in our Solar System, Q r~100 (Yoder Reference Yoder and Ahrens1995; Henning et al. Reference Henning, O'Connell and Sasselov2009), whereas giants have Q g=104−107 (Goldreich & Soter Reference Goldreich and Soter1966; Yoder Reference Yoder and Ahrens1995; Zhang & Hamilton Reference Zhang and Hamilton2008; Lainey et al. Reference Lainey2012), with a traditional value of 106 (Rasio et al. Reference Rasio, Tout, Lubow and Livio1996; Jackson et al. Reference Jackson, Greenberg and Barnes2008b, Reference Jackson, Barnes and Greenberg2009). For a typical star–planet configuration, tides usually damp both orbital eccentricity and semi-major axis, hence the different dissipation rates result in different damping timescales. In other words, terrestrial planets should circularize more quickly and/or at larger separations than gaseous planets. This discrepancy could reveal the value of R crit, and also Q r and Q g, all of which are still poorly constrained observationally. I find that the different orbital evolutions of rocky and gaseous exoplanets should be detectable, but is not apparent in publicly-available Kepler data.
In the next section, I describe the tidal models used in this study and the expected orbits of exoplanets in the Kepler field of view. In ‘The transit duration anomaly’ section, I describe how transit data can constrain eccentricity through the transit duration and basic orbital mechanics, including an extended discussion of observational biases. In ‘Methodology’ section, I describe how I produce hypothetical distributions of transiting planets that undergo tidal evolution. In ‘Results’ section, I present the results and find that the different orbital evolutions of rocky and gaseous exoplanets should be detectable, but are not apparent in publicly-available Kepler data. Those data contain a distribution of impact parameters that are inconsistent with isotropic orbits, suggesting a systematic error may be present. In the ‘Discussion’ section, I discuss the results and finally in the ‘Conclusion’ section, I draw my conclusions.
Tidal theory
For my calculations of tidal evolution, I employ ‘equilibrium tide’ models, originally conceived by Darwin (Reference Darwin1880). This model assumes that the gravitational potential of a perturber can be expressed as the sum of Legendre polynomials (i.e., surface waves) and that the elongated equilibrium shape of the perturbed body is slightly misaligned with respect to the line that connects the two centres of mass; see Fig. 1. This misalignment, which only occurs for non-zero eccentricity and obliquity, is due to dissipative processes within the deformed body and leads to a secular evolution of the orbit as well as the spin angular momenta of the two bodies. These assumptions produce six coupled, non-linear differential equations, but note that the model is, in fact, linear in the sense that there is no coupling between the surface waves which sum to the equilibrium shape. Considerable research has explored the validity and subtleties of the equilibrium tide model (e.g., Hut Reference Hut1981; Ferraz-Mello et al. Reference Ferraz-Mello, Rodríguez and Hussmann2008; Wisdom Reference Wisdom2008; Efroimsky & Williams Reference Efroimsky and Williams2009; Leconte et al. Reference Leconte, Chabrier, Baraffe and Levrard2010). For this investigation, I will use the models and nomenclature of Heller et al. (Reference Heller, Leconte and Barnes2011) and Barnes et al. (Reference Barnes2013), which are summarized below.
Fig. 1. Schematic representation of the position of a planetary tidal bulge due to a star. Dissipative processes inside the planet prevent the bulge, marked by the dotted line, from pointing towards the centre of the perturber, shown by the dashed line. Instead the bulge points away by an angle that is either inversely proportional to the tidal Q (for the CPL model), or proportional to the time lag τ (for the CTL model).
The constant–phase–lag model
In the ‘constant-phase-lag’ (CPL) model of tidal evolution, the angle between the line connecting the centres of mass and the tidal bulge is constant and the planet responds to the perturber like a damped, driven harmonic oscillator. The CPL model is commonly used in planetary studies (e.g., Goldreich & Soter Reference Goldreich and Soter1966; Greenberg Reference Greenberg2009). Under this assumption, and ignoring the effect of obliquity, the evolution is described by the following equations:
$$ \eqalign{ \displaystyle{{{\rm d}e} \over {{\rm d}t}}\;& = \; - \displaystyle{{ae} \over {8GM_1 M_2}} \mathop \sum \limits_{i = 1}^2 Z_i^' \left( {2{\rm \varepsilon} _{0,i} - \displaystyle{{49} \over 2}{\rm \varepsilon} _{1,i} + \displaystyle{1 \over 2}{\rm \varepsilon} _{2,i} + 3{\rm \varepsilon} _{5,i}} \right) \cr \displaystyle {{{\rm d}a} \over {{\rm d}t}}\;& = \; \displaystyle{{a^2} \over {4GM_1 M_2}} \mathop \sum \limits_{i = 1}^2 Z_i^' \;\cr &\times \left( {4{\rm \varepsilon} _{0,i} + e^2 \left[ { - 20{\rm \varepsilon} _{0,i} + \displaystyle{{147} \over 2}{\rm \varepsilon} _{1,i} + \displaystyle{1 \over 2}{\rm \varepsilon} _{2,i} - 3{\rm \varepsilon} _{5,i}} \right]} \right)}\hskip-6 $$where e is the eccentricity, t is the time, a is the semi-major axis, G is the Newton's gravitational constant, M 1 and M 2 are the two masses, and R 1 and R 2 are the two radii. The quantity Z ′i is
$$Z_i^' \equiv 3G^2 k_{2,i} M_j^2 (M_i + M_j )\displaystyle{{R_i^5} \over {a^9}} \; \displaystyle{1 \over {nQ_i}} $$where k 2,i are the Love numbers of order 2, n is the mean motion, and Q i are the tidal quality factors. The signs of the phase lags are
$$\matrix{ {{\rm \varepsilon} _{0,i} = & {\rm \Sigma} (2{\rm \omega} _i - 2n)} \cr {{\rm \varepsilon} _{1,i} = & {\rm \Sigma} (2{\rm \omega} _i - 3n)} \cr {{\rm \varepsilon} _{2,i} = & {\rm \Sigma} (2{\rm \omega} _i - n)} \cr {{\rm \varepsilon} _{5,i} = & {\rm \Sigma} (n)} \cr {{\rm \varepsilon} _{8,i} = & {\rm \Sigma} ({\rm \omega} _i - 2n)} \cr {{\rm \varepsilon} _{9,i} = & {\rm \Sigma} ({\rm \omega} _i )\;} \cr} $$where ωi is the rotational frequency of the ith body, which I force to be the equilibrium frequency, (1+9.5e 2)n. Σ(x) is the sign of any physical quantity x, and thus Σ(x)=+1, −1, or 0.
The constant–time–lag model
The ‘constant–time–lag’ (CTL) model assumes that the time interval between the passage of the perturber and the tidal bulge is a constant value, τ. This assumption allows the tidal response to be continuous over a wide range of frequencies, unlike the CPL model. But, if the phase lag is a function of the forcing frequency, then the system is no longer analogous to a damped driven harmonic oscillator. Therefore, this model should only be used over a narrow range of frequencies; see Greenberg (Reference Greenberg2009). Ignoring obliquity, the orbital evolution is described by the following equations:
$$\displaystyle{{{\rm d}e} \over {{\rm d}t}}\; = \; \displaystyle{{11ae} \over {2GM_1 M_2}} \mathop \sum \limits_{i = 1}^2 Z_i \left( {\displaystyle{{\,f_4 (e)} \over {{\rm \beta} ^{10} (e)}}\displaystyle{{{\rm \omega} _i} \over n} - \displaystyle{{18} \over {11}}\displaystyle{{\,f_3 (e)} \over {{\rm \beta} ^{13} (e)}}} \right)$$
$$\displaystyle{{{\rm d}a} \over {{\rm d}t}}\; = \; \displaystyle{{2a^2} \over {GM_1 M_2}} \mathop \sum \limits_{i = 1}^2 Z_i \left( {\displaystyle{{f_2 (e)} \over {{\rm \beta} ^{12} (e)}}\displaystyle{{{\rm \omega} _i} \over n} - \displaystyle{{f_1 (e)} \over {{\rm \beta} ^{15} (e)}}} \right)$$where
$$Z_i \equiv 3G^2 k_{2,i} M_j^2 (M_i + M_j )\displaystyle{{R_i^5} \over {a^9}} \; {\rm \tau} _i \; $$and
$$\matrix { {\beta (e) = \sqrt {1 - e^2}} \cr {\,f_1 (e) = 1 + \displaystyle{{31} \over 2}e^2 + \displaystyle{{255} \over 8}e^4 + \displaystyle{{185} \over {16}}e^6 + \displaystyle{{25} \over {64}}e^8} \cr {\,f_2 (e) = 1 + \displaystyle{{15} \over 2}e^2 + \displaystyle{{45} \over 8}e^4 \; \; + \displaystyle{5 \over {16}}e^6} \cr {\,f_3 (e) = 1 + \displaystyle{{15} \over 4}e^2 + \displaystyle{{15} \over 8}e^4 \; \; + \displaystyle{5 \over {64}}e^6} \cr {\,f_4 (e) = 1 + \displaystyle{3 \over 2}e^2 \; \; + \displaystyle{1 \over 8}e^4} \cr {\,f_5 (e) = 1 + 3e^2 \; \; \; + \displaystyle{3 \over 8}e^4} \cr} $$As in the CPL model, I force the rotational frequency to equal the equilibrium frequency, which is (1+6e 2)n in the CTL model.
There is no general conversion between Q p and τp. Only if e=0 (and the obliquity is 0 or π), when merely a single tidal lag angle exists, then
$$Q_{\rm p} {\rm \;} \approx {\rm \;} 1/(2|n - {\rm \omega} _{\rm p} |{\rm \tau} _{\rm p} )$$as long as n−ωp remains unchanged. Hence, the canonical values of the dissipation parameters for dry, rocky planets in the Solar System, Q=100 (Goldreich & Soter Reference Goldreich and Soter1966) and τ=638 s (Lambeck Reference Lambeck1977, Barnes et al. Reference Barnes2013), are not necessarily equivalent. Hence, the results for the tidal evolution will intrinsically differ between the CPL and the CTL model, even though both choices are common for the respective model. While tidal dissipation in super-Earths remains observationally unconstrained, here I will assume that it is similar to the rocky bodies in the Solar System.
Differential circularization of super-Earths and mini-Neptunes
As described above, most previous research predicts several orders of magnitude difference in Q or τ between terrestrial and giant planets. As an example in Fig. 2, I simulate the evolution of two hypothetical planets with different compositions that formed with identical orbits around identical stars. The line shows 10 Gyr of CPL tidal evolution of a 2 R Earth planet with a density of 1 g cm−3 (i.e., a gaseous 3.8 M Earth mini-Neptune) and a tidal Q of 106 (see, e.g., Goldreich & Soter Reference Goldreich and Soter1966; Jackson et al. Reference Jackson, Greenberg and Barnes2008b), while the filled circles are the orbit of a 2 REarth planet with a mass of 10 MEarth and a tidal Q of 100 (i.e., a super-Earth) every 100 Myr. The super-Earth circularizes in about 1 Gyr; the mini-Neptune barely evolves, even over 10 Gyr. This discrepancy is despite that the equilibrium tidal models predict evolution scales with planetary mass – the large difference between the Qs dominates.
Fig. 2. Comparison of the tidal evolution of a gaseous 2 REarth mini-Neptune (solid line; for 10 Gyr) and a rocky 2 REarth super-Earth (open squares; in 100 Myr intervals). Both planets begin with an orbital period of 5 days and eccentricity of 0.2. The mini-Neptune experiences little orbital evolution, but the super-Earth circularizes in about 1 Gyr. This discrepancy is due to the four orders of magnitude difference in tidal dissipation between gaseous and rocky planets.
The transit duration anomaly
In this section, I review and revise previous work on the ‘transit duration anomaly’ (TDA), defined here as the ratio of the observed transit duration to the duration if the orbit were circular. This parameter has gone by several names in the literature, such as the ‘transit duration deviation’ (Kane et al. Reference Kane, Ciardi, Gelino and von Braun2012), and the ‘photoeccentric effect’ (Dawson & Johnson Reference Dawson and Johnson2012). Here I use the name proposed by Plavchan et al. (Reference Plavchan2012)Footnote 1, as in celestial mechanics the term ‘anomaly’ refers to a parameter's value relative to pericentre, e.g., the ‘true anomaly’ is the difference between a planet's true longitude and its longitude of pericentre. The analogy is not perfect, as the TDA is not measured relative to the duration at pericentre, but rather to the duration due to acircular orbit. Nonetheless, ‘anomaly’ captures the fact that the duration is set by the longitude relative to pericentre, the true anomaly θ, as shown below.
In this section, I first review how the TDA can be used to determine the minimum eccentricity of an orbit. Then I review the biases implicit in TDA measurements and update previous results.
Determination of the minimum eccentricity
Transit data coupled with knowledge of semi-major axis enable the imposition of a lower bound on the eccentricity e min (Jason W. Barnes 2007; Ford et al. Reference Ford, Quinn and Veras2008). The determination of e min requires knowledge of both the physical size of the orbit, as well as a precise determination of the orbital period P, planetary and stellar radii (R p and R *), and the impact parameter b; see Fig. 3. If the transit is well-sampled, then these parameters can be obtained (Mandel & Agol Reference Mandel and Agol2002).
Fig. 3. Schematic representation of a planetary transit. The planet has a radius R p, and the star R *. The planet crosses the stellar disc a projected distance b from disc centre. The distance the planet travels during transit is D, the distance between the centres of the planet at first and fourth contacts.
The transit duration is the time required for a planet to traverse the disc of its parent star, and to first order is:
$$T = \displaystyle{D \over {v_{{\rm sky}}}} = \displaystyle{{2\left( {\sqrt {(R_* + R_p )^2 - b^2}} \right) \over {v_{{\rm sky}}}} $$where v sky is the azimuthal velocity, i.e., the instantaneous velocity of the planet in the plane of the sky. Although several different definitions of the duration are possible (Kipping Reference Kipping2010), I choose this definition to match the Kepler public data. On a circular orbit, the azimuthal velocity is constant and equal to the orbital velocity. Therefore, the duration for a circular orbit is
$$T_c = \displaystyle{{\sqrt {(R_* + R_p )^2 - b^2}} \over {{\rm \pi} a}}P$$where P is the orbital period. For an eccentric orbit the orbital velocity is a function of longitude (Kepler's Second Law), and is given by
$$v = v_c \sqrt {\displaystyle{{1 + 2e\,\cos \; {\rm \theta} + e^2} \over {1 - e^2}}} $$where e is the eccentricity and v c is the circular velocity. Finally, from classical mechanics, the azimuthal velocity is
$$v_{{\rm sky}} = v_c \displaystyle{{1 + e\,{\rm cos}\; {\rm \theta}} \over {\sqrt {1 - e^2}}} .$$From transit data alone, the value of θ is unknown, and hence so is e.
However, one can exploit the difference between T and T c to obtain a minimum value of the eccentricity, e min (J. W. Barnes 2007). The situation is somewhat complicated because T can be larger or smaller than T c depending on θ. If the planet is close to apoapse, T>T c, while near periapse T<T c. To derive e min, one must assume that θ=0 or π. While the velocity could be larger at some other position in the orbit, the maximum deviation from the circular velocity is at least as large as the measured velocity, and hence e must be at least a certain value. If I define the TDA as
$${\rm \Delta} \equiv T/T_c = v_c /v_{{\rm sky}} = \displaystyle{{\sqrt {1 - e^2}} \over {1 + { e}\,{\rm cos}\; {\rm \theta}}} $$then
$$e_{{\rm min}} = \left| {\displaystyle{{{\rm \Delta} ^2 - 1} \over {{\rm \Delta} ^2 + 1}}} \right|$$is the minimum eccentricity permitted by transit data.
Several studies invoked the orbital velocity instead of the sky velocity (e.g., Tingley & Sackett Reference Tingley and Sackett2005; Burke Reference Burke2008; Kipping Reference Kipping2010). However, J. W. Barnes (2007) correctly noted that T is actually a function of the azimuthal velocity, which equals the orbital velocity at pericentre and apocentre. For many other studies, the assumed form of the velocity is unclear. As transits are a photometric phenomenon, unaffected by the component of the velocity along the line of sight, the projection of the orbital velocity into the sky plane, v sky, is the appropriate choice. This implies that previous studies are only approximately correct. As I show in the next section, using v c will only amount to a small error.
The TDA has been used in several studies to constrain the eccentricity distribution, often with the assumption that the impact parameter is unknown, as proposed by Ford et al. (Reference Ford, Quinn and Veras2008). In that case, one can only use those systems in which T>T c for a central transit (b=0) to estimate e min. Moorhead et al. (Reference Moorhead2011) analysed the first three quarters of Kepler data and found that the KOIs appeared to be consistent with a mean eccentricity near 0.2. They also found that eccentricities appear to be large regardless of orbital period, and that small planets tend to have larger eccentricities. More recent work has failed to determine if the Kepler eccentricity distribution is consistent with the radial velocity planets (Kane et al. Reference Kane, Ciardi, Gelino and von Braun2012; Plavchan et al. Reference Plavchan2012). These studies were limited by the number of known candidates, as well as the relatively poor characterization of the transits themselves. Note that Kane et al. (Reference Kane, Ciardi, Gelino and von Braun2012) used the difference between T and T c, rather than the quotient, to model eccentricity. Dawson & Johnson (Reference Dawson and Johnson2012) demonstrated that in some cases careful statistical analyses of transit data can provide constraints on the actual eccentricity, especially if Δ deviates significantly from unity.
Observational biases in the transit duration anomaly
To further elucidate the TDA, as well as to revisit previous results that may have invoked an inappropriate definition, in this section I review the geometry and biases associated with it. In Fig. 4, I show the orbital and azimuthal velocities as a function of true anomaly for three different eccentricities. The differences between the panels are subtle at low e, but can be significant for large e.
Fig. 4. Velocity relative to the circular velocity as a function of true anomaly for three different eccentricities. The dotted line shows the circular velocity. Dashed lines show the values of true anomaly at which the velocity equals the circular velocity. Left: Orbital velocity. Right: Velocity in the plane of the sky.
The vertical lines in Fig. 4 show the values of θ at which the velocity is equal to the circular velocity. For the orbital velocity, the longitude where v=v c is given by
$${\rm \theta} _{\bi c}^{{\rm orb}} = \pm \cos ^{ - 1} ( - {e})$$and for the azimuthal velocity, it lies at
$${\rm \theta} _{\bi c}^{{\rm sky}} = \pm \cos ^{ - 1} \left( {\displaystyle{{\sqrt {1 - e^2} - 1} \over { e}}} \right).$$In Fig. 5, I show schematics for four orbits. In all cases an observer at x=+∞ views the transit such that either v=v c (left) or v sky=v c (right).
Fig. 5. Locations on an orbit where a velocity is equal to the circular velocity. In all cases, the orbit's semi-major axis is 1, and the observer is assumed to be located at x=+∞, as shown by the dashed line. The longitude of periastron is shown by the dotted line. In each panel, the eccentricity is listed, as is one of the values of the true anomaly at which the velocity is circular. In all cases, the orbit is rotated so that that longitude, θc, lies on the +x-axis. Left: The total velocity of the planet is equal to the circular velocity when it crosses the +x-axis. In these cases, the observer will see transit durations longer than that from an identical planet on a circular orbit. Right: The azimuthal velocity is equal to the circular velocity when it crosses the +x-axis. In these cases, an observer will see a transit duration equal to that of an identical planet, but on a circular orbit, i.e., e min=0.
As I am only interested in the azimuthal velocity, for the remainder of this paper I will assume that θc=θcsky and drop the superscript. For an eccentric orbit, the planet travels faster than the circular velocity for θc/π>0.5 of the orbit. There is therefore an observational bias, which I will call the ‘velocity bias,’ to observe Δ<1 (Tingley & Sackett Reference Tingley and Sackett2005; Burke Reference Burke2008; Plavchan et al. Reference Plavchan2012). The probability that T<T c due to this effect is just
$$p_{{\rm vel}} (T \lt T_c ) = \displaystyle{{{\rm \theta} _c} \over {\rm \pi}} $$and is shown by the dashed curve in Fig. 6.
Fig. 6. Probability of detecting a transit duration shorter than that of a planet on a circular orbit. The dashed curve is the ‘velocity bias’ that arises because more longitudes have v sky>v c. The solid curve is the ‘duration bias’ that also takes into account the increased observability of transits when the star–planet separation is small. The solid curve is the actual bias, which is not a straight line.
The bias towards small Δ is magnified by the geometrical bias towards transits occurring at smaller star–planet separations. As shown in J. W. Barnes (2007; see also Borucki & Summers Reference Borucki and Summers1984), the overall transit probability is
$$\eqalign{\,p_{{\rm transit}} = & \displaystyle{1 \over {4{\rm \pi}}} \displaystyle{{2R_*} \over {a(1 - e^2 )}}\int_0^{2{\rm \pi}} {(1 + e\cos {\rm \; \theta} )} \,d{\rm \theta} \cr = & \displaystyle{{R_*} \over {a(1 - e^2 )}}} $$and the probability to observe the transit when T<T c is
$$\eqalign{\,p_{{\rm short}} = & \displaystyle{1 \over {4{\rm \pi}}} \displaystyle{{2R_*} \over {a(1 - e^2 )}}\int_{ - {\rm \theta} _c} ^{{\rm \theta} _c} {(1 + e\cos {\rm \; \theta} )\,d{\rm \theta}} \cr = & \displaystyle{{R_* ({\rm \theta} _c + e\,{\rm sin}\; {\rm \theta} _c )} \over {a{\rm \pi} (1 - e^2 )}}} $$and thus the bias towards observing a transit duration shorter than the circular duration is the ratio of these two equations,
$$p_{{\rm duration}} = \displaystyle{{{\rm \theta} _c + e\,{\rm sin}\; ({\rm \theta} _c )} \over {\rm \pi}} $$which I will call the ‘duration bias.’ This effect is shown by the solid curve in Fig. 6, and is the actual likelihood to observe a transit with T<T c. As e→1, it becomes extremely unlikely to observe a long transit. This effect can make studies that rely on transit durations longer than that predicted for a central transit (b=0) unlikely to find high eccentricity objects (e.g., Moorhead et al. Reference Moorhead2011; Plavchan et al. Reference Plavchan2012).
In practice, this bias is not dramatic as most planets are not on very eccentric orbits. Burke (Reference Burke2008) used analytic fits to the then-current eccentricity distribution as determined from radial velocity planets, excluding those with orbital periods less than 10 days that may be tidally circularized, to find that the mean value of Δ should be 0.88. Recall that Burke (Reference Burke2008) used v instead of v sky in his definition of Δ, thus his mean should be slightly lower than the actual mean as v≥v sky.
To update the expectations of Burke (Reference Burke2008), I recompute the expected distribution of Δ from radial velocity planets. In the left panel of Fig. 7, I show the current distribution of e with the thick grey lineFootnote 2. I exclude those planets with a>0.1 AU leaving 362 planets. To evaluate the expected distribution of Δ, I created 107 synthetic systems consisting of a star with a radius between 0.7 and 1.4 Rʘ, and a planet whose radius I ignored. The orbit had a=0.05 AU, an eccentricity distribution given by the dashed histogram in the left panel of Fig. 7, and an isotropic distribution of orbits. I then calculated Δ for all transiting geometries with durations larger than 1 hour and show the resulting Δ distribution in the right panel of Fig. 7. As noted, 66% of cases have Δ<1, with a mean value of 0.90. If I instead use the circular velocity to calculate Δ, I find 68% have Δ<1, with a mean of 0.88, reproducing the Burke (Reference Burke2008) result. The difference between using v c and v sky to calculate the TDA is modest for the known radial-velocity-detected planets.
Fig. 7. Left: Eccentricity distribution of exoplanets. The solid grey line is the observed distribution (with bin size 0.05) of 362 radial velocity exoplanets with a>0.1 AU from exoplanets.org as of 4 June 2013. The dashed line is the distribution used in this study. Right: The TDA expected from the simulated data set assuming a random viewing geometry and that the transit duration is set by the azimuthal velocity, equation (12).
Methodology
In order to determine if the difference in Q values can permit the identification of R crit, I perform Monte Carlo simulations of both the CPL and CTL models. I then compare the results to publicly-available Kepler data to search for the predicted signal.
For my synthetic data, I created 25 000 star–planet configurations with initial semi-major axes uniformly in the range [0.01, 0.15] AU, planetary radii in the range [0.5, 10] REarth, stellar masses in the range [0.7, 1.4] MSun, a radius in solar radii equal to its mass in solar masses, and ages in the range [2, 8] Gyr. If the planetary radius is less than 2 REarth, then the mass is (R/REarth)3.68 MEarth (Sotin et al. Reference Sotin, Grasset and Mocquet2007), if larger, then I assume that the density is 1 g cm−3, similar to the planets in the Kepler-11 system (Lissauer et al. Reference Lissauer2013). The initial eccentricity is drawn from the currently observed distribution of distant planets (a>0.1 AU); see Fig. 7. For CPL runs I used 30≤Q r≤300, 106≤Q g≤107, and 106≤Q *≤107. For CTL runs I used 30≤τr≤300 s, 0.003≤τg≤0.03 s and 0.001≤τ*≤0.01 s. I then integrated the system forward for the randomly chosen age and assumed that I observed the system in that final configuration. In order to calculate e min, I choose a random value for θ that represents the direction of the observer, and an inclination i chosen uniformly in cos i, with i measured from the plane of the sky. I then calculate the separation between the star and planet using
$$r = \displaystyle{{a(1 - e^2 )} \over {1 + {e}\,{\rm cos}\; {\rm \theta}}} $$and determine the impact parameter, b=r tan (π/2−i). If b<R *, then the planet transits, and I calculate the transit duration. As pointed out in Burke (Reference Burke2008), short transit durations can be missed, and I therefore throw out transit durations that are less than 1 hour, which is the approximate minimum duration detectable by Kepler. The vast majority of the rejected transits are too short due to a large impact parameter; however, a few are due to large eccentricity and alignment of the longitude of pericentre with the line of sight. Thus, my estimates of the minimum eccentricity distribution are slightly biased towards lower values. From the remaining transits, I calculate the TDA and e min using equations (13–15).
I also compute e min for Kepler candidates that have all the requisite parameters presented in Kepler Planet Candidate Data ExplorerFootnote 3 (see also Batalha et al. Reference Batalha2013). These data do not contain error bars and assuredly contain some false positive, but the data set is uniform and sufficient for this proof of concept. I limit my sample to those with orbital periods less than 15 days, but will refer to this subsample as the ‘Kepler sample’ in the upcoming sections.
Results
Tides and the minimum eccentricity
I begin by considering an intermediate step: In Fig. 8, I show the average final eccentricity of my simulated planets as a function of planetary radius, R p, and orbital period, P. The paucity of eccentric orbits at low R p and P is due to the more effective circularization of rocky bodies. Furthermore, I can see the features that correspond directly to three parameters that are currently very poorly constrained: R crit via the rapid rise in <e> at the imposed value of 2 REarth; Q g (τg) via the rapid rise in <e> at 1 day above 2 REarth; and Q r (τr) via the rise over 2–10 days and below 2 REarth. Thus, despite the order of magnitude uncertainty I gave to each physical parameter, the large discrepancy between Q r (τr) and Q g (τg) does result in an important difference in the expected orbits of close-in planets of FGK stars. For the CTL model, 1955 planets merged with the star during the integrations; 1118 for CPL (see Jackson et al. Reference Jackson, Barnes and Greenberg2009; Levrard et al. Reference Levrard, Winisdoerffer and Chabrier2009).
Fig. 8. Average final eccentricity of a suite of ~23 000 systems of one star and one planet that experienced tidal evolution for 2–8 Gyr. Planets larger than 2 REarth are gaseous, and those smaller are rocky. The binsizes are 0.5 REarth in radius and 0.5 day in period. Gaseous planets retain a residual eccentricity if P>1.5 days, while rocky planets require P>3 days. The transition occurs at 2 REarth, which in this case is R crit. Left: Orbits evolve according to the CTL model. Right: Orbits evolve according to the CPL model.
Next I calculate the average minimum eccentricity <e min> for transiting geometries of simulated rocky and gaseous planets in 0.5 day orbital period intervals and plot <e min> as a function of orbital period for different radii as solid lines in Fig. 9. For the CTL model, I obtained 2 127 observable transits, and for CPL 2 151, about twice as many planets as in the same period range as the Kepler sample. Note that my synthetic data do not share several properties of the Kepler data, such as the planetary radius and period distributions. For both the CTL (left) and CPL (right) models, the trends are the same. For R<R crit, <e min>~0 up to about a 4–5-day period. However, for larger radii, circular orbits are only guaranteed for periods less than about 1.5–2 days. The rocky and gaseous distributions become about equal at P=13 days, albeit with considerable scatter. At large orbital periods, the simulated data become sparser as the transit probability is dropping (producing the apparent oscillations in <e min>), and those that do transit are more likely to be near pericentre of an eccentric orbit, causing the secular growth in <e min> with period. Note that similar variations are present in the Kepler sample.
Fig. 9. Average minimum eccentricities for transiting exoplanets as a function of period and radius. The solid curve and filled squares represent planets with radii below the selected critical radius of 2 REarth, while the dashed curve and open circles are larger planets. The lines are for the simulated data set incorporating the duration bias, symbols for KOIs. The latter have values of <e min> of 0.2–0.5 regardless of period, whereas model gaseous planets have non-zero eccentricities if the period is larger than 1.5 days, and rocky planets at larger than 4. If tidal dissipation is a function of exoplanet radius, it should be detectable. Left: Results for the CTL model. Right: Results for the CPL model.
Figure 9 also contains the values of e min provided by the Kepler team as squares. Solid squares correspond to R p<R crit, open to R p>R crit. Nearly all the observed data are above the predictions, in agreement with the results of Moorhead et al. (Reference Moorhead2011). Therefore, it does not appear that there is any signal of tidal evolution in the Kepler data, regardless of orbital period! This result is in stark contrast to radial velocity data that show clear signs of circularization at small P (e.g., Butler et al. Reference Butler2006). Moreover, the average eccentricity in the Butler et al. catalogue is ~0.25, which is lower than the vast majority of minimum eccentricities derived from Kepler transits in this study. As the radial velocity data are older and have been reproduced by multiple teams, the Kepler data are more likely to be incorrect. In the next section, I describe several plausible explanations for the discrepancy.
A closer look at the Kepler sample
The lack of evidence of tidal evolution in the KOIs suggests there is an issue with the interpretation of the light curves. In this section, I examine several features of the Kepler sample and conclude that the data suffer from a systematic bias. As described in Batalha et al. (Reference Batalha2013), transits are fit to the geometric, limb-darkened transit model of Mandel & Agol (Reference Mandel and Agol2002). This model can determine planetary radius and impact parameter, as well as other parameters that are not relevant to the current study. The publicly-available data are long cadence, and hence the transits are not well-sampled. This sparse sampling is most likely to affect the impact parameter, as the shape of the transit is crucial to its estimation. Below I show that the impact parameters do indeed appear to be suspicious.
A partial list of the Kepler data used in this study is shown in Table 1, with the full table available in the supplementary material. In Fig. 10, I plot Δ as a function of orbital period in the left panel. Although no trends are present, it does appear that most values of Δ are greater than 1. In the right panel, I bin the Δ values to confirm this impression. This distribution cannot be explained by orbital mechanics and isotropic orbits, which predict that Δ distributions can only be biased towards Δ<1. Instead I find 78% of KOIs have Δ>1, with a mean value of 1.38.
Fig. 10. The transit duration anomalies of the KOIs in our sample (which are short-period). Left: TDA as a function of orbital period. Right: The histogram of TDA values shows that most values are larger than 1.
Table 1. Properties of short period KOIs
To search for the source of this discrepancy, I considered relationships among the parameters that permit the calculation of e min. The durations increase monotonically with the orbital period, albeit with significant scatter, as expected (Kane et al. Reference Kane, Ciardi, Gelino and von Braun2012). Several studies have pointed out systematic errors in the stellar characterization (Dressing & Charbonneau Reference Dressing and Charbonneau2013; Everett et al. Reference Everett, Howell, Silva and Szkody2013). In particular, Everett et al. (Reference Everett, Howell, Silva and Szkody2013) studied 220 Kepler host stars and found the vast majority have larger radii than reported by the Kepler team, and that one-quarter are 35% larger than suggested by the Kepler team. Since Δ∝R *−1, such a revision could significantly lower Δ and potentially resolve the discrepancy. Since they ‘only’ examined 220 host stars, some of which are not known to host a close-in planet, I do not include their results here, so that my analysis is kept to a uniform sample.
Perhaps the biggest inconsistency in the Kepler data lies in the impact parameter, see Fig. 11. The distribution predicted by the radial velocity exoplanets beyond the reach of tides is shown with the dashed histogram and is taken from the same sample that produced the right panel of Fig. 7. It is approximately flat, with the slight rise towards small values due to isotropically distributed orbits favouring edge-on geometries. Instead, the Kepler sample, shown by the solid line, rises sharply to large values of b. This distribution hints that a systematic error may be present in the Kepler analysis, which manifests itself in my analysis into large values of Δ and e min. I conclude that the currently available Kepler data produce unreliable values of b and hence e min.
Fig. 11. Distribution of impact parameters. The expected distribution of the radial velocity exoplanets with no tidal circularization is shown by the dashed line. The short-period KOIs tend to have large impact parameters.
Discussion
My simulated data show that R crit is identifiable in transit data due to the difference in the tidal Qs of gaseous and rocky bodies, at least for my idealized model of exoplanetary properties. My analysis is somewhat circular as I split the data in Fig. 9 on the selected value of R crit. In reality, its value is unknown and must be searched for. However, given the large and approximately equal values of <e min> in the Kepler data, I did not perform that search. More accurate data are needed, and may be available as KOI parameters are refined. Resolution of transit ingress and egress may be possible with short cadence data (which are unpublished but assuredly a small fraction of the total number of KOIs), or by folding the hundreds of transits together (e.g., Jackson et al. Reference Jackson2013), potentially rectifying the discordance between the Kepler and radial velocity data.
In order to accurately determine the minimum eccentricities, one needs both reliable information for both R * and b, but they are not yet available. Although subsets of more reliable data are available for the former (e.g., Everett et al. Reference Everett, Howell, Silva and Szkody2013), the transit fits are still plagued by inaccurate calculations of the impact parameter. Determination of the impact parameters in short cadence data or by folding would require a new and comprehensive analysis of those light-curves and is beyond the scope of this study. After those data have been properly analysed, the technique described in this study should be re-applied in order to determine R crit, Q g and Q r.
Aside from systematic errors in the analysis of the light curves, physical effects can also impact the value of <e min>. First, I note that additional companions can pump eccentricity through mutual gravitational interactions, even if tidal damping is ongoing (Mardling & Lin Reference Mardling and Lin2002; Bolmont et al. Reference Bolmont2013). Therefore one must be cautious when interpreting Fig. 9, as additional companions, both seen and unseen, can maintain non-zero eccentricities. However, Bolmont et al. (Reference Bolmont2013) showed that planet–planet interactions cannot maintain the eccentricity of the hot super-Earth 55 Cnc e above 0.1. That system is particularly relevant as there are many close-in planets orbiting a typical G dwarf. Therefore, I conclude that eccentricity pumping can be significant, but cannot explain the discrepancy between the observed and simulated systems shown in Fig. 9.
Another possibility is that stellar winds and activity can strip an atmosphere, reducing the mass and radius, and potentially changing the planet from a mini-Neptune to a super-Earth (Jackson et al. Reference Jackson2010; Valencia et al. Reference Valencia, Ikoma, Guillot and Nettelmann2010; Leitzinger et al. Reference Leitzinger2011; Poppenhaeger et al. Reference Poppenhaeger2012). Recently, Owen & Wu (Reference Owen and Wu2013) argued that the Kepler sample is consistent with hydrodynamic mass loss, and that some low-mass planets could have formed with substantially more mass. Mass loss should increase the time to circularize the orbit, assuming that the radius doesn't become very large, which is unlikely after about 100 Myr (Lopez et al. Reference Lopez, Fortney and Miller2012). Therefore, mass loss could stall circularization for mini-Neptunes, but not for super-Earths. Although a few radial velocity measurements exist for planets with radii less than ~1.5 REarth, they have densities consistent with silicate compositions (Batalha et al. Reference Batalha2011). Thus, mass loss seems unlikely to explain the differences seen for the smallest candidates in the Kepler field.
Radial inflation by irradiation (e.g., Lopez et al. Reference Lopez, Fortney and Miller2012) or tidal heating (Bodenheimer et al. Reference Bodenheimer, Lin and Mardling2001; Jackson et al. Reference Jackson, Greenberg and Barnes2008c; Ibgui & Burrows Reference Ibgui and Burrows2009) also work to decrease e since the evolution scales as R p5. Hence, bloated planets should be found on circular orbits, but no such trend is observed in the Kepler candidates.
In this study, I used two qualitatively different equilibrium tidal models and standard assumptions for dissipation. However, different tidal models have been proposed (e.g., Ogilvie & Lin Reference Ogilvie and Lin2004; Henning et al. Reference Henning, O'Connell and Sasselov2009; Socrates et al. Reference Socrates, Katz, Dong and Tremaine2012; Makarov & Efroimsky Reference Makarov and Efroimsky2013) and could be applied to this problem. However, the trend I predict here holds unless the dissipation rates in gas giants and rocky planets are within 1–2 orders of magnitude of each other, rather 4–6. Recently, Storch & Lai (Reference Storch and Lai2013) have proposed just such a model in which all tidal dissipation in gas giants occurs in a rocky core. Should all planets show the same trend in e min, then that would be evidence in support of their model. Hence, even if the expectations laid out above prove to be incorrect, some tidal models could be rejected by the methodology used in this study.
I have focused on transiting planets, but an analogous study could be applied to radial velocity data. Those data may be more amenable to such a study as orbital eccentricity is a direct observable. The problem lies in the low reflex velocities induced by the small planets as well as the ambiguity in mass due to the mass-inclination degeneracy. Nonetheless, with enough objects and an accounting for the expected isotropy of orbits, it may be possible to determine tidal dissipation as a function of mass in radial velocity data.
The Kepler spacecraft was designed to discover a potentially habitable planet orbiting a solar-like star. Such a planet (m p≲10 MEarth; r p≲2 REarth; P≈1 year) would have an undetectable radial velocity signature, preventing a direct mass measurement. Thus, confirmation of that planet's rocky nature is daunting. However, Kepler data may also hold the key to a convincing solution to the problem, as shown in this study. The distribution of TDAs in conjunction with tidal theory suggests the value of R crit may be calculated from the close-in planets. Although these planets are not habitable, they may provide crucial information to assess the habitability of Earth-like planets that transit Sun-like stars.
Conclusions
I have shown that the expected difference in tidal dissipation between gaseous and terrestrial exoplanets should lead to tidal circularization at different orbital distances. I have also shown how transit data, namely the TDA, can be used to determine the critical radius between gaseous and terrestrial planets. Moreover, an analysis of e min can also constrain the tidal dissipation in exoplanets, an understanding of which is sorely needed. Using standard values for tidal parameters and the critical radius, I find that a large ensemble of transit data should identify the critical radius between rocky and gaseous exoplanets. My analysis of available Kepler data reveals that the theoretical expectations are not met. However, this discrepancy cannot be used to refute the hypothesis because the values of <e min> are inconsistent with radial velocity detected exoplanets, particularly where tidal damping has been observed.
I have also reviewed the derivation of the TDA, as well as the known biases towards small values. Previous studies have advocated different choices for the velocity of the transiting planet as a function of true anomaly. The transit duration is actually determined by the azimuthal velocity, which is the velocity in the sky plane (equation (12)). While it is not clear in all previous studies which velocity was used, those that used the orbital velocity will obtain slightly smaller expected values of Δ. I find that the current distribution of exoplanet eccentricities predicts that 66% of transit durations should be less than T c, with a mean of 0.9. In contrast, 78% of short-period KOIs have durations greater than T c with a mean of 1.38. This distribution is inconsistent with celestial mechanics and the expectation isotropic orbits, regardless of tidal damping.
As the Kepler data are refined, or as new data, e.g. from the TESS mission arrive, this hypothesis should be revisited. The value of R crit is crucial for the interpretation of the habitability of transiting Earth-sized planets orbiting Sun-sized stars. Moreover, planets in the habitable zones of M dwarfs are susceptible to tidal effects (Dole Reference Dole1964; Kasting et al. Reference Kasting, Whitmire and Reynolds1993; Jackson et al. Reference Jackson, Barnes and Greenberg2008a; Heller et al. Reference Heller, Leconte and Barnes2011; Barnes et al. Reference Barnes2013), so a determination of Q for terrestrial exoplanets is also crucial to assessing habitability of planets such as those orbiting Gl 581 (Udry et al. Reference Udry2007; Mayor et al. Reference Mayor2009; Vogt et al. Reference Vogt2010) and Gl 667C (Anglada-Escudé et al. Reference Anglada-Escudé2012, Reference Anglada-Escudé2013; Bonfils et al. Reference Bonfils2013). As missions like Kepler and TESS have been designed to find potentially habitable worlds, the determination of R crit through their data alone would be an important step forward in determining the occurrence rate of terrestrial planets in the HZ.
Supplementary material
Supplementary materials of this paper are available at http://journals.cambridge.org/IJA
Acknowledgements
I thank Andrew Becker, Eric Agol, Leslie Hebb, Jason Barnes, Brian Jackson and René Heller for helpful discussions. This work was supported by NSF grant AST-110882 and the NAI's Virtual Planetary Laboratory lead team. I also thank an anonymous referee and Nader Haghighipour for reviews that greatly improved the clarity and accuracy of this manuscript.
