Spin-orbital tidal dynamics and tidal heating in the TRAPPIST-1 multi-planet system
aa r X i v : . [ a s t r o - ph . E P ] A p r Spin-orbital tidal dynamics and tidal heatingin the TRAPPIST-1 multi-planet system
Valeri V. Makarov
US Naval Observatory, Washington DC 20392e-mail: valeri.makarov @ navy.mil
Ciprian T. Berghea
US Naval Observatory, Washington DC 20392e-mail: ciprian.berghea @ navy.mil
Michael Efroimsky
US Naval Observatory, Washington DC 20392e-mail: michael.efroimsky @ navy.mil
ABSTRACT
We perform numerical simulations of the TRAPPIST-1 system of seven exo-planets orbiting a nearby M dwarf, starting with a previously suggested stableconfiguration. The long-term stability of this configuration is confirmed, but themotion of planets is found to be chaotic. The eccentricity values are found to varywithin finite ranges. The rates of tidal dissipation and tidal evolution of orbits areestimated, assuming an Earth-like rheology for the planets. We find that underthis assumption, the planets b , d , e were captured in the 3:2 or higher spin-orbitresonances during the initial spin-down, but slipped further down into the 1:1 res-onance. Depending on its rheology, the innermost planet b may be captured in astable pseudosynchronous rotation. Nonsynchronous rotation ensures higher levelsof tidal dissipation and internal heating. The positive feedback between the viscos-ity and the dissipation rate — and the ensuing runaway heating — are terminatedby a few self-regulation processes. When the temperature is high and the viscosityis low enough, the planet spontaneously leaves the 3:2 resonance. Further heating isstopped either by passing the peak dissipation or by the emergence of partial meltin the mantle. In the post-solidus state, the tidal dissipation is limited to the levelssupported by the heat transfer efficiency. The tides on the host star are unlikely tohave had a significant dynamical impact. The tides on the synchronized inner plan-ets tend to reduce these planets’ orbital eccentricity, possibly contributing therebyto the system’s stability. Subject headings: planet-star interactions — planets and satellites: dynamical evo-lution and stability — celestial mechanics — planets and satellites: individual(TRAPPIST-1) 2 –
1. Introduction
Main-sequence M dwarfs are a numerically dominant stellar species in the Galaxy, andthe possibility of finding a terrestrial planet in their narrow habitable zones provides the bestchance of locating a world suitable for biological life in the relative vicinity of the Sun. Super-Earth exoplanets in the habitable zones of nearby M dwarfs are probably abundant, withan estimated rate of η Earth = 0 . +0 . − . per host star (Bonfils et al. 2013). Super-Earths aremore massive than the Earth, but smaller than 10 Earth masses. The recent discovery ofthe TRAPPIST-1 system of seven planets (Gillon et al. 2017) transiting a nearby, otherwiseunremarkable M dwarf (Reiners & Basri 2010; Filippazzo et al. 2015; Gillon et al. 2016) opensa new realm of extrasolar worlds – tightly packed systems of Earth-sized planets. The estimatedmasses of planets in this system range from 0 . . M Earth . With the currently outermostplanet orbiting the star in less than 19 days, we are challenged to explain how such compactsystems emerge and why they remain stable on the astronomical time scale.A planet’s rotation is one of the determinants for the presence and properties of its atmo-sphere, and for the distribution of the stellar irradiation on the surface. The climate patternmay be radically different on planets in the habitable zone, depending on the spin-orbit reso-nance they reside in (Wang et al. 2014). Physically justified, frequency-dependent tidal mod-els, such as those based on the Maxwell or combined Maxwell–Andrade rheologies (Efroimsky2012a,b), predict high probabilities of capture into higher-than-synchronous spin-orbit reso-nances for close-in exoplanets with finite eccentricities (Makarov 2012). This result holds for awider range of realistic models, such as the Newtonian creep model developed by Ferraz-Mello(2015). The orbital eccentricity and average viscosity are the two crucial parameters definingthese probabilities.The issue of dynamical stability and long-term survival of a tightly packed close-in systemsuch as TRAPPIST-1 has already been discussed in the discovery paper by Gillon et al. (2017).A system with the observed set of parameters breaks up within a relatively short time innumerical N -body integrations. Fine-tuning of initial orbital parameters is required to find along-term stable configuration within the observational uncertainties. This indicates that themulti-dimensional parameter space is shredded, with small islands of stable solutions separatedby areas of instability. How can a real physical system of high complexity converge on a smallisland of stability among an infinite number of possible unstable configurations? The stabilityislands should be attractors, i.e., there should exist a regulatory force driving the systemto one of such islands and damping small perturbations. Gillon et al. (2017) suggested thatconsiderable tidal effects can provide the required stabilization. However, tidal forces act muchtoo slowly even for such close planets, while the break-up of an unstable system happens wellwithin 1 Myr. Even though tidal effects are important on the lifetime scale, we still have toexplain how the system originally found a stable configuration. Tamayo et al. (2017) proposedthat a convergent migration in a protoplanetary disk can be such a regulatory mechanism.Subject to dynamic friction in the disk, young planets migrate generally inward and towardsmaller eccentricities. Still, something would halt this process when a planet reaches a stabilityniche, even before the disk dispersal. This may be possible when the stability niche is also 3 –associated with a certain mean-motion resonance (MMR), or better still, a system of MMRswith other planets, because an MMR is known to be an attractor. The surviving protoplanetsend up in intricate patterns of near-commensurability (Terquem & Papaloizou 2007), a classof N -body systems radically different from the Solar System. Luger et al. (2017) find that theTRAPPIST-1 system represents a chain of MMRs, where each triplet of adjacent planets islocked in a three-body resonance.In this paper, we begin with testing a long-term stable configuration suggested by Quarles et al.(2017). We also test the system for chaos, using the siblings method in numerical integrations.The resulting orbital eccentricities for the seven planets, termed for brevity as TR1-1, TR1-2,. . . , TR1-7 (instead of the traditional notation TRAPPIST-1 b , c , . . . , h ), are used to estimatethe likely spin-orbit end-states of the planets, assuming an Earth-like rheology described by theMaxwell model. We show that not all of the planets are expected to fall in the 1:1 resonance(synchronous rotation) if they remain cold and inviscid during the initial spin-down. The in-nermost planets, in particular, are more likely to be initially locked in higher resonances (2:1or 3:2) despite their low mean eccentricity. This gives a significant boost to the estimated rateof energy dissipation. The rate of tidal heating is so high for TR1-1 that one may wonder if theplanet should be completely molten. We further describe the mechanism of tidal heating self-regulation, which could prevent a planet with a rocky mantle from going past solidus. Besidesdefining the tidal heating rates, the rotation states of tidally deformed bodies influence theirorbital evolution. It should not be taken for granted that the tidal forces always circularizeand shrink orbits. We show that the tides in a fast-rotating host star may have the oppositeeffect on the outer planets’ orbits.Our research is based on analytical and semianalytical techniques. Applicable only touniform or two-layered bodies, the analytical and semianalytical techniques do not lead to aself-consistent detailed modeling of internal processes, such as the internal heat transfer or thesurface irradiation. In contradistinction, numerical simulations of tidal dissipation for specifictemplate cases provides certain advantages. They allow for a more detailed investigation intothe effects of a multilayered structure, and make it easier to explore the influence of partic-ular rheologies, such as Maxwell (Bou´e et al. 2016; Walterov´a & Bˇehounkov´a 2017), Burgers(Henning et al. 2009), and Kelvin–Voigt (Zlenko 2015; Frouard et al. 2016). However, theseadvantages of numerics are all but gone in view of the great uncertainties we meet when analyz-ing the tidal scenarios for exoplanets. We are limited to an approximate, order-of-magnitudeestimation based on a wide range of possible values for the most significant parameters, not tomention the remaining ambiguity of physical modeling of the layers. For this reason, useful in-sight can be gained through a low-cost analytical or semi-analytical approach; and hypothesesinstilled in this approach can be tested against the observational data. 4 –
2. The orbits of the TR1 planets
For our analysis of the orbital stability of the TR1 system, we used observational datasets from Gillon et al. (2017) and the configuration reported by Quarles et al. (2017, Table 3)as long-term stable. The integrations, which used the observed orbital parameters as initialconditions, invariably broke apart after 10 –10 years. However, the configuration presentedby Quarles et al. (2017) remained intact for the duration of integration (which was 10 Myrin our case). Among the initial configurations tested by those authors, 28.2% proved to bestable. We confirm that relatively small changes in the initial values of the orbital elementscan result in a stable system, suggesting the small size of the stability islands in the sea ofunstable configurations. Time (yr) −4 −3 −2 E cc e n t r i c i t y Fig. 1.— Numerically integrated behavior of planet h (TR1-7) eccentricity.We find that the eccentricity of each planet rapidly varies in an apparently stochastic way,but remains within a certain range for the entire duration. Figure 1 shows the result forthe outermost planet TR1-7. Other derived characteristics are given in Table 1, including themean, median, and robust standard deviation. The latter values, for example, are found in therelatively narrow range [0 . , . − . The absolute distance in astronomical units betweenthe siblings was then computed as a function of time. Chaotic motion manifests itself as anexponential divergence between the trajectories. We find such exponential divergence in thetrajectories of all seven planets. Figure 2 demonstrates this result for planet TR1-7. The rateof divergence allows us to estimate the Lyapunov time, which is approximately 1.2 yr for thisplanet. The conclusion is that even a tiny change or uncertainty in the orbital parametersTable 1. Numerically simulated statistics for orbital eccentricity of TRAPPIST-1 planets.Planet Range Mean Median StD Robust Mean Robust StD1 . . . . b . . . . . . c . . . . . . d . . . . . . e . . . . . . f . . . . . . g . . . . . . h . . . . . . Time (yr) −9 −7 −5 −3 −1 D i s a n c e b e w ee n s i b li n g s ( A U ) Lyapunov ime 1.2 yr
Fig. 2.— Simulated distance in semimajor axis a between two sibling trajectories for planetTR1-7. 7 –makes the long-term prediction of the future state impossible.In our previous simulations of exoplanet systems, we found clear signs of chaos for GJ581 (Makarov et al. 2012), but not for GJ 667C (Makarov & Berghea 2014). The former is acomplex system with a few planets interacting, whereas the latter is assumed to include onlytwo planets locked in the 2:1 MMR. Chaotic behavior may be common for systems with morethan two planets, even when some of them are locked in mutual resonances.
3. The Maxwell model and the tidal quality function k /Q In the Solar System, Mercury is the closest analogue to the TR1 planets, in the sense of or-bital separation and mass. Mercury is believed to have a massive molten core, whose existenceis revealed by the large amplitude of forced libration detected via ground-based radar measure-ments (Margot et al. 2007). The core should be nearly decoupled from the solid mantle (i.e.,it should have a low coefficient of friction) to enable the increased libration. Unlike the TR1planets, however, Mercury’s orbital eccentricity is quite high, which explains its permanent3:2 spin-orbit resonance. Early attempts to explain this state by using simplistic tidal models(such as the constant- Q or constant-time-lag models) met considerable difficulties, renderinglow values for the probability of any equilibrium state apart from the synchronous rotation,the 1:1 spin-orbit resonance (Correia & Laskar 2004). Noyelles et al. (2013) demonstrated thatthese difficulties are resolved within a physics-based – and, therefore, more realistic – descrip-tion of the rheological response. Based on a combined Maxwell–Andrade rheology introducedon physical grounds by Efroimsky (2012a,b), the tidal response demonstrated a pronouncedfrequency dependence. Armed with this model, and employing physically reasonable ranges ofthe parameters’ values, Noyelles et al. (2013) demonstrated that the 3:2 resonance is indeedthe most likely outcome of the tidal evolution of a homogeneous rotator, and that the captureis a very quick a process (less than 20 Myr). On the other hand, the presence of a massivedecoupled core made the higher resonances, such as 5:2 or 2:1, more probable end states. Theseresults led the authors to the conclusion that, in all likelihood, Mercury was captured in the3:2 spin-orbit state shortly after its accretion, i.e., well before it was heated up by radioactivity.The authors also acknowledged the possibility that Mercury might have been trapped in oneof such high resonances, before a massive impact during a period of relatively low eccentricitydrove it out of that state.The Maxwell model has long been in use in the studies of exoplanets and the solar-systembodies (e.g., Storch & Lai 2014; Correia et al. 2014). The parameters entering this model —rigidity and viscosity — cannot be directly constrained by observation. It is also difficult toestimate them theoretically, because these parameters are sensitive to the presence of partialmelt, to the pressure, and even more so to the temperature. The high temperatures andpressures required to emulate the conditions of a deep mantle are hard to achieve in laboratorymeasurements, and tidal manifestations remain the principal method for the indirect estimationof these parameters. Furthermore, the tidal response of terrestrial planets, satellites, and small 8 –bodies turns out to be less sensitive to the rigidity and much more sensitive to the viscosity(Efroimsky 2015). Figure 3 shows the dimensionless tidal quality function, k /Q , plottedagainst the excitation frequency for two values of the Maxwell time of the TR1-1 planet. Forthe shorter Maxwell time (implying a lower viscosity value), the dependence is almost linear,and the classic constant-time-lag model turns out to be a good approximation. For somewhatcolder mantles of higher viscosity, the function is remarkably nonlinear, and conclusions canbe drawn only from numerical simulations. τ M d τ M (cid:0) χ q / n qua li t y f a c t o r K ( χ ) = k / Q Fig. 3.— Tidal quality ratio K = k /Q for the planet TR1-1, as a function of the tidalexcitation frequency, plotted for two values of the Maxwell time, τ M = 0 .
02 and 0 .
4. Self-regulation of tidal heating
As the partners tidally interact with one another over billions of years, tides become animportant factor of the dynamical history of the system. The rate of tidal dissipation dependson a number of orbital and planetary parameters, including the eccentricity e , the masses of thestar and the planet ( M and M ), the semimajor axis a , the planet’s radius R , and its Maxwelltime τ M . The dependences of the tidal evolution rate and of the capture probabilities on theseparameters are quite strong and generally nonlinear, making the physical interpretation ofobservations intrinsically uncertain. The approach we are using here is to assume some of theparameters to be close to the Earth’s benchmark values when available, and to consider a widerange of values otherwise. The Maxwell time is the ratio of the mean shear viscosity of the mantle to its mean shear rigidity modulus, τ M = η/µ . dE/dt τ M W dTR1-1 3 . × . . × . . × . . × . . × . . × . . × . - Maxwell time d
Log ( ⅆ E / ⅆ t ) W Fig. 4.— Estimated rate of tidal energy dissipation for planet TR1-1 as a function of theMaxwell time τ M , for three possible spin-orbit resonances, 1:1, 3:2, and 2:1. 10 –Figure 4 displays the Maxwell-time dependence of the energy dissipation rate (i.e., of thepower exerted by the tidal friction) in the planet TR1-1, for a fixed orbital period P orb =1 .
51 d. We used the formula for the tidal dissipation rate derived by Efroimsky & Makarov(2014), which is a generalization of a result obtained earlier by Peale & Cassen (1978). WhilePeale & Cassen (1978) were dealing exclusively with a synchronous case, Efroimsky & Makarov(2014) extended their theory to an arbitrary spin state, resonant or nonresonant. Due to thisplanet’s considerable size and its proximity to the star, the tidal dissipation in this planet isvery intense. It surpasses by many orders of magnitude the tidal dissipation rates typical forthe solar-system bodies. The peak rate for the 1:1 spin-orbit resonance, which is the likeliestspin state (see Section 5), reaches 4 × W. This results in an equilibrium surface flux of6 . × W m − . By comparison, the thermal flux on Earth, caused by the tides from theMoon and the Sun, amounts to only ≃ .
001 W m − (Henning & Hurford 2014). As can beseen from the plot, the Maxwell time corresponding to the maximal power is 0.048 d. Themaximal dissipation rates and the corresponding values of τ M are given in Table 2 for the sevenplanets, under the synchronous rotation assumption.
5. The likely spin–orbit states of the TR1 planets
The rotation history of a tidally perturbed elongated body is governed by the torque ~ T ( TRI ) due to the body’s permanent triaxiality and the torque ~ T ( TIDE ) due to tides in the body (Danby1962; Goldreich & Peale 1966, 1968). At small obliquities (i.e., when the perturber’s orbitalplane is close to the equatorial plane of the tidally perturbed body), only the polar componentsof the torques matter, and the equation of motion acquires the form (cid:5) (cid:5) θ = T (TRI) polar + T (TIDE) polar C = T (TRI) polar + T (TIDE) polar ξ M R , (1)Here M , R and C = ξ M R are the mass, radius and the maximal moment of inertia ofthe body, while ξ is a dimensionless factor (equal to 2 / θ of the body is conventionally reckoned from the line of apsides to the largestelongation axis.The permanent-triaxiality-generated torque ~ T ( TRI ) is oscillating, while the tidal torque ~ T ( TIDE ) contains both a secular and an oscillating parts (Efroimsky 2012a). The long-termevolution and the likely end states of the rotator are defined mainly by the interplay of ~ T ( TRI ) and the secular part of ~ T ( TIDE ) .It is commonly accepted that a planet rapidly rotates at the time of formation, whetherin the prograde or retrograde sense. The gradual dissipation of the rotational energy by thetidal friction causes angular deceleration. For close-in planets, the energy damping rate is While the oscillating part of the tidal torque changes the fate of each particular trajectory, the overallstatistics of tidal capture in spin-orbit resonances is defined overwhelmingly by the secular part (Makarov et al.2012).
11 –intensive enough to make the spin-down time short (e.g., 10 −
20 Myr for Mercury). Still, itis not sufficiently intensive to cause a significant tidal heating, because the overall supply ofthe rotational energy is limited. Simple calculations show, for example, that the lengtheningof the solar day on the Earth results in a heating much smaller than the surface heating bythe solar irradiation or than the internal heating by the radioactive decay. So the planetquickly reaches a spin-orbit equilibrium state before having undergone a significant internalrestructuring. The subsequent history of the planet after the capture much depends on whichequilibrium spin-orbit state it has ended up in. p r obab ili t y o f c ap t u r e Fig. 5.— Probabilities of capture of planet TR1-1 into spin-orbit resonances 3:2, 2:1, and 5:2(left to right)It has often been assumed in the literature, incorrectly, that all close-in exoplanets aretidally synchronized. Facing the stars with one side for billions of years, such planets shouldbe quite limited in their ability to support a highly developed biological life, or even a stableatmosphere rich in volatiles. This paradigm has to be reconsidered, however. With the ex-ample of two nearby systems, GJ 581 and 667, harboring potentially habitable super-Earths,Makarov et al. (2012); Makarov & Berghea (2014) posited that the eccentricities can be largeenough for these planets to be captured in a 3:2 or even higher order spin-orbit resonance.This leads to a much more uniform distribution of stellar irradiance on the surface, but also toincreased tidal heating. Our task now is to estimate how likely it may be for the TR1 planetsto be captured in such Mercury-like resonances.An early formula for the probability of capture in a spin-orbit resonance was proposedby Goldreich & Peale (1968) for a CPL (constant phase-lag) model, i.e., for a model whereinthe intensity of tidal interaction was frequency-independent. The formula was generalizedby Makarov (2012) for a general frequency-dependent tidal model. Here we used the latterdevelopment to compute capture probabilities for a range of values of τ M , the most uncertain 12 –parameter. The rheology was presumed to be Maxwell — which should be adequate for close-inexoplanets of terrestrial composition. Figure 5 shows the results for the innermost planet TR1-1, which is subject to the strongesttidal interaction, versus eccentricity, with a fixed Maxwell time of 1 d. Depending on the valueof eccentricity at the time of capture, the planet might have been trapped in any of thehigher-than-synchronous resonances, 3:2, or 2:1, or higher. It is, however, doubtful that for adynamically stable configuration the eccentricity could be much larger than its present value(0.011). A tightly packed multiplanet system with initial high eccentricities would be unlikelyto survive in its present-day configuration. Therefore, using the fixed values emerging fromextended orbital integrations, e.g., the robust means shown in Table 1, is justified. Table3 renders the thus estimated probabilities of capture into the 3:2 and 2:1 spin-orbit states.We determine that planets TR1-2, TR1-5, TR1-6, and TR1-7 should almost certainly besynchronized. On the other hand, planets TR1-1, TR1-3, and TR1-4 are likely to have beeninitially captured in the 3:2 resonance, with a finite probability of an even higher spin state.
71 5 10 50 1000.00.20.40.60.81.0 Maxwell time d p r obab ili t y o f c ap t u r e i(cid:1)(cid:2)(cid:3) : Fig. 6.— Probabilities of capture of the TRAPPIST-1 planets in the 3:2 spin-orbit resonanceversus the Maxwell time.The largest uncertainty in this analysis comes from the assumed Maxwell time. Thisparameter, poorly known even for the solar-system bodies, is a strong function of the mantle’stemperature and pressure through its proportionality to viscosity (see Section 6 below). Theinner planets of the system receive vast amounts of tidal heating and should be hot, lowering For physical reasons, the Maxwell model should be substituted with a combined Maxwell-Andrade model,when the forcing frequencies are higher than a certain threshold (Efroimsky 2012a). Exponentially sensitiveto the temperature, the value of this threshold is of the order of yr − for the Earth, but is much higher forwarmer planets (Karato & Spetzler 1990, Eq. 17). So for close-in planets (presumed sufficiently warm) theAndrade mechanism may be neglected and the tidal reaction of the mantle may be presumed Maxwell.
13 –the Maxwell time. Figure 6 shows how the probability of capture in the 3:2 spin-orbit resonancedepends on the effective Maxwell time. Cold planets have higher viscosity and are more readilycaptured into this resonance. Planets TR1-1, 3, and 4 are certainly captured if their Maxwelltime is longer than ∼
40 d. The outermost planets 6 and 7 can be captured, in principle, ifthey were cold (like the Earth today) at the time of their initial spin-down. Planets 2 and 5cannot be in this resonance because of their low eccentricities.Apart from the strong dependence on the Maxwell time τ M and eccentricity e , theprobabilities of capture also depend on the triaxiality σ = ( B − A ) /C , which is unknown. Thisdependence, however, is relatively weak. Near-oblate bodies (i.e., those with a vanishinglysmall ( B − A ) /C ) are always captured into a resonance if in its vicinity the secular part ofthe tidal torque is positive (accelerating in the prograde sense). A rule to keep in mind isthat less triaxial planets are more readily captured into higher-than-synchronous resonances,with all other parameters being equal. Based on the data for the Solar System, larger planetsare less triaxial than the smaller planets or moons. The triaxiality value we assumed here,( B − A ) /C = 1 . × − , is close to that of the Earth (Lambeck 1980) and therefore is areasonable educated guess.High-energy impact is a stochastic factor that may also define the spin-orbit states ofexoplanets. Within some scenarios, the current 3:2 spin-orbit state of Mercury could havebeen the outcome of such impacts in the past (Correia & Laskar 2012; Noyelles et al. 2013).Speculatively, early collisions at the stage of planetary system formation could not have a long-lasting effect because the eccentricities have settled to finite values, and the tidal mechanismshad sufficient time to change the spins of the TR1 planets. If a late bombardment episodewere to happen in this settled state, fortuitous tangential hits could have temporarily drivensynchronously rotating planets into higher spin-orbit resonances, causing periods of increasedinternal heating.The inner TR1 planets may also be subject to considerable retrograde torques originatingfrom their hypothetical atmospheres. Owing to the friction coupling between the atmosphereand the planet’s surface (Dobrovolskis & Ingersoll 1980), the secular thermal torque on theatmosphere is transferred to the planet, counterbalancing the action of the lagging solid-bodybulge. This is believed to explain the retrograde rotation of Venus with its massive and thickatmosphere. Leconte et al. (2015) discussed that close exoplanets orbiting low-mass stars mayalso be dominated by thermal atmospheric tides (even if their atmospheres are thin) and maymaintain a stable asynchronous rotation. We defer analysis of this more complex configurationuntil future times when more knowledge is acquired about the atmospheres of close-in super-Earths. The existence of such atmospheres cannot currently be taken for granted even for thehabitable zones (Heng & Kopparla 2013). 14 –
6. What comes first, solidus or peak heating?
Investigating the dependence of the tidal heating rate on the average viscosity and tem-perature of Earth-like planets, Bˇehounkov´a et al. (2011) found that a positive feedback shouldlead to a runaway effect, i.e., to a situation where a higher mantle temperature causes evenmore heating. This can impact the habitability of close-in planets in the habitable zones oflow-mass stars. For synchronously rotating planets, the possibility of runaway heating dependson the orbital eccentricity value and on the presence of plate tectonics, which strongly enhancesthe heat transfer in the mantle. The feedback is seen in our Figure 4 as the descending branchof the function on the right-hand side. Under synchronous rotation, the excitation frequencyis (nearly) constant in time, but the effective viscosity and Maxwell time change with therise of the temperature. An initially cold inviscid planet with a large Maxwell time beginsto warm up if the heat transfer does not keep up with the energy input. Then, the viscositydecreases, and even more energy is dissipated. However, Figure 4 also shows that there is aself-regulation mechanism embedded in the Maxwell model. When the peak dissipation haspassed, the steady warm-up and the shortening of the Maxwell time lead to a reduction of thedamping rate. Using Figure 3, this can be visualized as the peak of the quality function shiftingto the right as the viscosity drops and the value of the quality function is rapidly decreasingat a fixed frequency.A more radical and universal stop-mechanism for the runaway heating is provided bythe phase transition of solid minerals when they begin to melt (Henning & Hurford 2014;Shoji & Kurita 2014). Both the viscosity and the rigidity modulus precipitously drop, and therate of dissipation is reduced by orders of magnitude. Even a limited fraction of partial meltis sufficient to bring about these drastic changes. This self-regulation mechanism is likely tohave radically slowed down the initial expansion of the Moon’s orbit by keeping the dissipationrate low despite the small orbital separation (Zahnle et al. 2015). We do not expect eventhe closest exoplanets to be completely molten to the surface (magma ocean worlds) — atleast, not by the tidal forces. However, since these processes can take place on potentiallyhabitable planets orbiting M-type stars, the important question is whether the point of solidusis reached sooner than the peak dissipation, i.e., which of the self-regulation mechanisms haltsthe runaway internal heating.We analyzed three different rheological models proposed in the literature in applica-tion to the TR1 planets. These models assume a steeper, Arrhenius-type dependence ofthe mineral viscosity on temperature than the Frank-Kamenetskii approximation used byBˇehounkov´a et al. (2011, ;) (the implications of this choice are discussed in, e.g., Stein & Hansen2013). We find that the model proposed by Doremus (2002), which is limited to the range oftemperatures [1000 , ◦ C , provides unrealistically low values of viscosity at higher tem-peratures that are expected in the deep layers of tidally heated exoplanets. The rheologicalmodel described by Stamenkovi´c et al. (2012) gives more reasonable results when applied tothe inner TR1 planets. For example, for the conditions expected at the bottom of the Earth’smantle ( T = 4000 K and P = 136 GPa), the estimated viscosity is η = 1 . × Pa · s. Thecorresponding Maxwell time is 264 d, assuming an effective rigidity of µ = 8 × Pa. We 15 –consider these values to be the benchmark solidus viscosity. In fact, the presolidus rigidity isalso a function of temperature and pressure, but this dependence is much weaker than that forviscosity.A surprising result emerging from this model is the long Maxwell time for the presolidusdeep mantle of the Earth and, consequently, a low rate of dissipation (see Section 3). However,this finding should be considered with due caution because of the large uncertainties associatedwith the physical parameters involved and the very strong dependences assumed in the model.For example, the estimated gradients of the Maxwell time at this (
T, P ) phase point are − . − and 14 d GPa − . The resulting value may be lower by orders of magnitude ifthe actual solidus temperature is somewhat higher or if the actual pressure is somewhat lower.Second, the colder outer layers of the mantle can provide most of the dissipation if the pressureis significantly lower. Without detailed profiles of the temperature and pressure with depth, it isnot possible to predict which layer is responsible for the bulk of the dissipation. Unfortunately,these conditions are not accurately known even for the Earth and are completely subject ofspeculative calculation for exoplanets. In particular, the solidus temperature versus pressurewas experimentally estimated by Hirschmann (2000), but the model is confined to the range[0 ,
10] GPa, which is hardly relevant to the tidally reacting deep layers of the mantle. A slightlymore elaborate, but very similar, solidus model is found in Takahashi (1990).Figure 7 shows the color-coded rate of dissipation as a 2D function of temperature andpressure inside an Earth-like silicate mantle. For comparison, the expected solidus points fromHirschmann (2000) are shown with filled circles. It appears that the model by Stamenkovi´c et al.(2012) definitely predicts that solidus occurs before the peak tidal dissipation. In other words,the mantle begins to melt at much lower temperatures than those required to reach the peakdissipation rates within the Maxwell model. In this scenario, the entire mantle should bepartially molten (perhaps, with a thin surface crust staying solid). This conclusion should bequalitatively valid for all the seven planets.A different picture may emerge from the model proposed by Papaloizou et al. (2017).Those authors argued that the Maxwell model is not realistic for terrestrial-type solids witha grainy structure, while the Andrade mechanism of pinning-unpinning of defects might be amore appropriate model of the tidal friction mechanics. The frequency-dependent part of theAndrade friction can be formally characterized by an “Andrade time” parameter (Efroimsky2012a). The internal friction takes place when the grains start to shift relative to each other,which requires the external stress frequency to surpass a certain threshold value, see footnote3 and a reference therein. The threshold frequency exponentially grows with temperature andalso depends on the characteristic grain size. The value of the threshold frequency remains apoorly constrained parameter. In the model used by Papaloizou et al. (2017), it is estimateddirectly from the ambient pressure and temperature, without separately modeling the viscosityand rigidity modulus (S. Padovan, 2015, priv. comm.). Our calculations show that it is possibleto achieve the peak dissipation before the solidus for the inner TR1 planets (see Figure 7, right)if a small grain size is assumed. For this calculation, we assumed a size of 1 µ m (1 micron),which is probably too small. Larger grains would also place the peak dissipation behind the 16 –Table 3. Probabilities of the TRAPPIST-1 planets’ capture in higher-than-synchronousspin–orbit resonances, provided all these planets’ Maxwell time and dynamical triaxialityassume the values τ M = 1 d and ( B − A ) /C = 1 . × − .Planet Probability3:2 2:1TR1-1 1 0 . .
02 0TR1-3 0 .
681 0 . .
693 0 . .
05 0TR1-7 0 .
017 0Fig. 7.— The rate of tidal dissipation as a function of temperature and pressure, for planetTR1-1. The incidence of solidus is shown with filled circles. Left: viscosity model ofStamenkovi´c et al. (2012). Right: Andrade-type model of Papaloizou et al. (2017). 17 –Table 4. Explanation of notationsNotation Description R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . the radius of a planet θ . . . . . . . . . . . . . . . . the polar rotation angle of a planet˙ θ . . . . . . . . . . . . . . . . . the polar rotation rate of a planet θ . . . . . . . . . . . . . . . . the polar rotation angle of the star˙ θ . . . . . . . . . . . . . . . . . . the polar rotation rate of the star M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . the mass of a planet M . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . the mass of the star a . . . . . . . . . . . . . . . . . . . . . the semimajor axis of a planet r the instantaneous distance of a planet from the star ν . . . . . . . . . . . . . . . . . . . . . . . the true anomaly of a planet e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .orbital eccentricity M . . . . . . . . . . . . . . . . . . . . . the mean anomaly of a planet i . . . the inclination of the orbit on a planet’s equator i . . . the inclination of the orbit on the star’s equator A ≤ B < C . . . . . . . . . . . . . . . . . the moments of inertia of a planet ξ . . . . . . . . . . . . . . . . . . . the moment of inertia coefficient A . . . . . . . . . . . . . . the dimensionless rigidity of a planet n . . . . . . . . the mean motion of a planet ( n = 2 π/P orb ) G . . the gravitational constant, = 66468 m kg − yr − η . . . . . . . . . . . . . . . the mean shear viscosity of a planet µ . . . . . . . the mean shear rigidity modulus of a planet τ M ≡ η/µ . . . . . . . . . . . . . . . . the mean Maxwell time of a planet ρ . . . . . . . . . . . . . . . . . . . . . . . the mean density of a planet 18 –solidus. We conclude from this analysis that the existing and sufficiently elaborate models ofexoplanet interiors suggest the termination of overheating by solidus, whence we expect theinner planets of TR1 to be partially molten.
7. Can the inner planets be pseudosynchronous?
The combination of the models by Stamenkovi´c et al. (2012) and Hirschmann (2000) sug-gests that for the boundary condition ( T = 2200 K, P = 10 GPa), the presolidus Maxwelltime is approximately 700 d. We do not know what happens in the deeper layers of the ex-oplanet’s mantle. Despite this long time estimate, it is possible to simulate scenarios wherethe inner four planets are “softened” by partial melt over the entire depths of their mantles.Indeed, if the heat transfer is inefficient due to the absence of plate tectonics, then — for aheat capacity of 1200 J kg − K − and a very conservative dissipation rate of 10 W for asynchronously rotating TR1-1 — the interiors should be warmed up by 50 K in just 1 Myr.The rise in temperature over 1 Gyr in this “stagnant lid” scenario is 220 000, 54, 920, and 522K for planets TR1-1 through 4; and is, at most, a few K for the others. Inevitably, the deeplayers of planets 1, 3, and 4 will be consumed by the liquid core, and the outer layer will bepartially molten. Such relatively inviscid, easily deformable planets are no longer governed bythe solid-mineral rheology, but perhaps are more accurately described by the Maxwell modelin the “semiliquid” regime, characterized by the finite rigidity of the outer shell and a shortMaxwell time (Makarov 2015). It can also be viewed as an approximation for a “solidus-state”terrestrial exoplanet marginally heated by tidal forces.Cold and inviscid Earth-like planets cannot be captured in the so-called pseudosyn-chronous equilibrium, because the tidal quality ratio k /Q declines in magnitude with in-creasing perturbation frequency in the vicinity of a spin-orbit resonance (Makarov & Efroimsky2013). On the other hand, within the Constant Time Lag model (where the quality ratio k /Q is linear in the perturbation frequency), pseudosynchronism is inevitable for any nonzero ec-centricity (Mignard 1980; Hut 1981; Murray & Dermott 1999; Williams & Efroimsky 2012;Rodr´ıguez et al. 2008). A semimolten state may be implemented between these two extremes,with the Maxwell time being short enough for the peak dissipation frequency to be higherthan the perturbation frequency. As was shown by Makarov (2015), the main condition forpseudosynchronism is that the magnitude of the secular acceleration caused by the residualtidal dissipation should be greater than the amplitude of free libration (taking the maximumin the vicinity of the 1:1 spin-orbit resonance). Practically, this means that the eccentricityshould be large enough and the triaxiality small enough for the planet to retain a slightly fasterequilibrium rotation rate than the synchronous value. The two conditions can be presented as τ M . n (1 + A ) , (2) e > . σ , (3)where σ = ( B − A ) /C is the triaxiality parameter computed from the principal moments of 19 –inertia. The dimensionless rigidity A is given by A = 578 π µG ρ R (4)where G is the Newton gravitational constant, while R , ρ , and µ are the radius, the meandensity, and the mean shear rigidity of the body, correspondingly.Condition (2) results in the upper limits of (0 . , . , . , . , . , . , .
82) d forplanets 1 through 7, respectively. Although much shorter than our estimated solidus Maxwelltime at 10 GPa (700 d, cf. Section 6), such values cannot be precluded for partially meltedmantles. Maxwell times shorter than 1 d can be expected for semiliquid ice satellites of theSolar System, e.g., Io or Enceladus.Condition (3) is more definitive. The moments of inertia are not known, but we cancompute the upper limits from the estimated eccentricities and compare with the Earth’svalue. We find that the only planets that pass this criterion are TR1-1 and TR1-4, althoughthe upper limit for TR1-3 (1 × − ) is only a half of the Earth’s value. The same threeplanets (1, 3, and 4) have the highest probability of capture into the 3:2 resonance (whichwould preclude the pseudosynchronous state). However, according to Figure 6, this can onlyhappen at much longer Maxwell times.Thus, planets TR1-2, 5, 6, and 7 are likely to reside in synchronous rotation states,remaining cold, inviscid, and undifferentiated. If the planets 1, 3, 4 were cold at the time ofinitial spin-down, they could be captured into a higher-than-synchronous resonance with theensuing rapid internal heating and partial meltdown.
8. Escape from the 3:2 resonance
What happens with a planet trapped in a 3:2 spin-orbit resonance, heated by tidal fric-tion? Noyelles et al. (2013) numerically investigated spin-orbit evolution scenarios for Mercuryunder variable eccentricity conditions and concluded that a planet practically never leaves aresonance, once captured. This observation was made for a “cold” Maxwell model with a fixedaverage viscosity. A close-in exoplanet locked in the 3:2 or higher resonance receives almost twoorders of magnitude more of tidal heating (Figure 4). The rise of internal temperature in com-bination with a stable pressure brings about a drop in viscosity of several orders of magnitude(Section 6). As the dominant quality peak corresponding to the 1:1 resonance shifts towardthe 3:2 resonance and the 3:2 kink spreads out, the necessary condition for the 3:2 circulationis no longer fulfilled. We then expect the planet to leave the 3:2 spin-orbit resonance by itself,without any additional perturbation, and start its descent to the 1:1 resonance.We now recall that a tidal Fourier mode ω lmpq excited in a rotating body is a linear com-bination of the mean motion n and the angular velocity (cid:5) θ of the body (Efroimsky & Makarov2013): ω lmpq ≈ ( l − p + q ) n − m ˙ θ , (5) 20 –while the tidal torque acting on the body is an infinite sum, over lmpq , of terms, each of whichis proportional to the quality ratio k l ( ω lmpq ) Q l ( ω lmpq ) . (6)This ratio, as a function of ω lmpq , has the shape of a kink, such as that seen in Figure 8around the ˙ θ = n resonance. The tidal quality ratio k l ( ω lmpq ) /Q l ( ω lmpq ) has this shape forany rheology. The shape is natural: on the one hand, the torque component must go smoothlythrough zero on crossing a spin-orbit resonance; on the other hand, the material’s reaction hassome delay and cannot respond instantaneously to forcing at the frequencies much higher thanthe Maxwell time — hence the fall-off at high frequencies (Efroimsky 2012a,b).Taking the weighted sum over all tidal modes, treating n as a slowly changing parameter,and expressing all the terms of the quality function as functions of (cid:5) θ , we obtain the overallquality function of frequency, which looks like a superposition of kinks of different amplitudesin Figure 8. In this figure, a smaller kink (corresponding to the term { lmpq } = 2201 ) isresiding on the slope of the principal, semidiurnal term. The capture into the 3:2 spin statetakes place when the small torque crosses the horizontal axis. Now consider tidal despinningof an oblate symmetrical planet with no triaxiality. If we start close to the 3:2 spin-orbitstate (but with a slightly faster initial spin), we are guaranteed to get trapped into this state,provided the smaller kink is crossing the horizontal axis. In this case, the small kink is actingas a trap, ensuring that the oblate rotator gets captured in the 3:2 spin-orbit resonance. As thetemperature increases, the Maxwell time becomes shorter, and both kinks spread out. At somepoint, the smaller kink “sinks” below the horizontal axis. The secular tidal torque becomesnegative everywhere in the vicinity of the 3:2 resonance, and the body begins despinning towardsynchronism (or pseudosynchronism, if semimolten).To confirm this scenario, we performed numerical simulations of planet TR1-1 capturedin the 3:2 spin-orbit resonance, subject to the regular triaxial torque and the tidal torque fromthe host star. The differential equation of motion for the spin of this planet was integratedfor 10000 days with a maximum step of 3 × − d. The shortening of the Maxwell time wasartificially accelerated, modeled as an exponential decline starting at τ M = 100 d at t = 0 andending at 0.1 d at t = 10000 d. All other essential parameters were kept fixed ( M = 0 . M E , R = 1 . R E , e = 0 . P = 1 . σ = 1 . × − ). Figure 9 shows the resultingtrajectory for the relative rotation velocity, ˙ θ/n . The planet spontaneously jumps out of the 3:2resonance and very quickly descends into the 1:1 resonance. Note that the rate of spin-downwas not artificially accelerated here, being defined by the magnitude of tidal dissipation.The planet breaks out of the 3:2 lock in the middle of this simulation, when the Maxwelltime is approximately equal to 8 d. This is still a very “cold” model, with the tidal qualityratio k /Q peaking at a low frequency (cf. Figure 3). As we noted before, when both theeccentricity and the triaxiality are small, a stable spin in a higher-than-synchronous resonanceis possible only when the secular tidal torque takes positive values in the close vicinity ofthe resonant frequency. The sign of the tidal torque near the 3:2 resonance is defined by the 21 – - - - r(cid:4)(cid:5)(cid:6) of rotation θ / n k l ( ω ) / Q ( ω ) Fig. 8.— Overall tidal quality function for the planet TR1-1 (b) presented as a functionof the spin rate ˙ θ treated as a slowly changing parameter. The vertical dashed lines markthe spin-orbit resonances. The principal, large kink corresponds to the semidiurnal term with { lmpq } = { } , while the small kink corresponds to the term with { lmpq } = { } . Atlower temperatures, the small kink is crossing the horizontal axis and is therefore acting as atrap, ensuring that the planet gets caught in the 3:2 spin-orbit resonance. As the temperatureincreases, both kinks spread out. As soon as the smaller kink “sinks” below the horizontalaxis, it can no longer act as a trap, and the body begins to decelerate toward synchronism (orpseudosynchronism, if it is semimolten). 22 – time (cid:9)(cid:10)(cid:11)(cid:12) θ / n Fig. 9.— Simulated rate of rotation of planet TR1-1 initially captured in the 3:2 spin-orbitresonance, with the Maxwell time exponentially dropping from 100 d to 0.1 d by the end ofthe time span.balance of the two leading inputs, those corresponding to the Fourier tidal modes ω and ω .Thus, the escape from the resonance happens when the peak acceleration related to the ω mode becomes smaller in magnitude than the secular acceleration associated with thesynchronous ω mode because of the decreasing Maxwell time. For a small e , this happens,to a good approximation, at K ( n ) = 32 A n τ M ( nτ M ) (1 + A ) + 1 = 147 A e
16 (1 + A ) . (7)The larger root of the emerging quadratic equation yields ( escape ) τ M = 4 + √ − e n e (1 + A ) . (8)This formula is only valid for e < . .
25 d.
9. Tidal dissipation in the host star
Most of the exoplanet hosts, being Gyr-old stars, rotate slower than the close-in planets(Matsumura et al. 2010). The tidal bulge raised on the surface of the star by a planet lagsbehind the direction to the planet. The situation is opposite to that in the Earth–Moon 23 –system, where the tidal bulge running across the surface of the Earth leads the direction to theMoon. The lagging tidal bulge inevitably accelerates the star’s spin and reduces the system’seccentricity and orbital distance.Photometric light curves obtained with different instruments lead to discrepant estimatesfor the TRAPPIST-1 star rotation period, such as 3.30 d (Vida et al. 2017; Luger et al. 2017),1.40 d (Gillon et al. 2016), and 0.819 d (Roettenbacher & Kane 2017). The shorter periodsmay be caused by a group of major photometric structures distributed in longitude. The longerperiod, on the other hand, seems to be in disagreement with the measurement of the projectedvelocity (Reiners & Basri 2010), which suggests a value not too different from 1 d. At anyrate, only the two inner planets TR1-1 and 2 may have a higher mean motion than the host’srotation rate. The other five planets move slowly, and their orbits may expand and becomemore eccentric. However, the action of these planets on the star is vanishingly small comparedto the inner planets. The direction of the stellar spin evolution is mostly defined by the tidalaction of the innermost planet. In the long term, this evolution will reflect on the orbits of theouter planets. Thus, the innermost planets change the orbits of the outer planets through themedium of the tidal response on the host star, and over a very long term, this action may bemore significant than the direct gravitational interaction between the planets.Using the potential expansion by Kaula (1964), we derived the orbital evolution equationsfor the special case discussed here: (1) a small eccentricity, e → i = 0 ; (3) the star’s rotation is not synchronized, ˙ θ = n .To the most significant power of e , the resulting differential equations are dadt ≃ − n a (cid:18) M M (cid:19) (cid:18) R a (cid:19) K (2 n − θ ) (9) dedt ≃ − n e (cid:18) M M (cid:19) (cid:18) R a (cid:19) × (cid:18) − K ( n − θ ) − K (2 n − θ ) + 14716 K (3 n − θ ) + 98 K ( n ) (cid:19) . (10)Here the subscript “1” refers to the star and “2” to the planet. Accordingly, M and M arethe stellar and planetary masses, while θ and ˙ θ denote the polar rotation angle and the polarrotation rate of the star. The notation K ( ω mpq ) ≡ (cid:20) k ( ω mpq ) Q ( ω mpq ) (cid:21) ( star ) (11)stands for the quality ratio k /Q of the star, which is a function of the tidal Fourier mode ω mpq excited in the star; see expression (5). In our approximation, we need to take into considerationthe inputs with ω , − = − n , ω = n , ω , − = n − θ , ω = 2 n − θ , ω = 3 n − θ .Expressions (9) and (10) are general and valid for an arbitrary rheological model. Thequality function (11) is an odd function of the Fourier mode,Sign[ K ( ω mpq ) ] = Sign( ω mpq ) .
24 –Also be mindful that, irrespective of the tidal model, K (0) = 0 . This condition warrants thecrossing of spin-orbit resonances without experiencing an abrupt change of the tidal torquevalue.We still do not know the sign of the principal (semidiurnal) tidal mode 2 n − θ . In otherwords, we do not know whether the star rotates faster or slower than the orbital motion ofTR1-1. Depending on that sign, the inner planet may spiral in or out, per equation (9). Sincethe other planets are locked in a complex chain of MMRs, they are expected to follow suit. Thecoefficients in front of the quality functions for planets TR1-1 and 2 are 0.11 and 0.03 m d − .Tidal dissipation is most efficient when K ( ω mpq ) attains its highest possible values, whichare of the order of unity. In this case, the shortest characteristic times of the distance evolution a/ | ˙ a | are of the order of 10 – 10 yr. Thus, in principle, the tides on the host star might haveimpacted the orbital configuration had the tidal quality ratio assumed sufficiently high values.It is, however, commonly believed that the values of k /Q for stars are orders of magnitudesmaller than unity. If this assumption is true, then the role of the stellar tides in the orbitaldynamics is small, if any.The situation with the sign of de/dt is even less certain, because in equation (10) fourtidal Fourier modes provide contributions of the order of e . The overall coefficients in frontof the linear combination of the quality terms amount to 1 . × − and 1 . × − d − for planets TR1-1 and 2, respectively. The corresponding characteristic times e/ | ˙ e | , to theprecision of the unknown quality factors, are also of the order of 10 yr. It is almost certainthat the sign of this derivative is negative, i.e., that the tidal action tends to circularize theorbit. We conclude again that the tidal effect on the host star is unlikely to impact much theorbital configuration, unless the value of the tidal quality function of the stellar material ishigh — which is against our expectations.Finally, the angular acceleration of the star due to the tidal dissipation in the star is, toa good approximation in this case, given by¨ θ ≃ n ξ (cid:18) M M (cid:19) (cid:18) R a (cid:19) K (2 n − θ ) , (12)with ξ denoting the moment of inertia coefficient (equal to 0.4 for a uniform sphere). Thecorresponding characteristic times n/ | ¨ θ | of the spin-down (or spin-up for n > ˙ θ ) are consider-ably longer, being of the order of 1 Gyr at a minimum. This simple calculation demonstratesthat the planets are too small to appreciably change the rate of rotation of the TRAPPIST-1star. The observed quick spin of the star therefore represents a puzzle, unless the star is quiteyoung.
10. Orbital evolution due to the tides in a planet
Following the same procedure, we start with the Kaula expansion for the disturbing po-tential, inserted in Lagrange orbital equations. This gives us the following approximations for 25 –the rates of a and e caused by the tidal dissipation in a planet: dadt ≃ − n a e (cid:18) M M (cid:19) (cid:18) R a (cid:19) K ( n ) (13) dedt ≃ − n e (cid:18) M M (cid:19) (cid:18) R a (cid:19) K ( n ) (14)As before, M and M are the stellar and planetary masses. The notations θ and ˙ θ stand forthe polar rotation angle and the polar rotation rate of the planet, while K ( n ) = k ( n ) /Q ( n )is the tidal quality ratio (the tidal quality function) of the planet, taken at the Fourier tidalmode equal to 1 n . These expressions can be used for the specific configuration when (1) theorbital eccentricity is small; (2) the orbital obliquity i on the planetary equator is negligible;(3) the planet is synchronized ( ˙ θ = n ). Note that the rate da/dt is now proportional to e ,because in the 1:1 spin-orbit resonance all the terms of the zeroth order in e become zero.The greatest uncertainty in applying this formula to the TR1 system comes from theunknown value of the tidal quality functions of the star and of the planets. To the precision ofthese unknown values, the characteristic times of the orbital decay, a/ | ˙ a | , amount to 0.3 Myr,2 Gyr, 0.1 Gyr, and 0.2 Gyr for planets TR1-1 through 4, respectively. It may seem puzzlingthat the tides on the innermost planet may be much more efficient at shrinking the orbit thanthe tides on the star, albeit the explicit proportionality is to the fifth power of the perturbedbody’s radius. The reason for this is the mass ratio, which is inverse and outweighs the impactof the smaller bulge. Qualitatively, this happens because a smaller perturber’s mass generatesa smaller bulge and tidal torque, but the angular momentum gets distributed to the totalmass of the system anyway. It should be noted here that in all the preceding derivations, weassumed M ≫ M . Unlike the case with stellar tides, a planet’s quality ratio k /Q can bequite large in this case, up to the absolute maximum of 3 / ? ), leading to rathercomplex dynamical scenarios, which are outside of the scope of this project.The situation is different for the rate of circularization caused by the tides on the synchro-nized TR1 planets. Even though the rate of dissipation is at a global minimum in this state(except for the possible pseudosynchronism of semiliquid bodies; see Levrard 2008; Wisdom2008), the derivative de/dt is still O ( e ), and takes rather large values for the inner planets.Table 5, third column, gives the characteristic times of circularization, computed for τ M = 1 d.The innermost planet, left to itself, would have completely circularized in a matter of 30 000years. We, however, know that this planet has a small but measurable eccentricity. Almostcertainly, the finite eccentricity of the TR1 planets is excited by the chaotic gravitational in-teractions boosted by the MMR chains. We have seen in our numerical simulations that thesystem breaks up when a planet abruptly increases the eccentricity and plunges into the star. 26 –Hence, the tidal damping of eccentricity may be pivotal in the long-term survival of the sys-tem, providing especially strong stabilization to the most vulnerable inner planets. Still, thelifetime trend for the system is to shrink. The relatively short decay time for the TR1-1 planetis probably not implemented in reality, because it is not possible for it to leave the 5:3 MMRwith TR1-2, and other, weaker resonances.
11. Direction for further research: eccentricities coupling and its possible rolein stabilizing the system
The formulae for da/dt and de/dt presented in the above two sections are those of atwo-body problem. In a multibody setting, however, various mutual interactions come intoplay. In the absence of mean-motion resonances, the long-term evolution of a conservativemultiplanetary system is described, in the first order, by a system of the Laplace–Lagrangelinear secular equations (Laskar 1990). Within this approximation, inclinations and eccen-tricities are decoupled. Orbital elements become coupled when dissipative phenomena (tides)are taken into account. In this situation, one should consider the eigenmodes of eccentricity,as was done in Greenberg & Van Laerhoven (2011) and Laskar et al. (2012). Things are ele-vated to a further level of complexity by MMRs. While on physical grounds we expect thatunder MMR the couplings are strengthened, mathematical derivation of the Laplace–Lagrangesystem becomes more difficult. Its integration will become problematic because, as we havefound, the behavior of eccentricities demonstrates chaos.While these effects are beyond the scope of our current paper, we wish to point them out,because they may provide a key to how the TRAPPIST-1 system ended up on its island ofstable solutions.
12. Discussion and conclusions
We have carried out the numerical integration of the long-term evolution of the TRAPPIST-1 system, starting with a previously suggested stable configuration. Having confirmed thelong-term stability of this configuration, we also have found that the motion of planets ischaotic. The eccentricity values vary within certain ranges in an apparently random way. Us-ing the robust mean values of eccentricity, the expected spin-orbit states have been computedin the probabilistic sense. The rates of tidal dissipation and tidal evolution of orbits have beenestimated, assuming an Earth-like rheology for the planets.The TRAPPIST-1 system may not be an oddball in the world of planetary systems, buta rather common and natural outcome of star formation and N -body dynamics, where tidalinteractions undoubtedly play a big role. The orbital motion of the tightly packed planetsis distinctly chaotic, with a “short-fused” Lyapunov time — and this picture may be normalfor exoplanet systems. The complex multiplet of MMRs, where each of the seven planetshas a commensurable orbital motion with each of the other six planets, suggests a highly 27 –organized structure unseen in the Solar System. The challenge is to explain how this degreeof organization and complexity can naturally emerge from the primary stochastic material,following the principles of classical mechanics and thwarting the flow of entropy. The MMRchains have a regularization effect on the orbital dynamics, damping the chaotically emergingperturbations and keeping eccentricity at low levels. This prevents the planets from comingat close range and directly interacting with each other. Are there other regularization agents?Could the secular tidal interaction have made the long-term orbital evolution more orderly andconstrained?Working within the Maxwell rheological model for the planets, we find that planets TR1-1,3, and 4, with initial fast rotation in the prograde sense and cold inviscid mantles, are certainto be originally captured into the 3:2 or even higher spin-orbit resonance. At the current valuesof their orbital eccentricities, planets 6 and 7 could have been captured too, if their averageviscosity was high enough. In higher-than-synchronous resonances, the rotational dynamicswith tides is quite complicated and nonlinear, with higher levels of forced libration involved(Frouard & Efroimsky 2018) and even with the possibility of additional secondary attractorsin the phase space (Bartuccelli et al. 2017). The rate of the kinetic-energy dissipation is O ( e )for all higher-than-synchronous resonances, which brings about high rates of internal heating.For TR1-1, for example, the heating rate in the 3:2 spin state is at least an order of magnitudehigher than in the synchronized state and can reach 10 . W when the Maxwell time valuecorresponds to the peak of the quality function (see Figure 3). The emerging equilibriumsurface flux amounts to a few MW m − . For comparison, it is equivalent to a small power-generating gas turbine on each square meter of the planet surface. Enceladus, one of the mostimpressive tidal power plants in the Solar System, generates ∼ W, which corresponds tothe mean surface flux of ∼ − W m − (Kamata & Nimmo 2017; Efroimsky 2018). The poweroutput on more distant planets of TRAPPIST-1 is less extreme, but quite significant, too. Theinternal temperature inevitably rises, as the tidal heat cannot be efficiently transferred to thesurface and radiated out.As the viscosity η of silicate minerals is a strong function of temperature, the planet’sinteriors have lower values of η — which further increases the dissipation rate. Runawayheating appears inevitable (Bˇehounkov´a et al. 2011), turning the initially rocky planets intoballs of liquid magma. We posit, however, that there are heating self-regulation mechanismsterminating the runaway process much sooner. First, at the low orbital frequencies observedfor TR1 system, and with presumably small degrees of triaxial deformation, the condition ofcapture into the 3:2 spin state vanishes abruptly, and the planet jumps out of the resonance. Itis bound to spin down very quickly, getting trapped in the permanent 1:1 spin-orbit resonance.The planets can be in the semiliquid state during this transition, so the pseudosynchronismequilibrium is possible in principle; however, only the inner planet passes the eccentricitycriterion of pseudosynchronism. The rate of dissipation drops by an order of magnitude (cf.Figure 4). This is still quite high for the inner four planets, and they continue to warm up. Theultimate stop to this process is set by the instance of solidus, when part of the material startsto melt. Both the viscosity and the rigidity modulus take a dive by many orders of magnitude, 28 –resulting in ultra-short Maxwell times and, hence, a much lower dissipation rate. According tothe existing models, the solidus point should be ∼ b , d , and e ) were captured in the 3:2 or higher spin-orbitresonances during the initial spin-down, the end product of this scenario is a TR1 system whereall the planets rotate synchronously with their orbital motion, minimizing the tidal heating.A possible exception may be the closest planet b , which may, depending on its rheology, betrapped in pseudosynchronism. The inner four planets are likely to have massive molten coresand hot mantles with a significant degree of melt. These conditions help them to generate asmuch tidal heat as can be transferred to the surface and irradiated into space. The limited tidaldissipation in the planets means a limited impact on the orbital dynamics and stability of thesystem. Owing to the tidal dissipation inside the star, the orbits are expected to keep shrinkingeven in the synchronous state, though at a moderate or low rate (unless the tidal quality ofthe star is much higher than what is commonly assumed). The eccentricity of the planets, onthe other hand, can go down quickly enough on the lifetime scale. We conclude that, for long-term stability to be maintained, the planet-planet gravitational interaction should keep theeccentricities’ values within finite ranges. We also conclude that the tides on the inner planetscan provide an additional regularization mechanism adding to the stability of the configuration. REFERENCES
Bou´e, G.; Correia, A.C.M.; and Laskar, J. 2016. “Secular and tidal evolution of circumbinarysystems.” CeMDA, 126, 31Bartuccelli, M.; Deane, J.; and Gentile, G. 2017. “Periodic and quasi-periodic attractors forthe spin-orbit evolution of Mercury with a realistic tidal torque.” MNRAS, 469, 127Bˇehounkov´a M.; Tobie, G.; Choblet, G.; and ˇCadek, O. 2011. “Tidally Induced ThermalRunaways on Extrasolar Earths: Impact on Habitability.” ApJ, 728, 89Bonfils, X., et al. 2013. “The HARPS search for southern extra-solar planets. XXXI. TheM-dwarf sample.” A&A, 549, A109Chambers J.E. 1999. “A hybrid symplectic integrator that permits close encounters betweenmassive bodies.” MNRAS, 304, 793Correia, A.C.M., & Laskar, J. 2004. “Mercury’s capture into the 3 / Fundamentals of celestial mechanics.
MacMillan, New YorkDelfosse, X., et al. 2013. “The HARPS search for southern extra-solar planets. XXXIII. Super-Earths around the M-dwarf neighbors Gl 433 and Gl 667C.” A&A, 553, A8Dobrovolskis, A.R., Ingersoll, A.P. 1980. “Atmospheric tides and the rotation of Venus. I -Tidal theory and the balance of torques.” Icarus, 41, 1Doremus, R.H. 2002. “Viscosity of silica.”
Journal of Applied Physics , 92, 7619Efroimsky, M. & Makarov, V.V. 2014. “Tidal dissipation in a homogeneous spherical body. I.Methods.” ApJ, 795, 6. ERRATA: ApJ, 819, 86 (2016).Efroimsky, M. & Makarov, V.V. 2013. “Tidal Friction and Tidal Lagging. Applicability Limi-tations of a Popular Formula for the Tidal Torque.” ApJ, 764, 26Efroimsky, M. 2012a. “Bodily tides near spin-orbit resonances.” CeMDA, 112, 283Efroimsky, M. 2012b. “Tidal dissipation compared to seismic dissipation: in small bodies,earths, and superearths.” ApJ, 746, 150. ERRATA: ApJ, 763, 150 (2013)Efroimsky, M. 2015. “Tidal Evolution of Asteroidal Binaries. Ruled by Viscosity. Ignorant ofRigidity.” AJ, 150, 98. ERRATA: AJ, 151, 130 (2016)Efroimsky, M. 2018. “Tidal viscosity of Enceladus.” Icarus, 300, 223Ferraz-Mello, S. 2015. “Tidal synchronisation of close-in satellites and exoplanets: II. Spindynamics and extension to Mercury and exoplanet host stars.” CeMDA, 122, 359Filippazzo, J.; Rice, E.L.; Faherty, J.K.; et al. 2015. “Fundamental Parameters and SpectralEnergy Distributions of Young and Field Age Objects with Masses Spanning the Stellarto Planetary Regime.” ApJ, 810, 158Frouard, J.; Quillen, A.C.; Efroimsky, M.; and Giannella, D. 2016. “Numerical simulation oftidal evolution of a viscoelastic body modelled with a mass-spring network.” MNRAS,458, 2890Frouard, J., & Efroimsky, M. 2017. “Tides in a body librating about a spin-orbit resonance:generalization of the Darwin-Kaula theory.” CeMDA, 129, 177Gillon, M.; Jehin, E.; Lederer, S. M.; et al. 2016. “Temperate Earth-sized planets transiting anearby ultracool dwarf star.” Nature, 533, 221 30 –Gillon, M., Triaud, A.H.M.J., Demory, B.-O., et al. 2017. “Seven temperate terrestrial planetsaround the nearby ultracool dwarf star TRAPPIST-1.” Nature, 542, 456Goldreich, P., and Peale, S. 1966. “Spin-orbit coupling in the Solar System.” AJ, 71, 425Goldreich, P., and Peale, S.J. 1968. “The Dynamics of Planetary Rotations.” ARA&A, 6, 287Greenberg, R., and Van Laerhoven, C. 2011. “Tidal Evolution of a Secularly Interacting Plan-etary System.” ApJ, 733, A8Hayes, W. 2007. “Is the outer Solar System chaotic?” Nature Physics, 3, 689Hayes, W. 2008. “Surfing on the edge: chaos versus near-integrability in the system of Jovianplanets.” MNRAS, 386, 295Heng, K., Kopparla, P. 2013. “On the stability of super-Earth atmospheres”. ApJ, 754, 60Henning, W.G., and Hurford, T. 2014. “Tidal Heating in Multilayered Terrestrial Exoplanets.”ApJ, 789, 30Henning, W.G.; O’Connell, R.J.; and Sasselov, D. 2009. “Tidally Heated Terrestrial Exoplan-ets: Viscoelastic Response Models.” ApJ, 707, 1000Hirschmann, M.M. 2000. “Mantle solidus: Experimental constraints and the effectsof peridotite composition.” Geochemistry, Geophysics, Geosystems, 1, article id.2000GC000070Hut, P. 1980. “Stability of tidal equilibrium.” A&A, 92, 167Hut, P. 1981. ‘Tidal evolution in close binary systems.” A&A, 99, 126Karato, S.-i., and Spetzler, H. A. 1990. “Defect Microdynamics in Minerals and Solid-StateMechanisms of Seismic Wave Attenuation and Velocity Dispersion in the Mantle.” Re-views of Geophysics, Vol. , pp. 399 - 423Kaula, W.M. 1964. “Tidal Dissipation by Solid Friction and the Resulting Orbital Evolution.”Reviews of Geophysics, 2, 661Kamata, S., and Nimmo, F. 2017. “Interior thermal state of Enceladus inferred from theviscoelastic state of the ice shell.” Icarus, 284, 387Lambeck K. 1980. The Earth’s variable rotation . Cambridge University Press, New YorkLaskar, J.; Bou´e, G.; and Correia, A. 2012. “Tidal dissipation in multi-planet systems andconstraints to orbit-fitting.” A&A, 538, A105Laskar, J. 1990. “The chaotic motion of the solar system: A numerical estimate of the size ofthe chaotic zones.” Icarus, 88, 266 31 –Leconte, J., Wu, H., Menou, K., Murray, N. 2015. “Asynchronous rotation of Earth-massplanets in the habitable zone of lower-mass stars.” Science, 347, 632Levrard B. 2008. “A proof that tidal heating in a synchronous rotation is always larger thanin an asymptotic nonsynchronous rotation state.” Icarus, 193, 641Luger, R.; Sestovic, M.; Kruse, E.; et al. “A seven-planet resonant chain in TRAPPIST-1.”2017, NatAs, 1, 129Makarov, V.V. 2012. “Conditions of passage and entrapment of terrestrial planets in spin-orbitresonances.” ApJ, 752, 73Makarov, V.V. 2015. “Equilibrium Rotation of Semiliquid Exoplanets and Satellites.” ApJ,810, 12Makarov, V.V., and Efroimsky, M. 2013. “No pseudosynchronous rotation for terrestrial planetsand moons.” ApJ, 764, 27Makarov, V.V.; Berghea, C.; and Efroimsky, M. 2012. “Dynamical Evolution and Spin-OrbitResonances of Potentially Habitable Exoplanets: The Case of GJ 581d.” ApJ, 761, 83Makarov, V.V., and Berghea, C. 2014. “Dynamical Evolution and Spin-orbit Resonances ofPotentially Habitable Exoplanets. The Case of GJ 667C.” ApJ, 780, 124Margot, J.-L.; Peale, S.J.; Jurgens, R.F.; Slade, M.A.; and Holin, I.V. 2007. “Large LongitudeLibration of Mercury Reveals a Molten Core.”
Science. , 316, 710Matsumura, S.; Peale, S.J.; and Rasio, F.A. 2010. “Tidal Evolution of Close-in Planets.” ApJ,725, 1995Mignard, F. 1980. “The Evolution of the Lunar Orbit Revisited. II.” M&P, 23, 185Murray, C.D., & Dermott, S.F. 1999.
Solar System Dynamics . Cambridge University Press,Cambridge UKNoyelles, B.; Frouard, J.; Makarov, V. V.; & Efroimsky, M. 2013. “Spin-orbit evolution ofMercury revisited.” Icarus, 241, 26 - 44Padovan, S.; Margot, J.-L.; Hauck, S.A.; Moore, W.B.; and Solomon, S.C. 2014. “The tides ofMercury and possible implications for its interior structure.” JGRE, 119, 850Papaloizou, J.C.B.; Szuszkiewicz, E.; Terquem, C.A. 2017, arXiv:1711.07932v1Peale, S.J. 2005. “The free precession and libration of Mercury.” Icarus, 178, 4Peale, S.J., and Cassen, P. 1978. “Contribution of tidal dissipation to lunar thermal history.”Icarus, 36, 245 32 –Quarles, B.; Quintana, E.V.; Lopez, E.; Schlieder, J.E.; and Barclay, T. 2017. “Plausible Com-positions of the Seven TRAPPIST-1 Planets Using Long-term Dynamical Simulations.”ApJ, 842, L5Reiners, A., and Basri, G. 2010. “A Volume-Limited Sample of 63 M7-M9.5 Dwarfs. II. Activity,Magnetism, and the Fade of the Rotation-Dominated Dynamo.” ApJ, 710, 924Rodr´ıguez, A.; Ferraz-Mello, S.; and Hussmann, H. 2008. “Tidal friction in close-in planets.”In: Y.S. Sun, S. Ferraz-Mello, and J.L. Zhou (Eds.)
Exoplanets: Detection, Formationand Dynamics. Proceedings of the IAU Symposium No 249, pp. 179 - 186Roettenbacher, R.M., and Kane, S.R. 2017. “The Stellar Activity of TRAPPIST-1 and Con-sequences for the Planetary Atmospheres.” ApJ, 851, 77Shoji, D., and Kurita, K. 2014. “Thermal-Orbital Coupled Tidal Heating and Habitability ofMartian-sized Extrasolar Planets around M Stars.” ApJ, 789, 3Stamenkovi´c, V.; Noack, L.; Breuer, D.; and Spohn, T. 2012. “The Influence of Pressure-dependent Viscosity on the Thermal Evolution of Super-Earths.” ApJ, 748, 41Stein, C., and Hansen, U. 2013. “Arrhenius rheology versus Frank-Kamenetskii rheology—Implications for mantle dynamics.” Geochemistry, Geophysics, Geosystems, 14, 2757Storch, N.I., and Lai, D. 2014. “Viscoelastic tidal dissipation in giant planets and formationof hot Jupiters through high-eccentricity migration.” MNRAS, 438, 1526Takahashi, E. 1990. “Speculations on the Archean mantle: Missing link between komatiite anddepleted garnet peridotite.” J. Geophys. Res., 95, 15941Tamayo, D.; Rein, H.; Petrovich, C.; and Murray, N. 2017. “Convergent Migration RendersTRAPPIST-1 Long-lived.” ApJ, 840, L19Terquem, C., and Papaloizou, J.C.B. 2007. “Migration and the Formation of Systems of HotSuper-Earths and Neptunes.” ApJ, 645, 1110Vida, K.; K˝ov´ari, Zs.; P´al, A..; Ol´ah, K.; and Kriskovics, L. 2017. “Frequent Flaring in theTRAPPIST-1 System—Unsuited for Life?.” ApJ, 841, 124Walterov´a, M.; and Bˇehounkov´a M. “Tidal effects in differentiated viscoelastic bodies: a nu-merical approach.” 2017, CeMDA, 129, 235Wang, Y., Tian, F., Hu, Y. 2014. “Climate Patterns of Habitable Exoplanets in EccentricOrbits Around M Dwarfs.” ApJ, 791, L12Williams, J.G., & Efroimsky, M. 2012. “Bodily tides near the 1:1 spin-orbit resonance. Cor-rection to Goldreich’s dynamical model.” CeMDA, 114, 387Wisdom J. 2008. “Tidal dissipation at arbitrary eccentricity and obliquity.” Icarus, 193, 637 33 –Zahnle, K.J.; Lupu, R.; Dobrovolskis, A. and Sleep, N.H. 2015. “The tethered Moon.” E&PSL,427, 74Zlenko, A.A. 2015. “A celestial-mechanical model for the tidal evolution of the Earth-Moonsystem treated as a double planet.” Astronomy Reports, 59, 72
This preprint was prepared with the AAS L A TEX macros v5.2.
34 –Table 5. Characteristic times of the orbital decay ( a/ | ˙ a | ) and circularization ( e/ | ˙ e | ) in Myr,assuming τ M = 1 d.Planet a/ | ˙ a | e/ | ˙ e | Myr MyrTR1-1 5 . . > . . . > > >>