Tidal Evolution into the Laplace Resonance and the
ICARUS 127, 93–111 (1997) IS965669 ARTICLE NO. Tidal Evolution into the Laplace Resonance and the Resurfacing of Ganymede ADAM P. SHOWMAN Division of Geological and Planetary Sciences 170-25, California Institute of Technology, Pasadena, California 91125, and Lunar and Planetary Institute, 3600 Bay Area Boulevard, Houston, Texas 77058 E-mail: firstname.lastname@example.org AND RENU MALHOTRA Lunar and Planetary Institute, 3600 Bay Area Boulevard, Houston, Texas 77058 Received January 10, 1996; revised December 5, 1996 v 1 / v 2 P 3/2 and v 1 / v 2 P 2 resonances can pump Ganymede’s free eccentricity up to p10 23 independent of Q9Gany /Q9J . We also show that Ganymede’s free eccentricity cannot have been produced by impact with a large asteroid or comet. 1997 Aca- We use the numerical model of R. Malhotra (1991, Icarus 94, 399–412) to explore the orbital history of Io, Europa, and Ganymede for a large range of parameters and initial conditions near the Laplace resonance. We identify two new Laplace-like resonances which pump Ganymede’s eccentricity and may help explain the resurfacing of Ganymede. Near the Laplace resonance, the Io–Europa conjunction drifts at a mean angular velocity v 1 ; 2n 2 2 n 1 , while the Europa-Ganymede conjunction drifts at a rate v 2 ; 2n 3 2 n 2 , where n 1 , n 2 , and n 3 are the mean motions of Io, Europa, and Ganymede. We find that Laplace-like resonances characterized by v 1 / v 2 P 3/2 and v 1 / v 2 P 2 can pump Ganymede’s eccentricity to p0.07, producing tidal heating several hundred times higher than at the present epoch and 2 to 30 times greater than that occurring in the v 1 / v 2 P 1/2 resonance identified previously by Malhotra. The evolution of v 1 and v 2 prior to capture is strongly affected by Q9Io /Q9J . (Here, Q9 5 Q/k is the ratio of the tidal dissipation function to second-degree Love number; the subscript J is for Jupiter.) We find that capture into v 1 / v 2 P 3/2 or 2 occurs over a large range of possible initial satellite orbits if Q9Io /Q9J # 4 3 10 24, but cannot occur for values $ 8 3 10 24. (The latter is approximately two-thirds the value required to maintain Io’s current eccentricity in steady state.) For constant Q/k, the system, once captured, remains trapped in these resonances. We show, however, that they can be disrupted by rapid changes in the tidal dissipation rate in Io or Europa during the course of the evolution; the satellites subsequently evolve into the Laplace resonance (v 1 5 v 2 ) with high probability. Because the higher dissipation in these resonances increases the likelihood of internal activity within Ganymede, we favor the v 1 / v 2 P 3/2 and 2 resonances over v 1 / v 2 P 1/2 for the evolutionary path taken by the Galilean satellites before their capture into the Laplace resonance. In addition to its surface appearance, Ganymede’s large free eccentricity (0.0015) has long been a puzzle. We find that the demic Press 1. INTRODUCTION The orbital resonances among the jovian moons Io, Europa, and Ganymede present a fascinating dynamical system. The strongest resonant interactions are those between Io and Europa and between Europa and Ganymede. The ratios of mean motions (i.e., mean orbital angular velocities) of these satellite pairs are both near 2:1, causing their successive conjunctions to occur near the same jovicentric longitude. This allows their mutual gravitational perturbations to add constructively and, as we shall see later, allows a secular transfer of energy and angular momentum from Io to Europa to Ganymede. As the ratio of mean motions is not exactly 2:1, the conjunctions between Io and Europa drift at a mean angular velocity g 1 ; 2n 2 2 n 1 , while the conjunctions between Europa and Ganymede drift at a rate g 2 ; 2n 3 2 n 2 , where n 1 , n 2 , and n 3 are the mean motions of Io, Europa, and Ganymede. The Io–Europa conjunction is locked to Io’s perijove and also to Europa’s apojove; the Europa– Ganymede conjunction occurs when Europa is near perijove. These pairwise resonances are described by the libration of the following resonance angles: u 11 5 2l 2 2 l 1 2 Ã 1 librates about 08, u 12 5 2l 2 2 l 1 2 Ã 2 librates about 1808, u 23 5 2l 3 2 l 2 2 Ã 2 librates about 08. 93 0019-1035/97 $25.00 Copyright 1997 by Academic Press All rights of reproduction in any form reserved. 94 SHOWMAN AND MALHOTRA Here l i and Ã i are the satellites’ mean anomalies and longitudes of perijove. In this paper, all subscripts i 5 1, 2, and 3 refer to Io, Europa, and Ganymede. (The Europa– Ganymede conjunction is not locked to either apse of Ganymede, so the other possible resonant variable, 2l 3 2 l 2 2 Ã 3 , circulates through all possible values.) The fourth major resonance, the Laplace resonance, is characterized by the libration of the following critical angle: f 5 2l 3 2 3l 2 1 l 1 librates about 1808. The Laplace resonance is a 1:1 commensurability between the rates of motion of the Io–Europa and Europa– Ganymede conjunctions (as opposed to the 2:1 commensurabilities between the satellites’ mean motions): the Io– Europa conjunction drifts at the same rate as the Europa–Ganymede conjunction, so that g 1 /g 2 5 1. Currently we have g 1 5 g 2 5 20.748 day 21. This is an extremely small value compared with the satellites’ mean motions, which range from approximately 508 day 21 for Ganymede to approximately 2008 day 21 for Io. These orbital resonances have a strong effect on the satellites’ thermal evolution. Io’s active volcanism and high thermal heat flux of p2 W m 22 (Smith et al. 1979, Veeder et al. 1994) are probably caused by tidal dissipation associated with its resonantly forced orbital eccentricity of 0.0044 (Peale et al. 1979). Europa’s tectonism possibly also results from tidal flexing (Malin and Pieri 1986). Although Ganymede’s eccentricity is currently too low for significant tidal heating, the ancient resurfacing on this satellite (McKinnon and Parmentier 1986) may be linked to higher tidal dissipation in the past. Especially for Ganymede, knowledge of past orbital history is critical for elucidating the thermal evolution. Yoder (1979) and Yoder and Peale (1981) constructed an analytical theory to explain the high rate of internal activity on Io as well as the origin of the Laplace resonance from initially nonresonant orbits. According to this scenario, tides raised on Jupiter by the satellites cause the satellite orbits to expand outward over time. As Io approaches the 2:1 resonance with Europa, g 1 approaches zero, forcing Io’s eccentricity, e 1 , to increase; however, tidal dissipation in Io (which increases with e 1 ) lowers Io’s semimajor axis and eccentricity. This counteracts the effects of Jupiter’s tides, which push Io outward into 2:1 resonance with Europa, as well as the resonant gravitational perturbations from Europa, which pump Io’s eccentricity. Thus, an equilibrium characterized by constant values of g 1 and e 1 is achieved, and the orbits of Io and Europa expand together while being locked in resonance. (This involves a secular transfer of orbital angular momentum from Io to Europa.) The equilibrium values of g 1 and e 1 estimated by Yoder and Peale are 21.28 day 21 and 0.0026, respectively. This is a metastable state, however, as Europa approaches the 2:1 resonance with Ganymede and the Europa–Ganymede resonant perturbations become significant. During this evolution, g 2 approaches g 1 and the satellites are captured into the Laplace resonance. Yoder and Peale calculate that the capture probability is high provided ug 1 u, ug 2 u & 28 day 21. As Jupiter’s tides continue to transfer angular momentum to the satellites (primarily Io), the Laplace resonance forces a secular transfer of orbital angular momentum from Io to Europa to Ganymede. A new equilibrium is reached in which g 1 stabilizes at a new value, and the resonantly forced eccentricities of all three satellites also reach constant values. If the present state of the system is at this equilibrium, this theory predicts Q91 /Q9J P 1.1 3 10 23. (Here Q9 ; Q/k is the ratio of the tidal dissipation function to the second-degree Love number, and subscript J refers to Jupiter.) Ganymede’s eccentricity remains low in this scenario. A second scenario was outlined by Greenberg (1982, 1987). He noted that the Yoder–Peale scenario was predicated on significant tidal dissipation within Jupiter, at a rate greater than any known physical mechanisms for tidal dissipation in gaseous planets. To circumvent this apparent difficulty, he suggested that Io, Europa, and Ganymede formed in orbits deep in resonance, with g 1 and g 2 closer to zero. Since satellite formation, dissipation in Io has decreased the satellite’s semimajor axis and increased ug 1 u; thus, Io has evolved away from the 2:1 resonance with Europa. Similarly, Europa and Ganymede were deeper in the 2:1 resonance in the past, so that Ganymede would have had a higher forced eccentricity. This scenario allows slightly more tidal heating in Ganymede than at present (eccentricity p0.003, as compared with the current free and forced eccentricities of 0.0015 and 6 3 10 24 ). However, recent theoretical work on the tidal Q of gaseous planets (Ioannou and Lindzen 1993) and estimates of low upper bounds for the tidal Q of other outer planets [Q , 39,000 for Uranus (Tittemore and Wisdom 1989) and Q , 3 3 10 5 for Neptune (Banfield and Murray 1992)] suggest that Q J was low enough for significant orbital evolution, increasing the plausibility of the tidal assembly of the Galilean resonances. Greenberg (1982) has also suggested the possibility of episodic tidal heating of Io, in which the Galilean satellites oscillate about the equilibrium point of the Laplace resonance, causing Io’s Q9 and resonantly forced eccentricity to vary periodically. This possibility was explored in some detail in Ojakangas and Stevenson (1986), and it remains a viable model for the present state of the system; however, it has not been shown to have significant import for Ganymede’s evolution. More recently, Malhotra (1991) showed that the evolutionary path described by the Yoder–Peale theory for the tidal assembly of the resonances is not unique. She found TIDAL EVOLUTION INTO LAPLACE RESONANCE that for a wide range of initial conditions, the satellites would have encountered and been temporarily captured in one or more ‘‘Laplace-like’’ resonances g 1 /g 2 P j/( j 1 1), j 5 1, 2, or 3, before evolving into the present state. (We define a ‘‘Laplace-like’’ resonance to be one in which the ratio of the mean conjunction drift rates, g 1 /g 2 , is that of two small positive integers.) Capture into any of these three resonances can occur at relatively high values of ug 1 u and ug 2 u (p7–88 day 21, before either pairwise resonance has achieved equilibrium), and is fairly likely. The 2:1 mean motion resonances then evolve in concert during passage through one of these Laplace-like resonances, as g 1 and g 2 continue to approach zero. At sufficiently small values of ug 1 u and ug 2 u, the g 1 /g 2 P j/( j 1 1) resonance is disrupted and the satellites are captured into the Laplace resonance. Of potentially great significance for Ganymede was the discovery that the g 1 /g 2 P 1/2 resonance pumps Ganymede’s eccentricity up to p0.01–0.03, possibly enough for internal activity and consequent resurfacing. For completeness, we mention here Tittemore’s (1990) proposal for the tidal heating of Ganymede. In this scenario, Europa and Ganymede pass through the pairwise 3:1 mean motion resonance which chaotically pumps up their orbital eccentricities to large values (e 2,max P 0.13, e 3,max p 0.06) before the satellites eventually disengage from that resonance. (They are presumed to subsequently evolve to their present 2:1 resonant orbits.) Tittemore argued that the extent of orbital evolution of Europa and Ganymede required in this scenario can be accommodated provided Io and Europa were locked in the pairwise 2:1 resonance early on. Tittemore’s numerical modeling of the Europa–Ganymede 3:1 resonance passage did not include the 2:1 resonant perturbations of the Io–Europa interaction, and also neglected tidal dissipation within the satellites, both factors that significantly affect the dynamical evolution of the system. It is possible to overcome these deficiencies and it would be worth reevaluating the 3:1 Europa–Ganymede resonance with a more complete numerical model; however, such a study is beyond the scope of the present work. We do not discuss this scenario further because it does not speak to the evolution of the satellites near the Laplace resonance. The three scenarios of Yoder and Peale, Greenberg, and Malhotra are best visualized by plotting the system’s path in g 1 –g 2 space. This will also prove useful for discussing our results. In Fig. 1 we depict the paths just discussed. The initial position in g 1 –g 2 space after satellite formation is completely unknown. Specifying a point on the plot is equivalent to specifying the ratios a 1 : a 2 : a 3 of the satellites’ semimajor axes. Consider the tidal assembly scenarios. Far from equilibrium of the pairwise 2:1 resonances and in the absence of any Laplace-like resonances, Io’s orbit expands much more rapidly than Europa’s, so g 1 increases faster than g 2 . Starting from its initial position, the system thus 95 FIG. 1. Several possible paths to the current state as proposed by previous authors, shown in g 1 –g 2 space. The current position is marked by a circle; the Laplace resonance, g 1 /g 2 5 1, and the Laplace-like resonance, g 1 /g 2 P 1/2, are shown by dotted lines. Dot–dashed line: Yoder and Peale (1981) scenario in which the Io–Europa 2 : 1 resonance equilibrates before capture into the Laplace resonance. Dashed line: Malhotra (1991) scenario for passage through g 1 /g 2 P 1/2 prior to capture into the Laplace resonance; in this scenario, neither 2 : 1 pairwise resonance equilibrates before capture into the Laplace resonance. Short solid line moving downward from near the origin: Greenberg (1987) scenario in which tidal dissipation in Jupiter is negligible. moves nearly horizontally to the right in g 1 –g 2 space. In the Yoder–Peale scenario, the system evolves unhindered by any Laplace-like resonance to equilibration of the Io– Europa resonance, and g 1 becomes constant (dash–dot line). As g 2 continues to increase, the system then moves vertically upward in g 1 –g 2 space. Capture into the Laplace resonance, g 1 /g 2 5 1, eventually occurs from below, i.e., from a smaller value of g 1 /g 2 . In contrast, Malhotra (1991) showed that the first metastable state in Yoder and Peale’s scenario was unlikely to be achieved as there was a high probability that the approach to this state would be interrupted by capture into a g 1 /g 2 5 j/( j 1 1) Laplace-like resonance. This is shown by the dashed line. Greenberg’s scenario is shown with a solid line. In this paper we use Malhotra’s (1991) numerical model to explore evolution into the Laplace resonance over a much wider range of conditions than she examined. Our main finding is that two other Laplace-like resonances above the Laplace resonance, g 1 /g 2 P 2 and g 1 /g 2 P 3/2, have high capture probability and can pump Ganymede’s eccentricity to p0.07. Capture into these resonances is possible only if Q91 /Q9J , 8 3 10 24. [This upper bound is 96 SHOWMAN AND MALHOTRA slightly smaller than the value needed to maintain the current orbital configuration in steady state, as discussed in Yoder and Peale (1981).] If this condition is satisfied, capture into one of these resonances is quite likely. We also find that if the Q/k are constant in time, the satellites do not evolve into the Laplace resonance from these resonances. We show that rapid changes in Q91 /Q9J or Q92 /Q9J can cause disruption followed by capture into the Laplace resonance. Several plausible mechanisms can easily produce the requisite time variability in Q9i /Q9J . The paper is organized as follows. In Section 2 we present our results. We begin with a brief description of the dynamical model and a discussion of the low-order perturbation theory for the Galilean resonances. (The details of the latter are given in the Appendix.) This is followed by a detailed description of some example runs. Next, we identify the conditions under which capture into g 1 /g 2 P 3/2 or g 1 /g 2 P 2 can occur, and we characterize the resonances by determining the eccentricities and dissipated power they produce under different conditions. We then explore the manner in which Q91 /Q9J or Q92 /Q9J must change to allow disruption of these new Laplace-like resonances and evolution into the Laplace resonance. We end the section with some additional results on the evolution of g 1 and g 2 toward equilibrium. In Section 3 we calculate the size of the cometary impactor necessary to excite Ganymede’s free eccentricity to its current value (0.0015); we show that the free eccentricity cannot have been produced by cometary impact. In Section 4 we summarize our results and conclusions. Implications for Ganymede’s thermal history are discussed in two companion papers (Showman et al. 1996, Showman and Stevenson 1996). 2. RESULTS 2.1. The Model The dynamical model we use is described in detail in Malhotra (1991). The model includes perturbations from Jupiter’s gravity field to order J4 and the mutual satellite perturbations to second order in the orbital eccentricities. The secular perturbations due to Callisto are also included. The effects of tidal dissipation in the planet as well as the satellites are parameterized by the tidal dissipation functions, Q9i , Q9J , and are included in the perturbation equations. The ratios Q9i /Q9J are free parameters which we specify as inputs to the model. As described in the Introduction, a particular solution for the evolution into the Laplace resonance admitted by this dynamical model was first found by Yoder and Peale (1981). This is not a unique solution, however, as was shown by Malhotra (1991), and there is a high probability that the evolutionary path described by the Yoder–Peale solution would be interrupted by capture into one of sev- eral possible Laplace-like resonances defined by g 1 /g 2 P j/( j 1 1), where j is an integer. These higher-order resonances arise from subtle three-body interactions that are difficult to extract analytically from the equations of motion, and the fidelity of standard perturbation theory to describe these resonances is difficult to establish. In fact, Yoder and Peale considered the possibility of one such higher-order Laplace-like resonance, g 1 /g 2 P 1/2. From their perturbation analysis, they concluded that this resonance was unstable, and passage through it would excite small free eccentricities on Io and Europa; at the lowest order, Ganymede’s eccentricity would not be perturbed by this resonance. They estimated higher order effects of this resonance on e 3 to be on the order of 10 24. Yoder and Peale were careful to note that these Laplace-like resonances are sufficiently subtle that their particular analysis may not have provided an adequate description. By taking the same perturbation approach as Yoder and Peale, we find that their analysis was incomplete in at least one respect: the g 1 /g 2 P 2 Laplace-like resonance appears in the same order in perturbation theory as the g 1 /g 2 P 1/2 resonance. Furthermore, this resonance does affect Ganymede’s orbital eccentricity in the lowest order, and is therefore of particular interest for the geophysical evolution of this satellite. We have obtained several analytical results for the g 1 /g 2 P 2 Laplace-like resonance that help to understand—and provide corroboration for—the numerical results. These calculations are given in the Appendix. Other Laplace-like resonances also are potentially significant for the problem of Ganymede, e.g., g 1 /g 2 P 3/2, but with the Yoder–Peale approach, these require going to the next order in perturbation theory. Higher-order perturbation theory in this particular context does not necessarily imply that the resonance is weak, because each new factor of the small parameter (e.g., perturbing satellite mass relative to that of Jupiter) is accompanied by at least one power of a small divisor. Many of the results of numerical simulations can be understood in light of this analysis, but not all. It appears likely to us that this particular analytical approach does not describe the complete picture for the three-body resonances, possibly because we are in a regime where the perturbation expansions are poorly convergent. (It should be emphasized that the problems lie not with the truncation of the disturbing function, but rather with the perturbative resonance analysis built on it.) A different approach is called for. We defer this to future work, and concentrate on the numerical solution of the perturbation equations. The differential equations are approximated by an algebraic mapping (details given in Malhotra 1991). This speeds up the numerical simulation by a factor of several hundred. Even so, it is necessary to artificially enhance the rate of orbital evolution to obtain results within reasonable computational time. We used Q J 5 100 in all our runs. TIDAL EVOLUTION INTO LAPLACE RESONANCE (The Q9i /Q9J are independently specified free parameters.) A typical run uses p8 hr of CPU time on an HP-735/99MHz workstation. To illustrate, suppose the ‘‘real’’ value of Q J is 10 5 ; then, using Q J 5 100 in the numerical simulation means that the entire evolution over the 4.5-byr age of the system is forced to occur over only 4.5 myr. Nevertheless, we expect that the qualitative features of the dynamics are not affected because tidal evolution with Q J 5 100 is slow enough to be adiabatic on the time scale of the gravitational perturbations. Quantitative confirmation of this assumption is discussed in Section 2.5. We note here that the model accounts for only the lowfrequency resonant perturbations, and so is valid only for sufficiently small ug 1 u and ug 2u. We restrict our calculations to ug 1 u, ug 2 u , 108 day 21. Although this range corresponds to only a few percent change in a 1 /a 2 and a 2 /a 3 , it is comparable to the extent of evolution expected over Solar System history for reasonable Q J values (p10 5 ). The satellites may thus have formed in the region of validity of the model for g 1 and g 2 . For most of the simulations discussed in this paper, Ganymede’s initial eccentricity was 0.001. Resonance encounter usually occurred (when at all) with an eccentricity somewhat smaller than the initial value. 2.2. Example Runs In this section, we describe in detail the evolution of the system in three different runs. These illustrate the range of possible orbital and dynamical histories of the Galilean satellites that we have found in more than 300 numerical simulations. The time evolution of several parameters in these runs is plotted in Figs. 2, 3, and 4. Panels (a), (b), and (c) in these figures depict the evolution of the orbital eccentricities of Io, Europa, and Ganymede; panel (d) shows the ratio g 1 /g 2 ; and (e) shows our assumed Q91 / Q9J . The time axis runs from 0 to 1.3 3 10 4 Q J years (assuming k J 5 0.38, following Gavrilov and Zharkov 1977). (In displaying the results of our simulations, we have factored out QJ on the time axis. This stresses the fact that QJ is unknown and that, for a given simulation with specified Q9i /Q9J over time, the timescale for the evolution scales linearly with QJ.) Thus, the evolution shown in these figures would occur in 4 byr if QJ 5 3 3 105, but only 400 myr if QJ 5 3 3 104. There is no special significance to the origin on the time axis. Note that the end state of the system in each of these runs is close to that observed at the present epoch: the satellites are trapped in the Laplace resonance, and the final orbital eccentricities are close to the observed forced eccentricities. The three runs differ in initial conditions and in the assumed tidal dissipation functions. Consequently, the runs differ in the sequence of Laplace-like resonances that the satellites encounter and temporarily enter before reaching the current state. 97 The initial values of the frequencies (g 1 ,g 2 ) in Figs. 2, 3, and 4 were (25.6,23.2), (26.2,22.6), and (24.7,28.0) degrees per day, respectively. All three runs begin with Q91 /Q9J 5 4 3 10 24. (Other parameter values are listed in the captions.) We also show the evolution for these runs on a g 1 –g 2 plot in Fig. 5. In Figs. 2 and 3, the initial value of g 1 /g 2 is greater than 1, whereas in Fig. 4 it is less than 1. In all three cases, g 1 and g 2 both increase toward zero over time, but g 1 increases more rapidly, so g 1 /g 2 initially decreases. In Fig. 2, g 1 /g 2 initially passes through 3/2 from above without entering this resonance. The rate of change of g 1 decreases, and at some point g 1 becomes almost constant at a value near 22.5 to 23.08 day 21 (see Fig. 5). (As in the Yoder–Peale scenario, competition between the effects of Jupiter’s and Io’s tides causes this effect, but here it occurs at a larger value of ug 1 u because Q91 /Q9J is one-third that in the Yoder–Peale scenario.) g 2 continues to increase, however, so g 1 /g 2 reverses direction and begins climbing. When g 1 /g 2 5 3/2 is reached from below, resonance capture occurs. The resonance excites Ganymede’s orbital eccentricity, and also causes large variations in Europa’s eccentricity. g 1 and g 2 continue to increase, and a new equilibrium is reached. At time t 5 10 4 Q J years, we abruptly increase Q91 /Q9J to 1.27 3 10 23. This change destabilizes the g 1 /g 2 P 3/2 resonance: a short time (p10 3 Q J years) later, the orbital parameters exhibit large fluctuations and the satellites enter the Laplace resonance (g 1 5 g 2 ). The evolution in Fig. 3 is qualitatively similar to that in Fig. 2: g 1 /g 2 initially decreases, minimizes, and then increases. The value of g 1 at the minimum is roughly the same as in Fig. 2 (see Fig. 5). Because we started with a greater initial value of g 2 , however, the minimum occurs at g 1 /g 2 . 3/2 in Fig. 3, rather than at g 1 /g 2 , 3/2 as in Fig. 2. Thus, in this case, the system never encounters the g 1 /g 2 5 3/2 resonance. When g 1 /g 2 approaches the value 2, the satellites are captured in this Laplace-like resonance. (Note again that this resonance capture occurs from below; the early encounter of g 1 /g 2 with the value 2 from above did not result in resonance capture.) This resonance also pumps up Ganymede’s eccentricity. At t 5 10 4 Q J years, we change Q91 /Q9J to 2.5 3 10 23, and the system jumps into the Laplace resonance. In our third example, shown in Fig. 4, the system first enters the g 1 /g 2 P 1/2 Laplace-like resonance; however, this resonance is soon disrupted, and g 1 /g 2 increases past 1 (without entering the Laplace resonance), and is next captured into the g 1 /g 2 P 3/2 resonance. When we change Q91 /Q9J to 1.9 3 10 23 at time at t 5 10 4 Q J years, the system jumps into the Laplace resonance. In this example there are two eccentricity pumping episodes, one for g 1 /g 2 P 1/2 and one for g 1 /g 2 P 3/2. If resurfacing is associated with each episode, two resurfacing events would occur. 98 SHOWMAN AND MALHOTRA Evolutionary paths as shown in Figs. 2, 3, and 4 represent plausible paths to the current state: in all three cases the system ends in the Laplace resonance. However, this has occurred only because we increased Q91 /Q9J by a factor of p3 during each run; we will see later that decreasing FIG. 2. First example run: The system was temporarily captured in the g 1 /g 2 P 3/2 resonance before evolving into the Laplace resonance. Shown are time evolution of the eccentricities of (a) Io, (b) Europa, and (c) Ganymede, and (d) the ratio g 1/g 2; panel (e) shows the assumed Q91 /Q9J over time. The time axes in (a)–(e) are the same and run from 0 to 1.3 3 104QJ years. In this run, the ratios of the tidal dissipation functions, Q9i /Q9J (Q9 ; Q/k), were Q93 /Q9J 5 0.127 and Q92 /Q9J 5 4.1 3 1023; Q91 /Q9J was initially set to 4 3 1024 and was changed at time 104QJ years to 1.27 3 1023, which caused the system to jump from g 1 /g 2 P 3/2 into the Laplace resonance (g 1 5 g 2 ) after rapid fluctuation of the variables. The state of the system at the end of the integration is close to that observed at present. FIG. 3. Second example run: The system was temporarily captured in the g 1 /g 2 P 2 resonance before reaching the current configuration. All panels are the same as in Fig. 2. The tidal parameters Q9i /Q9J are the same as for Fig. 2, except that at 104 QJ years, Q91 /Q9J was changed to 2.5 3 1023. This change disrupted the g 1 /g 2 P 2 resonance and allowed capture into the Laplace resonance. TIDAL EVOLUTION INTO LAPLACE RESONANCE 99 heat flow measured by Veeder et al. (1994) may require time variable Q91 /Q9J , given the lower bound on the timeaveraged Q J (Goldreich and Soter 1966, Yoder and Peale 1981). Furthermore, laboratory experiments have shown that terrestrial rocks have strongly temperature-dependent Q at high frequencies (Berckhemer et al. 1982). Although data at low frequencies are lacking, it is reasonable to expect temperature dependence at tidal frequencies also, especially at temperatures near the solidus. In addition, Io’s and Europa’s second-degree Love numbers depend on satellite structure: they are near p0.02 for a frozen interior but close to p1 for a massively molten interior. Changes in satellite internal temperature or structure could thus cause large variations in Q9i . Finally, processes may act in Jupiter to produce changes in Q J (e.g., Stevenson 1983, Ioannou and Lindzen 1993), possibly of large amplitude. Any of these mechanisms could produce time-variable Q91 /Q9J and may allow disruption of the g 1 /g 2 P 2 or g 1 /g 2 P 3/2 Laplace-like resonances followed by capture into the Laplace resonance. When Q91 /Q9J oscillates (as a sinusoid or square wave) with periods of p10 8 –10 9 years and amplitudes comparable to that used in Figs. 2–4, the Laplace-like resonances are generally disrupted near a maximum in the cycle (not necessarily the first). Subsequent minima of Q91 /Q9J , however, do not cause a second resonance capture into the Laplacelike resonance—the system generally remains in the La- FIG. 4. Third example run: The system evolved through both the g 1 /g 2 P 1/2 and g 1 /g 2 P 3/2 resonances in this run. All panels are the same as in Figs. 2 and 3. The ratios of tidal Q9i /Q9J are exactly the same as those used in Fig. 2, except for a change in Q91 /Q9J to 1.9 3 1023 at 104QJ years. The main differences are the initial values of g 1 and g 2. Q92 /Q9J by a factor of p100 also leads to disruption. Superficially, the discrete changes we make in the tidal Q’s in these numerical experiments may appear artificial to the reader; however, despite the fact that most analytical orbital modeling assumes constant Q, time variability of Q91 /Q9J and Q92 /Q9J is very likely. Indeed, the high ionian FIG. 5. Paths of the system in the three example runs of Figs. 2–4 are displayed on an g 1 –g 2 diagram. Note that all three paths change to high slope at g 1 P 238 day21. 100 SHOWMAN AND MALHOTRA place resonance itself. The Greenberg (1982)/Ojakangas and Stevenson (1986) model, in which the coupling between orbital dynamics and geophysics drives oscillation in Q91 /Q9J , thus constitutes a plausible mechanism for disruption of the g 1 /g 2 P 3/2 or 2 resonance and capture into the Laplace resonance. 2.3. Capture Statistics for the g 1 /g 2 P 3/2 and g 1 /g 2 P 2 Resonances We have seen that the g 1 /g 2 P 3/2 and g 1 /g 2 P 2 Laplace-like resonances excite Ganymede’s eccentricity to sufficiently high values that the consequent enhanced tidal heating could be geophysically significant. To estimate the viability of this scenario, we have to consider two issues: (1) What are the capture probabilities for g 1 /g 2 P 3/2 and 2 resonances assuming they are encountered? (2) What conditions allow g 1 and g 2 to evolve in such a manner that the resonances are encountered? To answer these questions, we made many numerical simulations with a large range of initial conditions. Our study covered the range (29,22) and (28,21) deg/day for the initial values of (g 1 , g 2 ), with initial g 1 /g 2 ranging from 0.2 to 3.1. We used a variety of Q91 /Q9J values and Q93 /Q9J values. Most runs used Q92 /Q9J 5 4 3 10 23 (which implies Q 2 5 100 for k 2 5 0.03, Q J 5 3 3 10 5, and k J 5 0.38); in a few runs we varied Q92 /Q9J by factors of p2. The Q9i /Q9J are constant in time for all runs discussed in this subsection. First consider capture probabilities. Of 88 runs encountering g 1 /g 2 P 3/2 from above, none was captured; however, of 64 runs encountering g 1 /g 2 P 3/2 from below, 60 were captured. For the g 1 /g 2 P 2 resonance, 0 of 36 runs were captured from above and 35 of 37 were captured from below. Thus, for both resonances, capture probabilities from below are very high, but capture from above apparently cannot occur. This behavior contrasts with that of the g 1 /g 2 P 1/2, 2/3, and 3/4 resonances, for which capture occurs from above (Malhotra 1991). This also differs from the Laplace resonance, for which capture from either above or below can easily occur. The runs encountering these resonances from above use Q91 /Q9J 5 (6–400) 3 1025, while those from below use Q91 /Q9J 5 (6–80) 3 1025. We noticed no dependence of capture probability on Q91 /Q J within that range. Next consider conditions leading to resonance encounter. From a wide variety of initial g 1 and g 2 , capture into g 1 /g 2 P 3/2 or 2 occurs commonly for Q91 /Q9J # few 3 1024, but never for Q91 /Q9J P 1023. This phenomenon results not from a different capture probability at higher Q91 /Q9J but because the resonances are never encountered from below for these large values. The runs displayed in Figs. 2, 3, and 4 suggest an explanation for this phenomenon, which we have confirmed with another p100 runs. As described by the Yoder and Peale (1981) scenario, g 1 initially increases much faster than g 2 , so the system moves almost horizontally (with low positive slope) across the g 1 –g 2 plot. Eventually, g 1 achieves equilibrium while g 2 continues to increase, so g 1 /g 2 minimizes and begins increasing, and the slope turns toward vertical (i.e., high positive slope). If Q91 /Q9J P 1023, as assumed by Yoder and Peale (1981), this happens at g 1 P 21.28 day21, as shown in Fig. 1. Unless g 2 is very close to zero, the minimum occurs at g 1 /g 2 , 1, so that the system encounters the Laplace resonance before encountering g 1 /g 2 P 3/2 or 2 from below. At these low values of ug 1 u and ug 2 u, capture into the Laplace resonance is ensured (Yoder and Peale 1981), so entry into g 1 /g 2 P 3/2 or 2 cannot occur. When Q91 /Q9J is lower, however, the equilibrium value of g 1 is more negative. This phenomenon has two effects which favor capture into g 1 /g 2 P 3/2 or 2. First, for given initial g 1 and g 2 , equilibration occurs at larger values of g 1 /g 2 than is possible for greater Q91 /Q9J . Thus, some runs minimize at g 1 /g 2 greater than 1, which is essentially impossible at Q91 /Q9J P 1023. These runs will never encounter the Laplace resonance, and will therefore be captured into g 1 /g 2 P 3/2 or 2 with near-unit probability if the minimum value of g 1 /g 2 is between 1 and 2. The fraction of initial g 1 and g 2 values for which such ensured capture occurs increases with decreasing Q91 /Q9J . Second, scenarios that encounter the Laplace resonance from below do so at larger negative values of g 1 , for which capture into the Laplace resonance is unlikely (Yoder and Peale 1981). There is thus a significant probability that the system will not enter the resonance, but instead that g 1 /g 2 will continue to increase. The system will then be captured into g 1 /g 2 P 3/2 or 2 with high probability. Our simulations confirm this picture. Of the 27 runs that minimize at g 1 / g 2 , 1 (for Q91 /Q9J P 4 3 1024), 11 entered g 1 /g 2 5 3/2 and 16 entered the Laplace resonance. (The run shown in Fig. 4 is one of the 11 that entered g 1 /g 2 P 3/2.) These runs encountered the Laplace resonance at g 1 and g 2 between 23 and 248 day21. These results are shown graphically in Fig. 6 which shows the capture probabilities into various resonances on an g 1 –g 2 diagram. For concreteness, we show probabilities for Q91 /Q9J 5 4 3 1024. The figure depicts the capture probabilities into (a) the Laplace resonance, (b) the g 1 / g 2 P 3/2 resonance, (c) the g 1 /g 2 P 2 resonance, and (d) other resonances. The shading at a given (g 1 , g 2 ) point gives the capture probability for the resonance in question, for evolution beginning at that point on the diagram. (The runs that enter g 1 5 g 2 often pass through g 1 /g 2 P 1/2, 2/3, or 3/4 first.) Dark hatching corresponds to unit probability, white to zero probability, and light hatching to intermediate (p50%) probability. The axes span approximately 0 to 258 day21. As can be seen, there are large portions of the diagram for which capture into g 1 /g 2 5 3/2 or 2 is TIDAL EVOLUTION INTO LAPLACE RESONANCE 101 FIG. 6. Schematic summary of capture probabilities as a function of initial g 1 and g 2 for the different resonances, shown on g 1 –g 2 plots. The probabilities shown hold for Q91 /Q9J P 4 3 1024. Capture probabilities from given initial (g 1, g 2) are 1 if the point is dark-hatched, 0 if the point is white, and intermediate for light hatching. The four panels show probabilities for capture into (a) the Laplace resonance, (b) g 1 /g 2 P 3/2, (c) g 1 /g 2 P 2, and (d) other states characterized by evolution to high g 1/g 2 values. The probabilities depend on Q91 /Q9J , as described in the text. ensured. More than half the area has moderate capture probability for g 1 /g 2 P 3/2 and for the Laplace resonance. Runs starting in the dark region indicated in Fig. 6d evolved toward greater g 1 /g 2 , stabilizing at g 1 /g 2 p 8–40, depending on initial conditions. We surmise that in these runs, where the system is initially very close to the Europa– Ganymede 2 : 1 resonance, equilibrium of the 2 : 1 pairwise resonances was achieved without the system ever becoming captured in a Laplace-like resonance. The probabilities shown in Fig. 6 depend on Q91 /Q9J . For smaller values of Q91 /Q9J , the pattern would be similar, but the axes would span a larger range of g 1 and g 2 and conversely for larger values of Q91 /Q9J . For example, for Q91 / Q9J p 7 3 1025, the axes span 0 to 2108 day21 in g 1 and g 2 . If Q91 /Q9J 5 1 3 1023 the picture is very different: capture into the Laplace resonance would occur from everywhere in the g 1 –g 2 diagram except possibly for very small initial values of ug 2 u; capture into g 1 /g 2 P 3/2 or 2 is impossible. The corresponding Fig. 6a would then be almost fully dark-hatched, and Figs. 6b and c would be fully white. Capture probability also depends on e 3 before resonance encounter. We conclude in the Appendix that capture into the g 1 /g 2 P 2 resonance is likely if the initial e 3 # 0.001. Because the damping time for the eccentricity is on the order of 108 years for reasonable (Q/k)3 , Ganymede’s eccentricity would likely have been negligible if resonance encounter occurred more than half a billion years after Solar System formation. In summary, for plausible values of Q91 /Q9J , capture into g 1 /g 2 P 3/2 or g 1 /g 2 P 2 is moderately probable. For a substantial fraction of plausible initial conditions, the satellites pass through the g 1 /g 2 P 1/2, 2/3, or 3/4 resonance. 2.4. Characteristics of g 1 /g 2 P 3/2 and 2 Resonances As mentioned in the Introduction, the g 1 /g 2 P 3/2 and g 1 /g 2 P 2 Laplace-like resonances provide potentially much greater tidal heating in Ganymede than the g 1 /g 2 P 1/2 resonance identified in Malhotra (1991). In Fig. 7, we display the maximum eccentricity and energy dissipation that Ganymede can achieve in the g 1 /g 2 P 3/2 and 2 resonances, and compare them with the maximum possible for the g 1 /g 2 P 1/2 resonance. As the dissipated power FIG. 7. (a) Maximum eccentricities and (b) the dissipated power associated with those eccentricities, for the g 1 /g 2 P 2, g 1 /g 2 P 3/2, and g 1 /g 2 P 1/2 resonances, as a function of Q93 /Q9J , the ratio of Q/k for Ganymede to that for Jupiter. For the g 1 /g 2 P 3/2 and 2 resonances, ‘‘maximum eccentricity’’ is the equilibrated (steady-state) eccentricity; for the g 1 /g 2 P 1/2 resonance, it is that at the peak of the largest possible eccentricity ‘‘bump.’’ (Runs for g 1 /g 2 P 3/2 and 2 used Q91 /Q9J 5 4 3 1024, and those for g 1 /g 2 P 1/2 used 1.2 3 1023. All runs used Q92 /Q9J 5 4 3 1023.) 102 SHOWMAN AND MALHOTRA depends on Q93 /Q9J , we plot the eccentricity and dissipation as a function of Q93 /Q9J . The eccentricity shown is the steady-state eccentricity occurring after equilibration of the three-body resonance. The dissipation is calculated from Q93 /Q9J , Q9J , and the equilibrated eccentricity using the formula for tidal dissipation in a homogeneous, synchronously rotating satellite in an eccentric orbit (Peale and Cassen 1978), Ė 5 21 k R 5GM 2p ne 2 , 2 Q a6 where Ė is the dissipated power, R is the satellite’s radius, a, e, and n are the orbital semimajor axis, eccentricity, and mean motion, Mp is the primary’s mass (Jupiter), G is the gravitational constant, Q is the satellite’s effective tidal dissipation function, and k is the satellite’s second-degree Love number; we use modern values a 5 1.07 3 109 m and n 5 1.0 3 1025 s21. For the highest eccentricities shown, the equilibrium eccentricity is reached 1–2 3 104 Q J years after resonance capture; however, for the lowest eccentricities, the equilibration time is only p103 Q J years. The maximum eccentricity and dissipation for g 1 /g 2 P 3/2 or 2 are significantly greater than for g 1 /g 2 P 1/2. Interestingly, as Q93 /Q9J R y, the eccentricity saturates at a finite value, and as Q93 /Q9J R 0, the eccentricity tends to zero. The dissipation requires more careful consideration. For a given Q9J , the dissipation saturates at a finite value as Q93 R 0 but tends to zero for Q93 R y; however, for a given Q93 , dissipation R y as Q9J R 0, while dissipation tends to zero for Q9J R y. Dissipation is thus not solely a function of Q93 /Q9J ; we plot it as such by incorporating Q J into the vertical axis, using k J 5 0.38. For Q J p 3 3 105, the maximum dissipation for these resonances is few 3 1011 W, about an order of magnitude lower than the primordial radiogenic heating rate (using carbonaceous chondritic radionuclide abundances for Ganymede’s rock). If the resonances are disrupted before equilibration occurs, the actual peak eccentricity and dissipation would be lower, as occurs, for example, in Fig. 4 during the g 1 /g 2 P 3/2 resonance. We find that the dissipation and eccentricity for g 1 /g 2 P 3/2 are roughly independent of g 1 or g 2 at capture. For a given Q93 /Q9J , runs achieving values of ug 1 u from 2 to 88 day21 at capture and from 1.5 to 4.38 day21 at equilibrium of the three-body resonance yielded steady-state eccentricities differing by only 4%. Similar results hold for g 1 /g 2 P 2. In contrast, Malhotra (1991) found that for g 1 /g 2 P 1/2 the maximum eccentricity depends strongly on g 1 at capture. This occurs because entry into the resonance at lower values of g 1 and g 2 allows a longer resonance lifetime, which leads to higher eccentricities; however, the g 1 / g 2 P 3/2 and 2 resonances are never disrupted by themselves, so the lifetimes are not short enough to prevent equilibration of the eccentricity. The power and eccentric- TABLE I Disruption from v1 / v 2 P 3/2 into Other Resonances Time of disruption (QJ years) (Q91 /Q9J )0 3 3 103 7 3 103 104 1.3 3 104 5 3 1024 8 3 1024 1.3 3 1024 No No g1 5 g2 No No g1 5 g2 No No g1 5 g2 No No g 1 /g 2 P 2 Note. Entries display the resonance entered after disruption. Runs labeled ‘‘no’’ were not disrupted from g 1 /g 2 P 3/2. ity for g 1 /g 2 P 1/2 shown in Fig. 7 are the maximum found by Malhotra; capture at other g 1 can produce a power for g 1 /g 2 P 1/2 up to p10 times lower. Two runs (with Q93 /Q9J P 1021 –1022) entering g 1 /g 2 P 3/2 yielded anomalously low eccentricities of p0.008. These values fall far below the curve in Fig. 7. We do not understand this phenomenon, but the probability appears small, since it occurred in only 3% of our runs for g 1 /g 2 P 3/2 and never occurred for g 1 /g 2 P 2. 2.5. Disruption of g 1 /g 2 P 3/2 and 2 Resonances For constant Q9i /Q9J values, none of the runs that entered g 1 /g 2 P 3/2 or 2 was ever disrupted from the resonance. Thus, if the Q/k are constant in time, the g 1 /g 2 P 3/2 and 2 resonances cannot be possible paths to the current state. We find that by changing Q91 /Q J or Q92 /Q9J during the run, however, the resonances can be disrupted and the satellites subsequently evolve into the Laplace resonance. We did a number of runs to characterize the conditions leading to disruption. First consider the g 1 /g 2 P 3/2 resonance. This set of runs was started with g 1 5 25.68 day21 and g 2 5 23.28 day21, with initial g 1 /g 2 5 1.8. The system thus encountered the g 1 /g 2 P 3/2 resonance from above, passed through it without capture, minimized at g 1 /g 2 Q 1.2, and as g 1/g 2 increased to 1.5, it entered the g 1 /g 2 P 3/2 resonance from below at g 1 5 22.58 day21. We used Q93 /Q9J 5 0.0127. The system entered resonance at t 5 2300Q J years, and we continued the integration until 1.7 3 104 Q J years (5 byr for Q J 5 3 3 105). During this initial phase of the run, we kept Q91 /Q9J constant at 4 3 1024. Then, at a time t 0 , we changed Q91 /Q9J to another value, (Q91 /Q9J)0 . We used several different values of t 0 [3300Q J , 6600Q J , 104 Q J , and 1.33 3 104 Q J years] and of (Q91 / Q9J)0 [5.0 3 1024, 7.6 3 1024, and 1.3 3 1023]. The results of these 12 runs are summarized in Table I. None of the runs for the smaller two values of (Q91 /Q9J)0 were disrupted, and all four runs for the larger value were disrupted. The conditions for disruption thus seem relatively insensitive to the time during resonance at which disruption occurs. Of the four runs that were disrupted, TIDAL EVOLUTION INTO LAPLACE RESONANCE three were captured into the Laplace resonance and one into g 1 /g 2 P 2. Disruption usually occurred p103 Q J years after t 0 rather than immediately; such delays are evident in Figs. 2–4 for our example runs. (When the system is disrupted from resonance, the new path taken is not predictable on a case-by-case basis. Thus, we cannot infer from Table I that capture probabilities into g 1 5 g 2 vs g 1 / g 2 P 2 vary with time of disruption from g 1 /g 2 P 3/ 2. Rather, our runs suggest only that disruption into the Laplace resonance is more likely than into g 1 /g 2 P 2 in an average sense.) We also tried two runs for smaller values of Q93 /Q9J , 1.27 3 1023 and 1.27 3 1024, with (Q91 /Q9J)0 5 1.3 3 1023 and Q92 /Q9J 5 4.1 3 1023. The time between t 0 and disruption increased substantially over the case with Q93 /Q9J 5 0.0127. The first run was disrupted, but the second run did not disrupt before the end of the simulation. The resonance appears to be slightly more stable at low Q93 /Q9J , presumably because of the lower eccentricities. Using a greater value (Q91 /Q9J)0 5 1.9 3 1023 allows disruption of the resonance almost immediately at these low Q93 /Q9J values, however. Thus, the dependence on Q93 /Q9J seems much weaker than that on (Q91 /Q9J)0 . We also made four runs for g 1 /g 2 P 3/2 in which we kept Q91 /Q9J constant over the run, but in which we increased or decreased Q92 /Q9J by one or two orders of magnitude at time t 0 . Our initial value was Q92 /Q9J 5 4.1 3 1023 (same as that in the runs discussed above). We found that abruptly increasing Q92 /Q9J by one or two orders of magnitude during the run produced no visible effect on the resonances. Decreasing by one order of magnitude destabilized but did not disrupt the resonance, and decreasing by two orders of magnitude knocked the system into the Laplace resonance. Next consider g 1 /g 2 P 2. We find that this resonance is somewhat more stable than g 1 /g 2 P 3/2, requiring an increase in Q91 /Q9J to (Q91 /Q9J)0 5 1.9 3 1023 for disruption (using the same Q93 /Q9J and Q92 /Q9J as the 12 runs in Table I). Changing Q92 /Q9J (but not Q91 /Q9J) during the runs had the same effect as for the g 1 /g 2 P 3/2 resonance: increasing Q92 /Q9J by one or two orders of magnitude had no effect, decreasing it by one order of magnitude destabilized but did not disrupt the resonance, and decreasing it by two orders of magnitude disrupted the resonance, with subsequent capture into the Laplace resonance. Both resonances are stable to large decreases in Q91 / Q9J from 4 3 1024 to 3 3 1025. The change simply shifted the three-body equilibrium from 21.5 to 288 day21. Finally, we mention two runs within the Laplace resonance, during which we changed Q91 /Q9J to the large value p1022. In these runs, Io’s eccentricity reached a maximum value of 0.012 after the change, soon followed by disruption of the Laplace resonance. The system then settled at g 1 / g 2 P 9/10. This result is consistent with the Yoder–Peale 103 theory, which predicts the Laplace resonance to be unstable for e 1 . 0.012. The above results can be understood as follows in light of the analysis given in the Appendix. Within any resonance, the equilibrium value of g 1 (and hence also the equilibrium value of Io’s forced eccentricity, e 12 ) is determined by the value of D1 Q 4.85(Q9J /Q91); e 12 is proportional to Q91 /Q9J [see Eqs. (A49) and (A60) in the Appendix]. For Q91 /Q9J 5 4 3 1024, we find from Eq. (A60) that e 12 Q 0.0020 in the g 1 /g 2 P 2 Laplace-like resonance; this is verified in the numerical simulations. (We have also estimated the equilibrium value for the g 1 /g 2 P 3/2 Laplacelike resonance: e 12 P 0.0023.) Note that the three-body resonant interactions also force large-amplitude oscillations about the mean forced eccentricity, so that the maximum eccentricity is significantly larger. When we change the value of Q91 /Q9J during the evolution, the equilibrium is shifted, and the system tries to evolve to the new equilibrium while maintaining the three-body resonance lock. A decrease in Q91 /Q9J shifts e 12 and g 1 to smaller values, and the system accordingly moves to the lower equilibrium point, with no loss of stability. An increase in Q91 /Q9J , on the other hand, raises these equilibrium values, and the system tries to move toward this higher equilibrium. In this case, there are several sources of instability. 1. The resonances are unstable when e 1 exceeds some maximum value. For the Laplace resonance, Yoder and Peale (1981) estimated the upper limit to be 0.012 from a linear stability analysis. For the other Laplace-like resonances, we expect a smaller upper limit for stability. 2. A second, qualitatively distinct source of instability is the overlap of all Laplace-like resonances sufficiently close to the origin in the g 1 , g 2 plane. The numerical explorations by Malhotra (1991) as well as those in the present work indicate that sufficiently close to the origin in the g 1 , g 2 plane, the phase space is dominated by the Laplace resonance. The Galilean satellites evolving along one of the Laplace-like resonances, g 1 /g 2 P j/( j 1 1) or ( j 1 1)/j, eventually exit that resonance when both ug 1 u and ug 2 u become sufficiently small; after a period of chaotic evolution, they eventually enter the Laplace resonance. 3. We see from Eqs. (A59) and (A60) that, after Io’s eccentricity has reached equilibrium, an increase in (Q91 / Q9J) by a factor greater than 55.2/21.0 Q 2.6 changes the sign of the tidal term in the pendulum-like equation for the g 1 /g 2 P 2 Laplace-like resonance. (This factor is slightly different for other Laplace-like resonances.) This means that the tidal torque on the equivalent pendulum would change direction, causing a slow increase of libration amplitude and eventual disruption of the resonance librations. All of the above three sources of instability come into play when we increase Q91 /Q9J during the course of the evolution. Which of them is the primary cause of resonance disruption 104 SHOWMAN AND MALHOTRA depends on the exact parameters and state of the system. But it is clear that even in the absence of the first two, a sufficiently large increase in Q91 /Q9J during the course of the evolution would certainly disrupt the resonance. In the numerical experiments where we start with Q91 / Q9J 5 4 3 1024, and later increase its value by a factor of about 3, we note that the g 1 /g 2 P 3/2 or 2 Laplace-like resonance is not immediately disrupted (cf. Figs. 2–4). Io’s eccentricity increases toward the new equilibrium value, while the ratio g 1 /g 2 is, on average, maintained at the resonant value. During this period, Ganymede’s forced eccentricity is also maintained at a high value, although it does not continue its previous increase; this is further evidence of the fact that the Laplace-like resonance remains in place with the new value of Q91 /Q9J , but there is no longer the same rate of secular transfer of orbital energy and angular momentum to Ganymede from Io and Europa. Eventually the g 1 /g 2 P 3/2 or 2 resonance does become unstable, when Io’s mean forced eccentricity approaches p0.005 (maximum e 1 exceeds p0.01) This is followed by a short period, p103 Q J years, of chaotic, large-amplitude fluctuations in all dynamical variables, culminating in the satellites’ entering the Laplace resonance. The Laplace resonance cannot maintain Ganymede’s high forced eccentricity, so e 3 plummets rapidly. It appears likely in this case that the parameters of the system conspire in such a way that all three causes of instability identified above occur nearly simultaneously. In the above numerical experiments, we have chosen to use a very simple step-function model for the time variation of Q91 /Q9J . This is admittedly not a physically realistic model; however, physically realistic models for the Q and k 2 of Io as well as the Q of Jupiter remain highly unconstrained at present. Given this, building a specific physical model for the evolution of these parameters and folding it in with the orbital dynamics model are poorly justified. The numerical experiments here do serve our limited purpose of investigating the response of the three-body Laplace-like resonances to changes in Q91 /Q9J . 2.6. Other Results We report in this section the results of a few runs in which the system entered other resonances. One entered g 1 /g 2 P 13/6, and another entered what was perhaps g 1 / g 2 p 7/3 or 11/5. Two runs spent a few percent of the integration time in g 1 /g 2 P 4/3 before disruption into another resonance (g 1 /g 2 P 2 in one case and g 1 /g 2 p 7/3 in another). None of these resonances pumps Ganymede’s eccentricity, and the combined probabilities for these paths appear low. These resonances are dynamically similar to the g 1 /g 2 P 3/2 and 2 resonances: capture occurred only from below, and, with the exception of g 1 /g 2 P 4/3, disruption never occurred; the system equilibrated at values of g 1 greater than that for the Io–Europa pairwise resonance alone. In contrast, resonances below g 1 /g 2 5 1 are always disrupted, usually into the Laplace resonance or g 1 /g 2 P 3/2 (depending on Q91 /Q9J); these resonances never achieved equilibrium. The resonances above the Laplace resonance thus seem to exhibit dynamics very different from those below the Laplace resonance. 3. GANYMEDE’S FREE ECCENTRICITY In this section, we consider the problem of Ganymede’s large free eccentricity, e free 5 0.0015. This has long been a puzzle because the tidal damping time for Ganymede’s eccentricity is 106 Q 3 , or about 108 years for a plausible Q 3 P 102 ; its free eccentricity should therefore have long damped by the present time. Although free eccentricity is often portrayed as a ‘‘remnant’’ primordial eccentricity, orbital resonances can pump the free as well as the forced eccentricity. (The free eccentricity is manifested in Figs. 2–4 as a rapid oscillation about the mean value.) Since the Laplace resonance does not pump Ganymede’s free eccentricity, e free may thus provide a useful constraint on past orbital evolution. We find that the free eccentricity Ganymede attains in the g 1 /g 2 P 3/2 or 2 resonance is greater for smaller Q91 /Q9J , possibly because smaller Q91 / Q9J leads to equilibration of the three-body resonances closer to the (g 1 , g 2 ) origin. For constant Q91 /Q9J P 4 3 1024, the g 1 /g 2 P 3/2 and 2 resonances pump Ganymede’s free eccentricity to p6–9 3 1024, independent of Q93 /Q9J . Although this is almost an order of magnitude higher than that typically achieved during g 1 /g 2 P 1/2, it is still a factor of p2 lower than that observed. We consider here the size of an asteroidal or cometary impactor required to excite Ganymede’s eccentricity to 0.0015, and compare this estimate with the impactor masses obtained from Ganymede’s largest craters. Let the comet strike Ganymede at an angle u tangent to Ganymede’s orbit, where u is measured in an inertial reference frame centered on Jupiter. Let u p 0 imply an overtaking collision and u p f imply a head-on collision. We assume the comet strikes through Ganymede’s center of mass, so no change in rotation occurs. We further assume an inelastic collision in which Ganymede retains the cometary mass. If the eccentricity is initially zero, then the eccentricity after impact is (to second order in m/M, which is taken to be a small quantity) e2 5 S DF m M 2 41 G v2c vc 2 cos u , 2 (1 1 3 cos u ) 2 8 vi vi where m and vc are the cometary mass and speed, and M and vi are Ganymede’s mass and speed, respectively. Both TIDAL EVOLUTION INTO LAPLACE RESONANCE speeds are measured relative to Jupiter. For vc 5 vi and u 5 0, the comet should ‘‘soft land’’ on Ganymede without perturbing the eccentricity. As required, the equation gives e 5 0 for this situation. A comet capable of pumping Ganymede’s eccentricity to 0.0015 must thus be of mass p1020 kg. Is this plausible? The largest fresh crater on Ganymede is Gilgamesh, which has a probable excavation diameter of 550 km (Shoemaker et al. 1982). We calculate the transient crater diameter using the relation from McKinnon and Schenk (1995), and then estimate the impactor mass from the scaling law in Chapman and McKinnon (1986, p. 502), m5 F G S D 3V 4f r 3 4fA 3/(32a) 3.22g u2 3a /(32a) , where V is the transient crater volume, r is the impactor density (assumed equal to the target density), u is impact speed, g is surface gravity, and A and a are dimensionless parameters. Using nominal values u 5 15 km sec21, g 5 1.5 m sec22, plus A 5 0.19 and a 5 0.65 appropriate to a solid ice target, we find that a 1017-kg bolide created Gilgamesh. Choosing smaller values for A and a allows a larger impactor for a given crater size. Using A 5 0.1 and a 5 0.5, which are plausible lower bounds, implies an impactor of mass p1018 kg, still far too small. Other large impact features on Ganymede, such as palimpsests, are so old that any eccentricity they produce would be damped by the present time. Thus, we conclude that Ganymede’s eccentricity cannot have been produced by an impact. Although the g 1 /g 2 P 3/2 and 2 resonances excite Ganymede’s free eccentricity to only half the current value in our runs, paths may exist that allow larger free eccentricities. Exploratory runs suggest that these resonances are slightly more stable to gradual billion-year changes in Q91 /Q9J than to abrupt changes (although disruption can still occur). Slow change to Q91 /Q9J 5 1023 leads to e free 5 1.0 3 1023 in cases when disruption does not occur, and larger e free may be possible if Q91 /Q9J rises even higher before disruption. 4. SUMMARY We have used the model of Malhotra (1991) to explore the orbital dynamics and tidal evolution of Io, Europa, and Ganymede near the Laplace resonance. Our principal results are those relevant to Ganymede’s thermal history. We have shown that if the satellites passed through a Laplace-like resonance characterized by either g 1 /g 2 P 3/2 or g 1 /g 2 P 2, Ganymede’s eccentricity could have risen as high as p0.07. These resonances produce a tidal dissipation rate several hundred times times higher than at the present epoch (for reasonable Q J ), and p2–30 times higher than 105 the maximum possible for the g 1 /g 2 P 1/2 Laplace-like resonance identified by Malhotra (1991). Capture probabilities for these two resonances are p0.9 if they are encountered from below (i.e., with g 1/g 2 increasing with time); the capture probability is negligible otherwise. We have found that resonance capture can occur with initial conditions in a substantial fraction of the g 1 –g 2 plane (i.e., for a large range of initial a 1 /a 2 and a 2 /a 3 values) provided Q91 /Q9J # 6 3 1024. This upper limit is slightly smaller than that needed to maintain Io’s current eccentricity in steady state, Q91 /Q9J P 1.1 3 1023 (cf. Yoder and Peale 1981). For the latter value, the properties of the evolution of g 1 and g 2 conspire with capture statistics for the Laplace resonance to prevent the system from ever encountering these new resonances from below. This explains why they were not seen by Malhotra (1991), who restricted her study to Q91 /Q9J P 1.1 3 1023. We note that the smaller values of Q91 /Q9J (which allow capture in the g 1 /g 2 P 3/2 and g 1 /g 2 P 2 resonances) are reasonable, and could in fact be close to the modern value because that required for equilibrium at present may be too large to maintain Io’s current heat flow (Veeder et al. 1994). These new resonances differ from g 1 /g 2 P 1/2 in one important respect: they achieve stable equilibria and, therefore, are never disrupted if the Q/k are constant in time; however, we find that increasing Q91 /Q9J by a factor of p3 or decreasing Q92 /Q9J by a factor of p100 disrupts the resonances, allowing capture into the Laplace resonance with high probability. We have verified many of the results of our numerical simulations by a perturbation analysis of the low-order resonant interactions. Just as in the previous study of Malhotra (1991), the net orbital expansion of the Galilean satellites required in these evolutionary paths is easily accommodated with a tidal dissipation function for Jupiter, Q J , of a few 3105. In two companion papers we discuss the implications of resonance passage for Ganymede’s thermal history. In Showman et al. (1996) we couple the orbital model to an internal model for Ganymede. In certain circumstances, nonlinear ‘‘thermal runaways’’ can occur within Ganymede, causing internal warming and melting of the ice mantle. In Showman and Stevenson (1996), we propose models of local, near-surface thermal runaways, and evaluate the efficacy of several resurfacing mechanisms. Because the greater dissipation increases the likelihood of resurfacing, we favor g 1 /g 2 P 3/2 or 2 over g 1 /g 2 P 1/2 for the orbital history of the Galilean satellites before evolution into the Laplace resonance. Since the current resonances do not explain Ganymede’s free eccentricity, e free is a remnant of Ganymede’s prior history. We have shown that the mass of a cometary or asteroidal impactor required to pump e free to its current value is p102 –103 times larger than that which produced Gilgamesh, the only candidate crater for such an impact. 106 SHOWMAN AND MALHOTRA We surmise that the current free eccentricity is thus remnant from an ancient resonance passage. The g 1 /g 2 P 3/ 2 and 2 resonances can pump e free up to two-thirds its modern value. Our models employ constant or step-function Q9i /Q9J , however. More realistic time variation in these parameters may yield larger e free . Other observed orbital parameters of the Galilean satellites also hold memory of their prior evolution and may provide further constraints. The libration amplitudes for the u 11 , u 12 , u 23 , and Laplace resonances (Sinclair 1975) contain information about the age of the resonances, the g 1 –g 2 path followed in the past, the Q9i /Q9J and their time histories, and other factors. A successful model of the orbital and thermal history of these satellites should, in principle, yield simultaneously the current values of all four libration amplitudes; this constraint could allow some orbital histories to be excluded. In practice, however, the factors affecting the amplitudes might be difficult to separate, preventing useful constraints on any individual parameter from being developed. Furthermore, passage through resonance inevitably involves the crossing of a chaotic region of phase space. This means that a certain degree of uncertainty in the final values of the orbital parameters is unavoidable. Finally, it is also worth keeping in mind that for a dissipative system, there are multiple, nonunique, paths to an equilibrium (or quasi-equilibrium) state; the final state of the system would retain only partial memory of initial conditions and intermediate states. tricity. [We also take account of three other sources of secular terms: in the Io–Ganymede interaction, the secular perturbations of Callisto on Io, Europa and Ganymede, and the effects of the oblateness of Jupiter. Details on these are not essential to the analysis here; the interested reader is referred to Malhotra (1991).] The expansion of R (12) is given by R (12) 5 f0(e 21 1 e 22 ) 1 f1 e 1 cos(u 1 2 Ã 1 ) 1 f2 e 2 cos(u 1 2 Ã 2 ) 1 f3 e 21 cos 2(u 1 2 Ã 1 ) 1 f4 e 22 cos 2(u 1 2 Ã 2 ) (A1) 1 f5 e 1 e 2 cos(2u 1 2 Ã 1 2 Ã 2 ) 1 f6 e 1 e 2 cos(Ã 1 2 Ã 2 ). The expression for R (23) is obtained by shifting the satellite indices by 1 (1 R 2, 2 R 3). Here Ã i are the longitudes of periapse; u 1 5 2l 2 2 l 1 and u 2 5 2l 3 2 l 2 are the pairwise 2 : 1 resonant combinations of mean longitudes l i ; and the fi are functions of Laplace coefficients, b s( j)(a) (a 5 a 1 /a 2 , 1 is the ratio of the semimajor axes). The expressions for the fi and their numerical values evaluated at exact resonance, a 5 (1/2)2/3 5 0.62996, are as follows: D d2 1 d (0) 1 a 2 2 b 1/2 2a (a) Q 10.39, 8 da da f1 5 2 1 d 41a b (2) (a) Q 21.19, 2 da 1/2 f2 5 1 1 d 31a b (1) (a) 2 2a Q 10.43, 2 da 1/2 f3 5 1 d2 1 d (4) 44 1 14a (a) Q 11.67, 1 a 2 2 b 1/2 8 da da f4 5 1 d2 1 d (2) 38 1 14a (a) Q 13.59, 1 a 2 2 b 1/2 8 da da f5 5 2 d2 1 d (3) 42 1 14a (a) Q 24.97, 1 a 2 2 b 1/2 4 da da f6 5 1 d2 1 d (1) 2 2 2a (a) Q 10.58. 2 a 2 2 b 1/2 4 da da APPENDIX In this Appendix we give a perturbative analysis of the evolution of the Galilean satellites near the current pairwise 2 : 1 mean motion resonances, the Laplace resonance and the lowest-order Laplace-like resonances. We generally follow the approach of Yoder and Peale (1981), so that the first part of the analysis is essentially similar to theirs; we provide it here only for completeness and easy access for the subsequent calculations. We carry the analysis further to isolate the g 1 /g 2 P 2 Laplace-like resonance that pumps up Ganymede’s eccentricity. Other Laplacelike resonances discussed in this paper require even higher order analysis, and we do not pursue that here. For the mutual perturbations of the inner three Galilean satellites, Io, Europa, and Ganymede, we have two disturbing functions, (Gm 1 m 2 /a 2 )R (12) and (Gm 2 m 3 /a 3 )R (23), where G is the universal constant of gravitation, m i are the satellite masses, and a i are the orbital semimajor axes; the indices 1, 2, and 3 refer to Io, Europa, and Ganymede, respectively. We use a series expansion of the pairwise disturbing function in powers of the orbital eccentricities, retaining only the terms that describe the secular and 2 : 1 resonant interactions between Io and Europa and between Europa and Ganymede, up to the second order in eccen- S S S S S S S f0 5 1 D D D D D D (A2) The equations for the perturbations are most conveniently written in terms of the Poincaré variables, (k i , h i ) 5 (e i cos Ã i , e i sin Ã i ), and the angles, u 1 and u 2 . Then, k̇ 1 5 2e 2 an 1 (12) R , h 1 ḣ 1 5 1e 2 an 1 (12) R , k 1 (A3) 107 TIDAL EVOLUTION INTO LAPLACE RESONANCE k̇ 2 5 2e 1 n 2 (12) (23) R 2 e 3 an 2 R , h 2 h 2 ḣ 2 5 1e 1 n 2 (12) (23) R 1 e 3 an 2 R , k 2 k 2 k̇ 3 5 2e 2 n 3 (23) R , h 3 (23) R , ḣ 3 5 1e 2 n 3 k 3 ṅ 1 5 13e 2 an 21 (12) R , u 1 (12) (23) ṅ 2 5 26e 1 n 22 R 1 3e 3 an 22 R , u 1 u 2 ṅ 3 5 26e 2 n 23 (23) R , u 2 variations. We thus obtain the following first order ‘‘forced’’ variations of k i , h i : SD S D de 2i dt T T 5 2c i (1 2 (7Di 2 12.75)e 2i ), 5 2Sd c i (7Di 2 4.75)e 2i . (A15) d1 k3 5 e32 cos(u2 1 «3), (A5) d1 h3 5 e32 sin(u2 1 «3), (A16) (A6) where, to lowest order, we have (A7) (A8) u 1 5 g 1 t 1 u 10 , (A17) u 2 5 g 2 t 1 u 20 , (A18) g 1 5 2n 2 2 n 1 (A19) g 2 5 2n 3 2 n 2 (A20) and are the unperturbed frequencies of u 1 and u 2 . The forced eccentricities are given by (A9) e 12 5 n 1 e 2 af1 , g1 2 g1 (A21) (A10) e 21 5 n 2 e 1 f2 , g1 2 g2 (A22) e 23 5 n 2 e 3 af1 , g2 2 g2 (A23) e 32 5 n 3 e 2 f2 , g2 2 g3 (A24) k̇ 1 5 2Jd c 1 D1 k 1 2 g 1 h 1 2 n 1 e 2 af1 sin u 1 , (A11) The phase lags in Eqs. (A14)–(A16) are given by k̇ 2 5 2Jd c 2 D2 k 2 2 g 2 h 2 2 n 2 e 1 f2 sin u 1 2 n 2 e 3 af1 sin u 2 , ḣ 2 5 2Jd c 2 D2 h 2 1 g 2 k 2 1 n 2 e 1 f2 cos u 1 1 n 2 e 3 af1 cos u 2 , (A12) k̇ 3 5 2Jd c 3 D3 k 3 2 g 3 h 3 2 n 3 e 2 f2 sin u 2 , ḣ 3 5 2Jd c 3 D3 h 3 1 g 3 k 3 1 n 3 e 2 f2 cos u 2 , (A14) d1 h2 5 e21 sin(u1 1 «21) 1 e23 sin(u2 1 «23), Now, the lowest-order equations for the (k, h) are as follows: ḣ 1 5 2Jd c 1 D1 h 1 1 g 1 k 1 1 n 1 e 2 af1 cos u 1 , d1 h1 5 e12 sin(u1 1 «1), d1 k2 5 e21 cos(u1 1 «21) 1 e23 cos(u2 1 «23), where the n i are the mean motions of the satellites, and e i 5 m i /M are the satellite masses relative to the mass of Jupiter. Note that the lowest-order terms for k̇ i , ḣ i are O (e), whereas for ṅ i , the lowest-order terms are O (ee). Following Yoder and Peale [1981, their equations (3)– (6)], we include the effects of tidal dissipation according to ṅ i ni d1 k1 5 e12 cos(u1 1 «1), (A4) (A13) where g i are the precession rates induced by the secular perturbations alone; in Yoder and Peale’s (1981) notation, g i 5 Ã̇ is . We can obtain an approximate solution of these equations by assuming that u̇ 1 5 2n 2 2 n 1 and u̇ 2 5 2n 3 2 n 2 are slowly varying quantities and neglecting their time «1 5 7 c 1 D1 , 3 g1 2 g1 (A25) « 21 5 7 c 2 D2 , 3 g1 2 g2 (A26) « 23 5 7 c 2 D2 , 3 g2 2 g2 (A27) «3 5 7 c 3 D3 . 3 g2 2 g3 (A28) The above approximate solution is valid in the regime where ug 1 u and ug 2 u are not too small compared with 108 SHOWMAN AND MALHOTRA pn(ee)1/2, a condition that is satisfied for most of the range of evolution of the system that we have explored in this paper. To obtain the lowest-order perturbations on the mean motions, we retain only the first-order eccentricity terms in the disturbing function and substitute the forced solutions [Eqs. (A14)–(A16)] into the right hand side of Eqs. (A6)– (A8). We also include the dissipative terms [Eq. (A9)]. This yields the following equations that describe the lowest-order three-body interaction: ṅ 1 5 13n 21 e 2 a[ f1 e 12 « 1 1 f2 e 21 « 21 1 f2 e 23 sin(f 1 « 23 )] 1 (ṅ 1 )T , (A37) ṅ 2 5 2A2 n 22 sin f 2 10.31n 1 c 1 D1 e 212 , (A38) ṅ 3 5 2A3 n 22 sin f, (A39) where A1 5 23e 2 e 3 a 21f1 f2 A2 5 13e 1 e 3 af1 f2 (A29) ṅ 2 5 26n 22 e 1[ f1 e 12 « 1 1 f2 e 21 « 21 1 f2 e 23 sin(f 1 « 23 )] 1 3n 22 e 3 a[ f1 e 23 « 23 1 f2 e 32 « 3 2 f1 e 21 sin(f 2 « 21 )] 1 (ṅ 2 )T , (A30) ṅ 3 5 26n 23 e 2[ f1 e 23 « 23 1 f2 e 32 « 3 2 f1 e 21 sin(f 2 « 21 )] 1 (ṅ 3 )T , ṅ 1 5 2A1 n 22 sin f 2 n 1 c 1(1 2 14D1 e 212 ), r5 (A32) 1 f2 e 21 « 21 ] SD (A33) a (A42) n 2 e 2 f2 e 21 g1 2 g2 (A43) ġ 1 5 2(2A2 2 A1 )n 22 sin f 1 n 1 c 1(1 2 34.6D1 e 212 ), (A44) sin f 1 (A45) A2 )n 22 10.3n 1 c 1 D1 e 212 , and the dynamical equation for f is f̈ 5 2(A1 2 3A2 1 2A3 )n 22 sin f 2 n 1 c 1(1 2 44.6D1 e 212 ). (A46) Q 7n 1 c 1 D1 e 212 1 9.50n 2 c 2 D2 e 221 , where we have set n 1 /n 2 Q 2, a 5 0.63, and used the numerical values of e i . Similar simplifications follow for the other terms. We then obtain ṅ 1 5 2n 1 c 1[1 2 (7e 21 1 7e 212 )D1 ] 1 9.50n 2 c 2 D2 e 221 1 3n 21 e 2 af2 e 23 sin(f 1 « 23 ), g1 2 g2 . g2 2 g2 From (A37)–(A39), the lowest order perturbation equations for g 1 and g 2 are as follows: ġ 2 5 2(2A3 2 n 1 e 2 af1 e 12 5 7n 1 c 1 D1 g1 2 g1 e2 n1 e1 n2 (A41) (A31) The terms involving the phase lags in the above equations can be simplified as follows: 1 7n 2 c 2 D2 (A40) and f 5 u 2 2 u 1 5 l 1 2 3l 2 1 2l 3 . 2 n2 (1 1 2r), g1 2 g2 n2 3 A3 5 2 e 1 e 2 f1 f2 , 2 g1 2 g2 where f is the Laplace resonance angle: 3e 2 n 21 a[ f1 e 12 « 1 n2 , g2 2 g2 (A34) ṅ 2 5 2n 2 c 2[1 2 (7e 22 1 7e 223 2 14e 221 )D2 ] 2 10.31n 1 c 1 D1 e 212 1 54.56n 3 c 3 D3 e 223 2 6n 22 e 1 f2 e 23 sin(f 1 « 23 ) (A35) 2 3n 22 e 3 af1 e 21 sin(f 2 « 21 ), ṅ 3 5 2n 3 c 3[1 2 (7e 23 2 14e 232 )D3 ] 2 1.80n 2 c 2 D2 e 223 (A36) 1 6n 23 e 2 f1 e 21 sin(f 2 « 21 ). Now, following Yoder and Peale (1981), we retain only the dissipative effects in Io, and we drop the distinction between e 1 and e 12 . Then, We note that the coefficient of sin f on the right-hand side of (A46) is positive. This means that the center of libration is near f 5 f, as is observed at the present epoch. It is useful to make some numerical estimates relevant to our simulations. For the value Q J 5 100 adopted in the numerical simulations, the magnitude of the tidal dissipation term in Eq. (A46) is p6 3 10210 n 22 , while the ‘‘strength of the resonance’’ represented by the coefficient of sin f has a magnitude p4 3 1026 n 22 (where we set n 2 /ug 1 2 g 2 u P 100, typical value that occurs when the Laplace resonance is encountered). In other words, the tidal torque is about four orders of magnitude weaker than the strength of the Laplace resonance, even with the ‘‘artificially enhanced’’ tidal torque that is necessary to do the numerical simulations in reasonable computer time. Therefore, we conclude that the adiabaticity of the real system near the Laplace resonance is preserved in our simulations. TIDAL EVOLUTION INTO LAPLACE RESONANCE Eliminating sin f from Eqs. (A44) and (A46), we obtain F f̈ 5 2bġ 1 1 (b 2 1)n 1 c 1 1 2 G 34.6b 2 44.6 D1 e 212 , b21 (A47) where b52 5 A 1 2 3A 2 1 2A 3 2A 2 2 A 1 3 1 (e 2 /ae 3 ) 1 6(1 1 (2ae 2 /3e 1 ))r 2 1 4(1 1 (ae 2 /e 1 ))r (A48) of the equations for n i and for (k i , h i ) [Eqs. (A3)–(A8)], this time retaining terms of the second order in eccentricity in the disturbing functions. By expanding the resulting expressions consistent to the next order in perturbation theory, we find new terms with arguments corresponding to the g 1 /g 2 5 j/( j 1 1) and g 1 /g 2 5 ( j 1 1)/j Laplacelike resonances. Omitting the details of the long and tedious algebra necessary to obtain the complete expansions, we reproduce below only the slow components corresponding to Laplacelike resonances: d 2 ṅ 1 5 13n 21 e 2 a[(2 f5 e 23 2 As f1B1 )e*1 sin(2u 1(0) 2 u 2(0) 2 Ã*1 ) 1 (4f4 e 23 2 As f2B1 )e*2 sin(2u 1(0) 2 u 2(0) 2 Ã*2 )], r 1 0.48 5 1.37 . r 1 0.37 (A53) Note that b is greater than 1 for all positive values of r and that its numerical value varies in the small range (1.4, 1.65) for r in the range p1/5 to p5. For the Laplace resonance, we recover Yoder and Peale’s (1981) results by setting r 5 1; then b 5 1.48, and we have f̈ 5 21.48ġ 1 1 0.48n 1 c 1(1 2 13.7D1 e 212 ). (A49) Furthermore, since f librates within the Laplace resonance, we have kf̈l 5 0. Then Eq. (A49) indicates that the equilibrium value of Io’s forced eccentricity in the Laplace resonance is given by e 12 5 (13.7D1 )1/2 where ġ 1 (and ġ 2 ) vanish. Away from the Laplace resonance, Eqs. (A44)–(A46) give the following first-order periodic components in the perturbations of u i : (2A2 2 A1 )n 22 d 1 u 1 5 B1 sin f, B1 5 , (g 2 2 g 1 )2 (A50) (2A3 2 A2 )n 22 . (g 2 2 g 1 )2 (A51) d 1 u 2 5 B2 sin f, B2 5 109 To obtain the slow frequency terms describing Laplacelike resonances, g 1 /g 2 P j/( j 1 1) or ( j 1 1)/j, we note that such terms can arise in higher-order perturbation theory through combinations of the lower-order solutions. Accordingly, we assume a solution of the form (k i , h i ) 5 (k*i 1 d 1 k i , h*i 1 d 1 h i ), (A52) where d 1 k i , d 1 h i are the first-order ‘‘fast’’ perturbations [Eqs. (A14)–(A16)], and k*i , h*i refer to the slow components corresponding to a Laplace-like resonance. We substitute the lowest-order solutions into the right-hand sides d 2 ṅ 3 5 26n 23 e 2[(4f3 e 21 1 As f1B2 )e*2 sin(2u 2(0) 2 u 1(0) 2 Ã*2 ) 1 (2 f5 e 21 1 As f2 B2 )e*3 sin(2u 2(0) 2 u 1(0) 2 Ã*3 )], (A54) d 2 ṅ 2 5 22 n 22 e 1 1 n 22 e 3 a d ṅ 2 d 2 ṅ 3 , 2 1 n 21 e 2 a 2 n 23 e 2 (A55) where u i(0) refers to the zeroth-order solution given in Eqs. (A17) and (A18). We see that there are two Laplace-like resonances, g 1 /g 2 P 1/2 and g 1 /g 2 P 2, that arise at this order in the perturbation analysis. Furthermore, there are two resonant arguments for each of these. Yoder and Peale (1981) had discussed the former resonances (i.e., the g 1 / g 2 P 1/2), but not the latter. Here we consider only one of the g 1 /g 2 P 2 resonances, namely, c 3 5 2u 2 2 u 1 2 Ã 3 , (A56) which is of interest for Ganymede. [In principle, we should first establish the conditions under which the two g 1 /g 2 P 2 Laplace-like resonances are decoupled and can be analyzed separately. Here we note only that to the order considered, there is no direct coupling between the c 3 and c 2 5 2u 2 2 u 1 2 Ã 2 resonances, i.e., we find no terms with arguments that are linear combinations of c 2 and c 3 , and use this as justification for analyzing the c 3 resonance in isolation.] From (A53)–(A55) we obtain the following perturbation equations for the c 3 components in g 1 and g 2 , ġ*1 5 16n 22 e 3 a(2 f5 e 21 1 As f2 B2 )e*3 sin c *3 1 n 1 c 1(1 2 34.6D1 e 212 ), (A57) ġ*2 5 2(12n 33 e 2 1 3n 22 e 3 a)(2 f5 e 21 1 As f2 B2 )e*3 sin c *3 (A58) 1 10.3n 1 c 1 D1 e 212 , and the pendulum-like equation for c *3 , 110 SHOWMAN AND MALHOTRA c̈ *3 Q 2(24n 23 e 2 1 12n 22 e 3 a)(2 f5 e 21 1 As f2 B2 )e*3 sin c *3 2 n 1 c 1(1 2 55.2D1 e 212 ), (A59) where we have neglected Ã̈ *3. Equation (A59) shows that center of libration of c *3 is near 0. Comparison of Eq. (A59) with Eq. (A46) shows that the c 3 resonance is weaker than the Laplace resonance by a factor of pe*3. This means that for the value Q J 5 100 adopted in our numerical simulations, the tidal dissipation term in Eq. (A59) is about two orders of magnitude weaker than the strength of the c 3* resonance [see discussion following Eq. (A46)]. Therefore, we conclude that the adiabaticity of the evolution of the real system near this resonance is preserved in our simulations. Eliminating e*3 sin c 3* from (A57) and (A59), we have c̈ *3 5 22.51ġ 1 1 1.51n 1 c 1(1 2 21.0D1 e 212 ). (A60) When the satellites are locked in this resonance, we have kc̈ 3*l5 0. Equation (A60) then indicates that within the c 3 resonance Io’s forced eccentricity e 12 reaches an equilibrium value at e 12 5 (21D1 )21/2, where ġ 1 vanishes; however, this is not a true equilibrium of the system because Ganymede’s forced eccentricity continues to grow. The forced eccentricity of Ganymede within the c 3 resonance is obtained by considering the equations for k*3 , h*3 , including terms up to the second order in eccentricity in the disturbing function, using the first-order solutions [(A14)–(A16), (A50), (A51)] to obtain the (2u 2 2 u 1 ) component in the perturbation equations. We find k̇*3 5 2e 2 n 3( f5 e 21 1 As f2 B2 ) sin(2u 2(0) 2 u 1(0) ), ḣ*3 5 1e 2 n 3( f5 e 21 1 As f2 B2 ) cos(2u (0) 2 2u (A61) (0) 1 ), which yield the following approximate solution, k*3 5 e*32 cos(2u 2(0) 2 u 1(0) ), h*3 5 e*32 sin(2u 2(0) 2 u 1(0) ), (A62) where e*32 5 n 3 e 2( f5 e21 1 As f2 B2 ) . 2g 2 2 g 1 (A63) Substituting the expressions for e 21 [Eq. (A22)] and B2 [Eq. (A51)] in (A63), setting 2g 2 5 g 1 everywhere except in the combination 2g 2 2 g 1 , and using numerical values of e i and the coefficients fi , we have e*32 Q 1.26 3 1029 F S DS S DG n2 2g 2 2 g 1 1 2 1.67 3 1024 n2 g1 n2 g1 2 g2 D (A64) 2 . We note that as ug 1 u decreases, the first two factors contribute to increasing Ganymede’s forced eccentricity, while the factor in square braces diminishes it. Thus the general behavior is that e*32 at first increases, reaches a maximum, then decreases and vanishes (and the c 3 resonance is disrupted) as ug 1 u decreases. [Numerical experiments revealed similar behavior of the forced eccentricity of Ganymede as well as Europa in other Laplace-like resonances; see Malhotra (1991).] The factor in square braces vanishes when g 1 /n 2 Q 0.013. The forced eccentricity of Io corresponding to this value is P0.003; however, long before g 1 and e 12 reach these values, we expect the above solution to break down. The approximate ‘‘forced’’ solution for k*3 , h*3 given in Eqs. (A62) and (A63) is valid provided 2g 2 2 g 1 varies slowly and its magnitude is not too small compared with nuee 21 e 3 u 1/2, the quantity that describes the approximate strength of the c 3 resonance [see Eq. (A59)]. (More precisely, use of a Hamiltonian formalism yields the critical value, e*3,crit , for the c 3 resonance which provides an equivalent estimate of the regime of validity of the approximate forced solution. We find e 3,crit P 1023.) Thus, we conclude that the above solution breaks down by the time Ganymede’s eccentricity is excited to values p1023 after the satellites are captured in the c 3 resonance. We also conclude that capture into the c 3 resonance has high probability if Ganymede’s initial eccentricity is smaller than 1023. We do not pursue this analysis further, for the degree of difficulty increases greatly at this stage. Many uncontrolled approximations are necessary to push the analysis further, and it has proven difficult to control the proliferation of terms that occurs in a higher-order analysis. The above analysis does show that (1) low-order Laplace-like resonances exist that can excite Ganymede’s eccentricity, and (2) the g 1 /g 2 P 2 resonance is sufficiently strong for resonance capture to occur and for the c 3 libration to be maintained. The analysis increases our confidence in the validity of the numerical results. ACKNOWLEDGMENTS We thank Paul Schenk for useful discussions. A.S. gratefully acknowledges the Lunar and Planetary Institute Visiting Student and National Science Foundation graduate fellow support. This research was done while one of the authors (R.M.) was a Staff Scientist at the Lunar and Planetary Institute, which is operated by the Universities Space Research Association under contract NASW-4574 with the National Aeronautics and Space Administration. TIDAL EVOLUTION INTO LAPLACE RESONANCE REFERENCES BANFIELD, D., AND N. MURRAY 1992. A dynamical history of the inner neptunian satellites. Icarus 99, 390–401. BERCKHEMER, H., W. KAMPFMANN, E. AULBACH, AND H. SCHMELING 1982. Shear modulus and Q of forsterite and dunite near partial melting from forced-oscillation experiments. Phys. Earth Planet. Inter. 29, 30–41. CHAPMAN, C. R., AND W. B. MCKINNON 1986. Cratering of planetary satellites. In Satellites (J. A. Burns and M. S. Matthews, Eds.), pp. 492–580. Univ. of Arizona Press, Tucson. GAVRILOV, S. V., AND V. N. ZHARKOV 1977. Love numbers of the giant planets. Icarus 32, 443–449. GOLDREICH, P., 375–389. AND S. SOTER 1966. Q in the Solar System. Icarus 5, GREENBERG, R. 1982. Orbital evolution of the Galilean satellites. In Satellites of Jupiter (D. Morrison, Ed.), pp. 65–92. Univ. of Arizona Press, Tucson. GREENBERG, R. 1987. Galilean satellites: Evolutionary paths in deep resonance. Icarus 70, 334–347. HENRARD, J. 1983. Orbital evolution of the Galilean satellites: Capture into resonance. Icarus 53, 55–77. IOANNOU, P. J., AND R. S. LINDZEN 1993. Gravitational tides in the outer planets. II. Interior calculations and estimation of the tidal dissipation factor. Astrophys. J. 406, 266–278. MALHOTRA, R. 1991. Tidal origin of the Laplace resonance and the resurfacing of Ganymede. Icarus 94, 399–412. MALIN, M. C., AND D. C. PIERI 1986. Europa. In Satellites (J. A. Burns and M. S. Matthews, Eds.), pp. 689–717. Univ. of Arizona Press, Tucson. MCKINNON, W. B., AND E. M. PARMENTIER 1986. Ganymede and Callisto. In Satellites (J. A. Burns and M. S. Matthews, Eds.), pp. 718–763. Univ. of Arizona Press, Tucson. MCKINNON, W. B., AND P. M. SCHENK 1995. Estimates of comet fragment masses from impact crater chains on Callisto and Ganymede. Geophys. Res. Lett. 22, 1829–1832. 111 OJAKANGAS, G. W., AND D. J. STEVENSON 1986. Episodic volcanism of tidally heated satellites with application to Io. Icarus 66, 341–358. PEALE, S. J., AND P. M. CASSEN 1978. Contribution of tidal dissipation to lunar thermal history. Icarus 36, 245–269. PEALE, S. J., P. CASSEN, AND R. T. REYNOLDS 1979. Melting of Io by tidal dissipation. Science 203, 892–894. SHOEMAKER, E. M., B. K. LUCCHITTA, D. E. WILHELMS, J. B. PLESCIA, AND S. W. SQUYRES 1982. The geology of Ganymede. In Satellites of Jupiter (D. Morrison, Ed.), pp. 435–520. Univ. of Arizona Press, Tucson. SHOWMAN, A. P., AND D. J. STEVENSON 1996. Resurfacing of Ganymede and other icy satellites by passage through orbital resonance. Icarus, submitted. SHOWMAN, A. P., D. J. STEVENSON, AND R. MALHOTRA 1996. Coupled orbital and thermal evolution of Ganymede. Icarus, submitted. SINCLAIR, A. T. 1975. The orbital resonance amongst the Galilean satellites of Jupiter. Mon. Not. R. Astron. Soc. 171, 59–72. SMITH, B. A., L. A. SODERBLOM, T. V. JOHNSON, A. P. INGERSOLL, S. A. COLLINS, E. M. SHOEMAKER, G. E. HUNT, H. MASURSKY, M. H. CARR, M. E. DAVIES, A. F. COOK II, J. BOYCE, G. E. DANIELSON, T. OWEN, C. SAGAN, R. F. BEEBE, J. VEVERKA, R. G. STROM, J. F. MCCAULEY, D. MORRISON, G. A. BRIGGS, AND V. E. SUOMI 1979. The Jupiter system through the eyes of Voyager 1. Science 204, 951–972. STEVENSON, D. J. 1983. Anomalous bulk viscosity of two-phase fluids and implications for planetary interiors. J. Geophys. Res. 88, 2445–2455. TITTEMORE, W. C. 1990. Chaotic motion of Europa and Ganymede and the Ganymede–Callisto dichotomy. Science 250, 263–267. TITTEMORE, W. C., AND J. WISDOM 1989. Tidal evolution of the uranian satellites: II. An explanation of the anomalously high orbital inclination of Miranda. Icarus 78, 63–89. VEEDER, J. G., D. L. MATSON, T. V. JOHNSON, D. L. BLANEY, AND J. D. GOGUEN 1994. Io’s heat flow from infrared radiometry: 1983–1993. J. Geophys. Res. 99, 17,095–17,162. YODER, C. F. 1979. How tidal heating in Io drives the Galilean orbital resonance locks. Nature 279, 767–770. YODER, C. F., AND S. J. PEALE 1981. The tides of Io. Icarus 47, 1–35.