Submitted to the Astrophysical Journal Orbital Evolution and Migration of Giant Planets: Modeling Extrasolar Planets D. E. Trilling Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721 8 9 W. Benz 9 Steward Observatory, University of Arizona, Tucson, AZ 85721 and Physikalisches Institut, 1 Universitaet Bern, Sidlerstrasse 5, CH-3012 Bern, Switzerland n a J T. Guillot 8 University of Reading, Department of Meteorology, PO Box 243, Reading RG6 6BB, UK and 2 Observatoire de la Coˆte d’Azur, BP 4229, 06304 Nice Cedex 04, France 1 v J. I. Lunine and W. B. Hubbard 2 Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721 9 2 1 A. Burrows 0 Department of Astronomy and Steward Observatory, University of Arizona, Tucson, AZ 85721 8 9 / h ABSTRACT p - o Giantplanetsincircumstellardiskscanmigrateinwardfromtheirinitial(formation) r t s positions. Radial migration is caused by inward torques between the planet and the a disk; by outward torques between the planet and the spinning star; and by outward : v torques due to Roche lobe overflow and consequent mass loss from the planet. We i X present self-consistent numerical considerations of the problem of migrating giant r a planets. Summing torques on planets for various physical parameters, we find that Jupiter-mass planets can stably arrive and survive at small heliocentric distances, thus reproducing observed properties of some of the recently discovered extra-solar planets. Inward migration timescales can be approximately equal to or less than disk lifetimes and star spindown timescales. Therefore, the range of fates of massive planets is broad, and generally comprises three classes: (I) planets which migrate inward too rapidly and lose all their mass; (II) planets which migrate inward, lose some but not all of their mass, and survive in very small orbits; and (III) planets which do not lose any mass. Some planets in Class III do not migrate very far from their formation locations. Our results show that there is a wide range of possible fates for Jupiter-mass planets for both final heliocentric distance and final mass. Subject headings: extra-solar planets, orbital migration, mass loss – 2 – 1. Introduction The recent discoveries (Mayor & Queloz 1995; Marcy & Butler 1996; Butler & Marcy 1996; Butler et al. 1997; Cochran et al. 1997; Noyes et al. 1997) of extra-solar planets have revitalized discussions on the theory of planetary system formation and evolution. In particular, several of these planets are found to be on the order of a Jupiter-mass (1M = 2×1030g) or greater, and J in very close proximity to their central stars (Table 1). Previous to the discoveries of planetary companions in very small orbits, it was predicted that Jupiter-type planets would form (and, by implication, exist) only at or outside of the ice line (3 - 5AU) (Boss 1995). In addition, although close giant planet formation may be theoretically possible (Wuchterl 1993, 1996), this requires the initial formation of a solid core of at least 5 to 10 Earth masses which may be difficult to achieve very close to the parent star. It is therefore most likely that Jupiters cannot form at small heliocentric distances (see also Boss 1995; Guillot et al. 1996), which thus leaves the question: How did the observed massive close companions to several stars get to their current locations? Jupiter-mass planets can migrate inward from their formation locations (see, for example, Goldreich & Tremaine 1980; Ward & Hourigan 1989; Lin & Papaloizou 1986; Lin & Papaloizou 1993; Lin, Bodenheimer, & Richardson 1996). We consider under what physical conditions planets may migrate inward, and how, where, why, and when a planet will stop its inward migration. Rather than consider catastrophic interactions among more than one planet (Weidenschilling & Marzari 1996; Rasio & Ford 1996), we model isolated massive planets which smoothly migrate as a result of the net torques on the body. We find three broad classes for massive planets: (I) planets which migrate inward very quickly and disappear due to mass loss from Roche lobe overflow; (II) planets which migrate inward, lose some but not all of their mass, and stably survive at small heliocentric distances; and (III) planets which do not lose any mass during migration. Some planets in Class III move radially only a small distance from their formation locations. In the context of this picture, massive close companions can result from a variety of initial conditions, and form a population which overlaps with and includes the detected close companion population. Our distribution of final masses and heliocentric distances predicts that massive planets can be present at any heliocentric distance between their formation locations and extremely small orbits, and we predict that as detection sensitivities increase, massive planets will be found to have an almost continuous distribution of heliocentric distances, from a fraction of an AU all the way out to the ice line or beyond. Our general model is to sum up the torques on a planet and find its radial motion in the circumstellar disk. We compute, in a one dimensional model, torques on a planet due to disk interactions (inward for the parameters which we adopt); torques on a planet due to interactions with the rapidly-spinning star (outward); and torques on the planet due to mass loss onto the star (outward for the conservative mass transfer case). Sections 3, 4, and 5, respectively, describe the calculation of these torques and the resulting radial motions. These torques must be calculated numerically in order to combine the results in a self-consistent model. In Section 6, we present the results of adding these torques together and calculating the orbital evolution of massive planets. – 3 – Our work is the first to quantitatively and self-consistently consider these three torques and the problem of extrasolar planets. Lastly, in Section 7 we discuss our results. Section 8 provides conclusions and predictions. 2. Planet Formation and Assumptions in the Model Our model begins, at time = 0, with a fully formed gaseous giant planet of a given mass at a given heliocentric distance. The circumstellar disk is assumed to have smooth, power-law radial density and temperature profiles, and to be axisymmetric. Jupiter-mass planets may require most of the lifetime of the disk to accrete (106 to 107years) (Pollack et al. 1996; Zuckerman, Forveille, & Kastner 1995). However, we incorporate this possibility in our model by allowing some disks to dissipate during the planet’s migration, thus letting the effective disk lifetime (τ ) be the actual disk disk lifetime minus the time to form the planet, and consider this quantity to be the limiting time constraint. When an accreting planet has sufficient mass, it will form a gap in its disk; gap formation terminates the accretion process. In all interesting cases in this study, the gap forms quickly. Therefore, starting with a fully formed giant planet in a disk with an initial smooth density distribution is adequate and appropriate. We solve for the one-dimensional (radial) movement of massive planets. The orbits are assumed to be Keplerian and circular at all times. 3. Disk Evolution and Disk Torques We model a thin axisymmetric disk with radial temperature and initial surface density profiles given by power-laws with exponents -1/2 and -3/2, respectively (Takeuchi, Miyama, & Lin, 1996). Our nominal circumstellar disk has the following physical parameters: α equal to 5×10−3, where α is the usual disk viscosity parameter defined by α = νΩ/c2, in which ν is the viscosity in the s disk, Ω is the Keplerian angular velocity at a given heliocentric distance, and c is the sound s speed (Shakura & Sunyaev 1973); Mdisk = 1.1×10−2M⊙ (Beckwith & Sargent 1993); and scale height equal to 0.05AU at 5.2AU (Takeuchi et al. 1996). We take the nominal effective lifetime of the disk to be 107years, based on observations (Zuckerman et al. 1995), with a range of 106 to 107years. All of these physical parameters are varied from their nominal values to model disks with different properties. A massive planet and the circumstellar disk interact tidally which results in angular momentum transfer between the disk and the planet (e.g., Goldreich & Tremaine 1980; Lin & Papaloizou 1986, 1993; Ward & Hourigan 1989; Takeuchi et al. 1996; Ward 1997a; Ward 1997b). The planet’s motion in the disk excites density waves both interior and exterior to the planet. These waves create a gap in the disk as the planet clears material from its orbit. The size of the gap depends on the viscosity of the disk and inversely on the mass of the planet (Lin & Papaloizou – 4 – 1986; Takeuchi et al. 1996; see eqs. 1 - 3, below). The magnitude of the torque depends on the amount of disk material present near the planet, and thus on the size of the gap. Torques on the planet are caused by gravitational interactions between the planet and density waves which occupy Lindblad resonances in the disk. We simplify the problem by adopting the impulse approximation of Lin & Papaloizou (1986), in which dissipation of the density waves is assumed to be local and angular momentum is deposited close to the planet through dissipative phenomenon such as shocks. The local dissipation approximation is reasonable because the disk which is nearest the planet has the largest influence in causing a torque on the planet. The accuracy and validity of this approximation is discussed at the end of this section. The radial motion of a planet in a circumstellar disk, due to interactions with the disk, is given by 1 da a 2 4π Rout = − RΛ(R)Σ(R)dR (1) dt (cid:18)GM⋆(cid:19) Mp!ZRin where a is the heliocentric distance of the planet; M and M are the masses of the central star ⋆ p and the planet, respectively; R is the radial coordinate, with R and R are the inner and in out outer boundaries of the disk, respectively; Λ is the injection rate of angular momentum per unit mass into the disk due to interactions between the disk and the planet; and Σ is the surface density of the disk (Lin & Papaloizou 1986). Note that the radial motion of the planet is inversely proportional to the mass of the planet, so that more massive planets move less rapidly. The angular momentum injection rate for the impulse approximation with local dissipation is also taken from Lin & Papaloizou (1986): 4 fq2GM R ⋆ Λ(R) = sign(R−a) (2) 2R |∆p|! where f is a constant of order unity, q is the mass ratio between the planet and the star (M /M ), p ⋆ and ∆ is equal to the greater of H or |R−a|, where H is the scale height of the disk. We solve p for Σ in equation 1 with a fully implicit solution to the continuity equation for the disk rewritten as (Lin & Papaloizou 1986) ∂Σ 1 ∂ ∂ 2ΛΣR3/2 = 3R1/2 νΣR1/2 − . (3) ∂t R∂R " ∂R (GM )1/2# (cid:16) (cid:17) ⋆ We take the density equal to zero in the innermost zone of the disk to represent material from the inner boundary of the disk falling onto the star. Thus, we solve for Σ(R) and Λ(R) and then compute the integral in equation 1 to find the radial motion of the planet. The gap formed by the planet in the disk is crucial in determining the behavior of the system. Figure 1 shows gap size in the disk as a function of initial planetary mass. The gap is defined as – 5 – the region in which the surface density is less than half of what it would be if there were no planet in the disk. This figure compares our results (circles), using the local dissipation approximation of Lin & Papaloizou (1986), against results using the more computationally intensive WKB approximation of Takeuchi et al. (1996) (shown as squares). For each model, α’s of 5×10−3 and 10−2 are shown. The impulse approximation which we use closely reproduces gap formations found with the more complex scheme; therefore our simplifying assumption is valid. The time “snapshot” shown in Figure 1 is 104years, a time after the gap has fully formed but before the planet has started to migrate significantly from its initial location. The disk continues to evolve after this time, and as the planet moves inward, the gap created by the planet moves inward as well. The inner edge of the gap is continually eroded, and the outer edge fills in viscously, as the planet moves inward. Note that larger viscosities cause smaller gaps and therefore faster inward migration; larger diskmasses also cause smaller gaps and subsequently more rapid inward migration. We assume that the circumstellar disk dissipates after 107years, based on observational evidence (Zuckerman et al. 1995). In our model, we do not attempt to model the exact physics of disk dissipation. Instead, we simply assume that, as a local phenomenon, the disk has disappeared from regions close to the planet after this time. In terms of the behavior of a planet, only the local disk is important, and disk-clearing processes elsewhere do not affect the behavior of the planet. However, in future work, we intend to model disk dissipation explicitly. 4. Torques from the Spinning Star As a migrating planet gets close to its central star, tidal bulges raised on the star by the planet become important to the net radial motion of the planet. Because stars are dissipative, the stellar tidal bulge is not aligned with the line of the centers of mass of the star and the planet (Goldreich & Soter 1966; Stacey 1977; Hubbard 1984). In most cases, our model planets do not migrate to separations smaller than the co-rotation point with the rapidly rotating star; therefore, we assume that in all cases, the star is rotating faster than the orbital period of the planet, so that the tidal bulge leads the line of centers. The torque on the planet in this case is in the outward sense, and the planet slows the rotation rate of the star, as is the case with the Earth-Moon system. As energy is dissipated within the star and the star slows down, the planet must move outward to conserve angular momentum. The radial motion of the planet due to spin torque interaction with the star is given by da 9 R 4 R ⋆ ⋆ = Ω q (4) p dt 2 a Q (cid:18) (cid:19) ⋆ where Ω is the Keplerian angular velocity of the planet, R is the stellar radius, and Q is the p ⋆ ⋆ tidal dissipation factor of the star (Goldreich & Soter 1966; Lin et al. 1996). The numerical – 6 – coefficient is from the expression for body of uniform density; however, departures from this can be incorporated into uncertainties in Q (see below). Including the effects of infall of material onto ⋆ the star and the release of energy from gravitational contraction, the stellar radius decreases from early times until 107years with the following function: −1 28πσT4t 1 3 R (t)= + (5) ⋆ " GM⋆2 ! Ro3# where σ is the Stefan-Boltzmann constant, T is the surface temperature of the star (taken to be 4600K), t is elapsed time since the star completed its accretion, and R is the star’s initial radius, o taken to be four times the sun’s current radius (Cameron 1995). The dissipation factor Q describes how efficiently rotational energy is dissipated by friction within an object (Hubbard 1984; Rasio et al. 1996). We use a value of 1.5×105 for Q (Lin et al. ⋆ 1996), and a range of Q = 1.5×104 to 1.5×106, for a pre-main sequence star. By comparison, ⋆ Jupiter has Q in the range 104 to 105 (Hubbard 1984), and a main sequence star has Q ≥ 106. J ⋆ For smaller Q ’s, the star dissipates energy less efficiently, and the (outward) tidal torque on the ⋆ planet is greater than is the case for larger Q ’s. ⋆ We adopt a spin-down time for the star of 108years, since after this time, the star’s rotation rate has decreased by roughly an order of magnitude from its rapidly rotating state (Skumanich 1972). The effect of different spindown timescales, from 107 to 108years, on the results of the model is small. The tidal torque is essentially off after 108years due to the star’s higher Q and ⋆ slower rotation rate. In fact, all systems with close companions (with the exception of the tidally locked τ Boo system) are still evolving dynamically due to tidal torques between the bodies, but the orbital decay time for 51 Peg b, for example, is more than 1012years (Rasio et al. 1996), much longer than the main sequence lifetime of the star, so we do not include this late orbital evolution in our model. Tidal heating of the planet is caused by tidal bulges raised on the planet due to the star, the same mechanism which produces outward tidal torques on the planet. Significant heating of a planet due to external torques would alter its internal structure and behavior and therefore be important in the mass overflow regime. The tidal heating rate of the planet is given by Lunine & Tittemore (1993) as dE = −k2M⋆n3Rp5 21e2+ 3θ2 (6) dt a3Q 2 2 p (cid:18) (cid:19) where E is the energy input to the planet due to tidal heating, k is the planet’s Love number 2 (k ∼ 0.5 (Hubbard 1984)), n is the orbital mean motion, e is the eccentricity of the orbit, and 2 θ is the obliquity. For a planet with a four-day period and e = 0.1 and θ = 0, the heat energy input is ∼ 1026ergs/sec. Since a Jupiter-mass planet in such a close orbit would have luminosity – 7 – ≥ 1028ergs/sec (see Figure 6), we find that tidal heating of a close planetary companion is negligible. 5. Mass Loss and Conservation of Angular Momentum When a migrating planet gets sufficiently close to its primary star, the planet’s radius can exceed its Roche radius. When this occurs, mass transfer from the planet to the star takes place. During transfer, the planet moves outward to conserve the angular momentum of the system. In the case of stable mass transfer, the planet moves to a distance at which its planetary radius is equal to the Roche radius. Therefore, the distance to which the planet will move is determined by the planetary radius, and hence is a function of the intrinsic properties of the planet: its age, temperature, and mass. The planets in our model are not point masses, but have radii and internal structures that are calculated at each location and time through the planet’s evolution. 5.1. Internal Structure of Model Planets A grid of quasi-static evolution models was calculated as described in Guillot et al. (1996). For a given composition, the radius of a planet (R ) is a function of its mass M , its p p equilibrium temperature T and time t. The grid was calculated for 0.4M ≤ M ≤ 10M , eq J p J 0K ≤ T ≤ 2000K, and for a few billion years, starting from an initial extended planet of about eq 16 Jupiter radii (R ) at t = 0. J When the incoming stellar heat is redistributed equally over the entire atmosphere of the planet, the equilibrium temperature is related to the orbital distance by: L 4πσT4 = (1−A) ⋆ (7) eq 4a2 where A is the Bond albedo of the planet and L is the stellar luminosity. For simplicity, we ⋆ assume A = 0 and L⋆ = L⊙ (1L⊙ = 3.9×1033erg/sec). Of course, any uncertainty on A is equivalent to an uncertainty on L (except in the limit A → 1). The effect of those is discussed ⋆ below. The effective temperature of the planet is defined by the equation 4πR2σT4 = L +4πR2σT4 (8) p eff p p eq where L is the intrinsic luminosity of the planet, and is a function of M , T and t. The time p p eq variable is inappropriate as a conserved quantity during mass loss calculations because the interior – 8 – of the planetand theorbital evolution are notcoupled, andbecausea timeorigin cannot bedefined in an absolute way. Heliocentric distance is not a useful conserved quantity either, since a changes while the planet moves. Therefore, although the radius of the planet is, in general, R (M ,a,t), p p when considering continuous mass loss from a planet, we use R (M ,T,S), where S is the specific p p entropy (entropy per unit mass). Planetary mass, temperature, and entropy are continuous variables over which the grid of models can be both interpolated and differentiated; temperature and entropy are both conserved during mass loss events. The grid of models has minimum mass of 0.4M ; therefore, planets are considered to have lost all their mass and disappeared when J M < 0.4M . Note that none of our planets therefore ever enter a regime in which the planetary p J core is involved in mass transfer. The core would be made of solid (rocky) material, and behave very differently in a mass transfer regime; we intend to study this behavior in future work. The conservation of the specific entropy at the center of the planet during mass-loss events is exact in the case of a fully convective planet, as long as the atmosphere can adjust to the new equilibrium on a time scale much shorter than the mass loss time scale. This is generally the case since for T ≈ 1000K and pressure ≈ 1bar the heat capacity is c ≈ 108erg/Kg and the Rosseland p mean opacity is κ ≈ 10−2cm2/g, so that the radiative diffusivity can be estimated as: R 1 16σT3 K = ≈ 109cm2/s. ν ρc 3ρκ p R Therefore, the corresponding time scale is τ ≈ L2/K ≈ 0.3years, using L ≈ 108cm (about 1% ν ν of the planetary radius). This is smaller than the mass loss time scale since it takes about 103 to 104years to lose about 1% of the planetary mass by Roche lobe overflow (see below). However, planets in very close orbit do not generally stay fully convective but develop an inner radiative zone due to a strong (and unavoidable) decrease of their internal luminosity (see Guillot et al. 1996, 1997). The entropy is therefore larger at the top than at the bottom of the radiative zone. As a result, not only the atmosphere but also the inner radiative zone have to adjust to the new equilibrium. The previous time scale estimation, now estimated using twice the temperature and a mean pressure of 100bar, yields a τ which is ∼300 times larger (assuming the ν same characteristic opacity), or about 100years. This is close to the mass loss characteristic time scale: non-equilibrium effects may therefore be significant. We will neglect them in the present work, thereby somewhat underestimating the extent of mass loss by Roche lobe overflow. Other sources of uncertainty are expected to be more significant, as discussed hereafter. The largest source of uncertainty in the calculation of the internal structure and evolution of the planets and therefore on the presence and magnitude of mass loss torques is due to the inaccurate representation of the atmospheres of these objects. In our calculations, mass loss always occurs when the planet is very close to its parent star. The incoming stellar heat flux is then very significant, and the intrinsic luminosity of the planet comparatively very small. Interior and atmospheric models are coupled using a relation linking the effective temperature to the – 9 – temperature at a given pressure level (see e.g. Saumon et al. 1996). This relation is based on atmospheric models calculated assuming no incoming stellar flux, and is equivalent to assuming that the stellar flux (mostly emitted at short wavelength) is absorbed deep in the atmosphere. Although it is not totally unrealistic (and it is the best that can be done at the moment), the resulting R (M ,a,t) relation is very uncertain (i.e., the accuracy on the partial derivative dR/dt p p is probably only about a factor two), thus affecting when mass loss occurs and how much mass is effectively lost by the planet. It is also important to emphasize that the radii are defined at the 10bar level, and that “hot Jupiters” might have extended atmospheres, so that the effective radius of the planet is slightly larger than estimated here, although the amount of mass contained in the extreme upper atmosphere is small. Finally, the uncertainty in A and L will tend to shift the orbital distances at which mass ⋆ loss occurs. For example, a stellar luminosity of 2L⊙ causes a planet to lose mass at a larger heliocentric distance than for the nominal case. The result is that in the 2L⊙ case, the planet’s final mass is smaller than for the same initial mass planet in the nominal disk. However, we stress that refinements in the treatment of planetary atmospheres will not affect the qualitative behavior of the system, although the exact values of mass and distance at which mass loss occurs may change with improvements in planetary atmosphere models. 5.2. Planetary Mass Loss and Consequent Radial Motion The Roche lobe overflow regime begins at the distance where the planet’s radius is equal to or greater than its Roche lobe radius. The Roche radius is given by the following (Eggleton 1983): 0.49q2/3a R = . (9) L 0.6q2/3 +ln 1+q1/3 (cid:0) (cid:1) For stable mass transfer in which material from the planet is transferred inward onto the central star, the mass loss rate due to Roche lobe overflow of a planet is determined by and balances the net inward torque, and is given by dM dlnJ 1dlnR p = C−1M p − p (10) p dt dt 2 dt (cid:18) (cid:19) where J is the orbital angular momentum of the planet, and where p 1 1 C = (1−q)− β(1+q)+ α (11) p 2 2 (Cameron & Iben 1986; Benz et al. 1990). The derivative of the planetary radius with respect to time in equation 10 is at constant mass, but allows heliocentric distance and age to change. Since temperature is a function of distance and entropy is the time variable, we have – 10 – dlnR ∂lnR dS ∂lnR dT da p p p = + . (12) dt ∂S dt ∂T da dt The planetary radius partial differentials come from our grid of atmosphere models; dS/dt comes from the evolution models as well; dT/da comes from equation 7; and da/dt is the planet’s radial motion from the sum of disk and tidal torques. The remaining parameters in equation 11 are dlnR L β = (13) dlnq (cid:18) (cid:19)a for R from equation 9; and the parameter α , which is the exponent in the mass-radius relation, L p and should not be confused with the disk viscosity α. The mass-radius relation exponent α , p which is the derivative of planetary radius with respect to planetary mass at constant temperature and entropy, is given by dlnR p α = . (14) p dlnMp! T,S We calculate α from our grid of atmosphere evolution and structure models. Note that α in this p p work is defined differently than in Cameron & Iben (1986) and Benz et al. (1990), by a minus sign. In the mass-age-temperature regime in which we consider mass loss (planets with 0.5M < M < 10M , temperature > 1000K, and age less than 108years), a planet J p J instantaneously losing mass at constant T and S expands. In order to conserve the angular momentum of the system, a planet which loses mass to its primary must move outward. A planet losing mass is pushed outward while R > R , and the planet will lose mass until p L ′ ′ ′ ′ ′ R = R for M < M , a > a , and R > R (15) L p p p p p where the primes refer to values at the planet’s new mass and heliocentric distance. Since the Roche radius is proportional to heliocentric distance (eq. 9), any subsequent inward motion of the planet decreases the Roche radius further, which results in more mass loss from the planet. Since the circumstellar disk provides a torque which pushes the planet inward, as long as the disk is present, mass loss proceeds continuously at the planetary Roche radius, once it has started. The planet will continue to lose mass and move to larger and larger heliocentric distances to conserve angular momentum. The timescale for mass loss from the planet onto the star can be quite short, a few times 106 years or less. However, when the disk dissipates, the inward torque goes to zero, and the mass loss stops. If a planet is in the process of losing mass when the disk dissipates, mass loss will cease, and the planet will be stranded at a distance where R = R , with M < M . The L p p i model stops during the mass loss regime when the planet’s mass becomes less than 0.4 M , or J when the model stops converging during mass loss, typically as the mass approaches 0.4 M . J