The determination of planetary structure in tidally relaxed inclined systems

[Abridged] The recent discovery of a transiting planet on a non-circular orbit with a massive highly eccentric companion orbiting HAT-P-13 offers the possibility of probing the structure of the short-period planet. The ability to do this relies on the system being in a quasi-equilibrium state in the sense that the eccentricities are constant on the usual secular timescale, and decay on a timescale which is much longer than the age of the system. Since the equilibrium eccentricity is effectively a function only of observable system parameters and the unknown Love number of the short-period planet, the latter can be determined with accurate measurements of the planet's eccentricity and radius. However, this analysis relies on the unlikely assumption that the system is coplanar. Here we generalize our recent analysis of this fixed-point phenomenon to mutually inclined systems and show that the fixed point of coplanar systems is replaced by a limit cycle, with the average value of the eccentricity decreasing and its amplitude of variation increasing with increasing mutual inclination. This behaviour significantly reduces the ability to unambiguously determine the Love number of the short-period planet if the mutual inclination is higher than around 10^o. We show that for Q-values less than 10^6, the HAT-P-13 system cannot have a mutual inclination between 54 and 126^o because Kozai oscillations coupled with tidal dissipation would act to quickly move the inclination outside this range, and that the behaviour of retrograde systems is the mirror image of that for prograde systems. We derive a relationship between the equilibrium radius of the short-period planet, its Q-value and its core mass, and show that given current estimates of e_b and the planet radius, the HAT-P-13 system is likely to be close to coplanar [...]

, decreasing and its amplitude of variation increasing with increasing mutual inclination. This behaviour significantly reduces the ability to unambiguously determine the Love number of the short-period planet if the mutual inclination is higher than around 10 o . (2) We show that for Qvalues less than 10 6 , the HAT-P-13 system cannot have a mutual inclination between 54 and 126 degrees because Kozai oscillations coupled with tidal dissipation would act to quickly move the inclination outside this range, and (3) that the behaviour of retrograde systems is the mirror image of that for prograde systems in the sense that (almost) identical limit cycles exist for a given mutual inclination and π minus this value. (4) We derive a relationship between e (av) b , the equilibrium radius of the shortperiod planet, its Q-value and its core mass, and show that given current estimates of e b and the planet radius, as well as the lower bound placed on the Q-value by the decay rate of e (av) b , the HAT-P-13 system is likely to be close to prograde coplanar, or have a mutual inclination between 130 o and 135 o . Lower rather than higher core masses are favoured. (5) An expression for the timescale for decay of the mutual inclination is derived, revealing that it evolves towards a non-zero value as long as e b > 0 on a timescale which is much longer than the age of the system. (6) We conclude with a scattering scenario for the origin of the HAT-P-13 system and show that almost identical initial conditions can result in significantly different outer planet eccentricities, stellar obliquities and planet radii. The implications for systems with high stellar obliquities such as HAT-P-7 and WASP-17 are briefly discussed.

INTRODUCTION
Transiting systems offer the opportunity to determine a wide variety of system parameters, including the mass of the transiting planet, its orbital eccentricity, the inclination of its orbit to the line of sight, its radius and hence mean density, and the sky projection of the angle between its orbit normal and the stellar spin axis. Many other system parameters are potentially measurable (Winn 2009), one of them being the Love number of the transiting planet if the system parameters are favourable (Wu & Goldreich 2002). The recent discovery of the HAT-P-13 system (Bakos et al. 2009) provides us with such a system as was recently pointed out by Batygin, Bodenheimer & Laughlin (2009). A reliable estimate of a planet's Love number in turn allows one to say something about the presence or otherwise of a planetary core, and hence about the formation mode; in the former case the planet is likely to have been formed via core accretion of solid material and subsequent accretion of a massive atmosphere (Pollack et al. 1996), while the absence of a core supports the gravitational collapse hypothesis of giant planet formation (Boss 1997). The HAT-P-13 system consists of a 0.85 Jupiter-mass planet (planet b) in an almost circular 2.9 day orbit about a 1.2 solar-mass star, and a companion with a minimum mass of 15.2MJ in a 428 day highly eccentric orbit (planet c). Table 1 lists the relevant parameters of the system, with data taken from Bakos et al. (2009) who performed a simultaneous fit of HATNet and KeplerCam photometric data and Keck spectroscopic data, a process they refer to as "global" modelling. Here P b and Pc are orbital periods, i los is the inclination of the orbit normal of planet c to the line of sight and ω b and ωc are periapse arguments. Note in particular the non-zero estimate of 0.021 ± 0.009 for the inner planet's eccentricity. The main aim of this paper is to explore the extent to which information can be gleaned from such a non-zero measurement, whether or not we know the mutual inclination of the orbits of the two planets.
The ability to determine the Love number of HAT-P-13b relies on the system being in a quasi-equilibrium state in the sense that after an initial transient phase of rapid tidal evolution, the eccentricities change on a timescale much longer than ⋆ E-mail: mardling@sci.monash.edu.au the age of the system and the apsidal lines of the inner and outer planets are aligned or anti-aligned. Moreover, the mass and eccentricity of HAT-P-13c ensure that the quasi-equilibrium eccentricity of the inner planet is significant and measurable. Some of the theory for this is developed in Wu & Goldreich (2002) for the HD83443 system, while a general theory is presented in Mardling (2007). Both of these studies, however, are for the coplanar case only. Using the HAT-P-13 system as illustration, we generalize the theory to non-coplanar systems in which the outer body dominates the total angular momentum and show that the quasi-relaxed state no longer corresponds to a fixed point in e b − η space, where η = ̟ b − ̟c with ̟ b and ̟c the longitudes of periastron, but rather corresponds to a limit cycle in this space (see, for example, Jordan & Smith 1999), with the average value of e b decreasing and its amplitude of variation increasing with increasing mutual inclination. As in the coplanar case, the relaxation timescale is around three times the circularization timescale. However, unlike the coplanar case, the rate of change of the argument of periastron of the inner planet plays an important role, with the limit cycle frequency equal to twice this quantity. In fact it is the appearance of this additional frequency which prevents the system from evolving to a fixed point, with terms which depend on it effectively acting as external forcing terms. Note that our analysis takes into account the fact that the actual mass of the outer planet increases with increasing mutual inclination, given the observed minimum mass determined via radial velocity measurements (Bakos et al. 2009).
The plan for this paper is as follows. Section 2 reviews the theory for the long-term tidal evolution of coplanar systems, Section 3 generalizes this to mutually inclined systems, including limit cycle behaviour of prograde systems (Section 3.1), the Kozai regime (Section 3.2), limit-cycle behaviour of retrograde systems (Section 3.3), the relationship between e (av) b , the equilibrium radius of planet b, its Q-value and its core mass (Section 3.4), and the timescale for decay of the mutual inclination in a relaxed system (Section 3.5). Section 4 presents scattering scenarios for the origin of HAT-P-13-like systems, including a discussion of stellar obliquity in two-planet systems (Section 4.1). Section 5 presents a summary, and Appendix A presents the orbit-averaged equations of motion for a two-planet Newtonian point-mass system up to octopole order, correct to leading order in the inner planet's eccentricity, the ratio of semimajor axes, and the sine of the inclination of the outer orbit relative to its initial orbit.

LONG-TERM TIDAL EVOLUTION OF COPLANAR SYSTEMS
In Mardling (2007), the long-term tidal evolution of short-period planets with single companions is studied. There it is shown that such systems evolve on three distinct timescales, an illustration of which is given in Figures 3 and 4 of that paper. As long as the eccentricity of the outer planet is non-zero, the eccentricities of both planets will initially execute anti-phased secular oscillations until a non-zero quasi-equilibrium value is reached, with maxima and minima of the inner planet's eccentricity given by expressions (20), (27) or (28) in Mardling (2007), the choice of which depends on the system parameters, and the corresponding values of the outer planet's eccentricity given by expression (29). The equilibrium value of the eccentricity depends on all contributions to the rate of apsidal motion of the inner planet. In Mardling (2007) only the contributions from the outer planet and the post-Newtonian terms in the star's potential were taken into account, however, as Ragozzine & Wolf (2009) point out, the contribution of the tidal bulge of a short-period planet like HAT-P-13b is significant and is included here (as are the contributions from the spin bulge of the planet and the tidal and spin bulges of the star). The equilibrium eccentricity for a coplanar system is given by where e b , a b , and m b are respectively the innermost planet's eccentricity, semimajor axis and mass, and ec, ac and mc are the corresponding values for the outer planet. Here εc = √ 1 − e 2 c and γ = γ GR * to first-order in the inner eccentricity. Here n b is the mean motion of the inner planet, c the speed of light, R b and k b the radius and Love number of the inner planet, m * , R * , k * and Ω * the star's mass, radius, Love number and spin frequency respectively, and synchronous rotation of the planet is assumed in the expression for γ spin b . The γ's are the ratios of the various contributions to the apsidal motion of the inner planet to the coplanar contribution of the outer planet. These are listed in Table 1 for the HAT-P-13 system (for mc = m min c ), together with those for the star assuming a spin period of 25 days. The Love number for the star (equal to twice its apsidal motion constant) is taken to be that for an n = 3 polytrope (Sterne 1941).
If the angle between the apsidal lines of the two planetary orbits, η, circulates rather than librates, the average eccentricity of the inner planet during the initial oscillatory phase will decrease on its own circularization timescale until its minimum eccentricity is zero, with the amplitude of oscillation remaining constant. This timescale is given to second order in the eccentricity by where Q b is the Q-value of the inner planet. Once the minimum of the inner eccentricity reaches zero, η will librate about zero or π (the choice of which depends on the system parameters), and the oscillation amplitude will reduce to zero until the inner eccentricity reaches a non-zero quasi-equilibrium value. The librating phase occurs on a timescale of 2τcirc.
Once the oscillatory phase is over, the system will evolve to the doubly-circular state on the approximate timescale where e * c is the value of ec at the beginning of this phase, and Note, however, that if the inner semimajor axis evolves appreciably on this timescale, (6) represents an upper limit only. Note also that τc is independent of the circularization timescale of the outer planet; if the latter is comparable to or shorter than τc then (6) again represents an upper bound. For the HAT-P-13 system, we have τcirc = 4 × 10 7 (Q b /10 5 ) yr and τc = 2.5 × 10 5 τcirc = 10, 000(Q b /10 5 ) Gyr, while the orbital decay timescale for planet b (Yoder & Peale 1981) is τa = e −2 b τcirc ≃ 2500 τcirc = 100(Q b /10 5 ) Gyr < τc, the latter two being much greater than the age of the system. A direct coplanar integration using the averaged code presented in Mardling & Lin (2002) (using a constant value of the radius of planet b) gives τa = 79 (Q b /10 5 ) Gyr and τc = ec/ėc = 9200 (Q b /10 5 ) Gyr, while e (eq) b decays on a timescale τe = 15 (Q b /10 5 ) Gyr. 1 The latter is consistent with the estimate obtained from (1) with ∆0 ≃ 1 + γε 3 c . The numerically determined timescales are listed in Table 1, together with the timescale for decay of the mutual inclination to its equilibrium value, τi, for the case that it is initially 30 o (see Section 3.5).
The estimate for τe places a lower bound on the value of Q b such that Q b /10 5 > (age/τe,5)/ ln[e b (0)/e b (age)], where τe,5 = τe(Q b = 10 5 ), e b (0) is the "initial" value of e b and e b (age) is its value now. Given the lower bound for the age of the system, using τe,5 = 15 Gyr and taking e b (0) = 1 gives Q b > 7000, while e b (0) = 0.1 gives Q b > 1.8 × 10 4 . Note that our estimate for τe is more than twice that of Batygin, Bodenheimer & Laughlin (2009) who estimate τe ≃ 6(Q b /10 5 ) Gyr, making their lower bound for Q b a factor 15/6 higher.
In general, coplanar systems for which one may say something about the structure of the short-period planet will have the following characteristics: (1) The circularization timescale of planet b will be (considerably) less than one third the age of the system, ensuring that the system is sufficiently relaxed; (2) The timescales on which the system becomes doubly circular and the orbit decays, τc and τa respectively, will both be longer than the age of the system, and (3) The equilibrium eccentricity and the radius of the short-period planet will be measurable.

INCLINED SYSTEMS
In this Section the analysis of Mardling (2007) is generalized to non-zero mutual inclination. While the HAT-P-13 system is used for illustration, the analysis may be applied to any system for which e b ≪ 1, a b /ac ≪ 1 and most of the angular momentum of the system resides in the outer orbit, the latter ensuring that the (sine of the) angle between the invariable plane and the outer orbit remains small.
In order to highlight the differences between coplanar and inclined systems, we begin by presenting the results of an integration of a HAT-P-13-like system for which the mutual inclination is 30 o , and the mass of the outer body is taken to be the observed minimum mass divided by the cosine of the mutual inclination. The stellar obliquity is such that the star's spin 1 Note that Mardling & Lin (2002) code assumes a forcing-frequency-independent Q-value, a legacy of the pioneering work of Goldreich & Soter (1966) whose supporting argument was based on the constancy of the Q-value over a wide range of forcing frequencies for the Earth, ie., for a solid body. For gaseous bodies, it seems more reasonable to use the "constant time-lag" concept originally devised by Darwin in which the lag angles of individual tidal components are proportional to their forcing frequencies. The latter is equivalent to the concept of fluid stress and its associated dissipation. Note, however, that for small eccentricities (which is the case in the present study), there is little difference between the formulations. axis is aligned with planet c's orbit normal. The integration is done using the "averaged" code presented in Mardling & Lin (2002) which includes accelerations due to spin and tidal bulges of both the star and inner planet (both quadrupole and dissipative) as well as the relativistic potential of the star. Orbit-averaged expressions are used for the evolution of the inner orbit, while the outer orbit is integrated directly. No assumptions are made about the magnitude of any quantity except the ratio of semimajor axes which is assumed to be small. For now we suppress the evolution of the planetary radius and assume it is constant at its currently observed value. Figure 1 shows the analogue of Figure 3 in Mardling (2007). The initial values of η = ̟ b − ̟c and e b are 70 o and 0.05 respectively. Rather than relaxing to a fixed point in e b − η space on a timescale of 3τcirc, the system relaxes to a limit cycle on the same timescale. Note that the Q-value of planet b is given an artificially low value of 50 in order to clearly demonstrate capture onto the limit cycle. While the modulation period is correct, the relaxation process would normally take Q b /50 as long as shown here. The circulatory phase is shown in black in each panel, the pre-capture libratory phase is shown in green and the limit-cycle phase is shown in red. Panel (c) shows an enlargement of the limit cycle centred on (e b , η) = (e , 22π), together with the theoretical limit cycle derived below (blue circle). Here e (av) b is the inclined-system analogue of the fixed-point value of e b for coplanar systems, e (compare equations (1) and (19)). While only around 1.5 relaxation timescales (ie., 4.5τcirc) are shown here, the system was integrated for 6 relaxation timescales during which e (av) b and i b decreased by less than 0.0005 and 0.05 o respectively. Note that 20 τcirc = 2.4(Q b /10 5 ) Gyr, comparable to the estimated age of the system for realistic values of Q b .
In order to guide the following analysis, we now present the results of several integrations of relaxed HAT-P-13-like systems with identical initial conditions except that the mutual inclination is varied between 0 o and 50 o (higher inclinations will be discussed in Sections 3.2 and 3.3). The initial value of e b is taken as e , and the apsidal lines are taken to be aligned. Initial angles are measured with respect to the initial orbit of planet c, that is, the relative inclination and longitude of planet b are specified, with the zero in longitude coinciding with the apse of planet c. Since planet c's orbit contains 98% of the total angular momentum, it coincides approximately with the invariable plane and there is very little change in its inclination as Figure 3(e) shows. Following Batygin, Bodenheimer & Laughlin (2009) we take the Love number of planet b to be 0.3, representative of a range of planetary structures according to their Table 1, and its radius of gyration, r gb , to be , given by equation (19). The artificially low Q-value allows the transient behaviour to die away quickly. Note that the modulation frequency and amplitude are independent of Q b and are given by 2Wω and A/2Wω, with Wω and A defined in equations (16) and (15) (19) is not reliable for this case. 0.26R b , appropriate to an n = 1 polytrope. For this set of experiments, an artificially low Q-value of 10 is used to hasten the evolution towards the relaxed state. The stellar obliquity relative to the invariable plane normal, θ * , is set to zero initially, and the stellar spin period is 25 days. The stellar Love number and radius of gyration rg * are taken to be 0.03 and 0.076R * respectively, the latter appropriate to an n = 3 polytrope.  Figure 3 shows the evolution of the other orbital elements for each value of i b . Panel (a) shows the evolution of the argument of periastron, each offset by 360 o for clarity. Note in particular that for these initial inclinations,ω b is approximately constant, and is given by (16). Panel (b) shows the libratory behaviour of η. Unlike the coplanar case which evolves to a constant value such that sin η = −WT /|Wq| (and hence is proportional to Q −1 b ), the quasi-relaxed state for systems with non-zero mutual inclination is oscillatory (although the offset is still Sin −1 [−WT /|Wq|]). The amplitude of oscillation is given by A/2Wω (times 180/π), with A and Wω defined in equations (15) and (16) respectively.
Panel (c) shows the evolution of the longitude of the ascending node of planet b for all inclinations. Note that while planet b's precession rate is proportional to mc cos i b (equation (A15)), mc = m min c /| cos i b | with m min c held constant so that the precession rates are all identical. The precession period is 7281 years. Panels (d) and (e) show the evolution of the inclinations of orbits b and c relative to the invariable plane, i b , and ic respectively, justifying our assumption in the analysis below that they are constant and that sin ic ≪ 1 (see Section 4.1 for analysis of the influence on i b of non-zero θ * ). Not shown is the variation in planet c's eccentricity; this is effectively constant with an amplitude of variation of 0.003.
Finally, panel (f) shows the stellar obliquity relative to orbit b (the angle between the star's spin axis and planet b's orbit normal), ψ * b . In Section 4.1 we show that the maxima and minima of ψ * b are |i b ± θ * |. In fact, θ * varies (the star nutates) due to torques between the misaligned stellar spin bulge and b's orbit, consistent with the variations seen in panel (f).
In light of the recent discoveries of high (sky-projected) stellar obliquities in the systems HAT-P-7 (Narita et al. Winn et al. 2009) and WASP-17 (Anderson et al. 2010), we discuss further the dynamics of stellar obliquity in Section 4.1.

Limit-cycle behaviour of prograde orbits
In Appendix A we give the orbit-averaged disturbing function for a two-planet Newtonian point-mass system up to octopole order for arbitrary planet b inclination, and from this derive the equations governing the secular evolution of the elements in the absence of perturbations. Only leading order terms in e b , a b /ac and sin ic are retained, each of which are of order 0.01 for the HAT-P-13 system whether or not the reference plane is taken as the initial plane of planet c or the invariable plane. Equations (A12), (A14), (A16) and (A18) should be compared with equations (4)-(7) of Mardling (2007) for the coplanar case, noting that here the inclination functions are such that fn(0) = 1, n = 1, 2, 3 and g2(0) = 0.
Our aim now is to write down equations governing the dominant long-term behaviour of the two-planet system under the action of tidal dissipation, spin-orbit coupling and the relativistic potential of the star. In particular, we wish to study the behaviour of the system in the e b − η plane in order to determine whether a relaxed system with non-zero mutual inclination is able to tell us anything about the internal structure of planet b. Noting from Figure 3 (and the following analysis) thaṫ ω b =̟ b −Ω b is approximately constant, and that as long as η librates, the argument ζ ≡ ̟ b + ̟c − 2Ω b = 2ω b − η so thaṫ ζ ≃ 2ω b , the equations governing e b and η are approximatelẏ witḣ and we have taken the time origin to coincide with ω b = 0. Here WT = τ −1 circ is the inverse of the tidal circularization timescale, and both Wq and Wo have analogues in the coplanar theory 2 and are such that and equal to minus the precession frequency divided by cos i b (see equation (A15)). Parameters which are zero for coplanar systems are while does not appear explicitly in the coplanar theory. Note that since a b /ac, i b , ic and ec are all approximately constant over a few times the circularization timescale (at least for small e b ), each of WT , Wq, Wo, WΩ,ω b , A and B are also approximately constant on this timescale. Referring to Figure 2 showing the relaxed behaviour of e b for different i b , our aim now is to characterize the limit cycle demonstrated in Figure 1 by determining the average value of e b as well as its amplitude as a function of the mutual inclination. Equations (9) and (10) together represent a damped nonlinear non-homogeneous system of ordinary differential equations with periodic coefficients for which the existence of periodic (limit-cycle) solutions are suggested by Figures 2 and 3 (see, for example, Jordan & Smith (1999) for a discussion of such systems). Assuming a limit-cycle frequency of 2ω b and zero phase, 3 that η librates rather than circulates (this is true for i b < ∼ 33 o for the HAT-P-13 system as determined from numerical integrations), and that relaxed values for e b and η are π/2 out of phase (a reasonable assumption given the form of (9) and (10)), the amplitude of variation of e b and η on the limit cycle can be estimated by putting and where e (av) b and ηav are defined in the following procedure. Substituting these expressions into (9) and (10), assuming A lc ≪ 1 and η (lc) ≪ 1, matching constant terms and those with phase 2ωt and neglecting terms with phase 4ωt, we obtain where ∆ is defined in (12), together with which reduces to (1) when i b = 0, and Expressions (19) and (20) should be compared with equations (36) and (46)  (1 ± A lc ) .
In fact the true amplitude tends to be around 2 3 A lc , perhaps due to the neglect of terms proportional to cos(4ω b ) and sin(4ω b ) (these are associated with the shape of the true limit cycle -see Figure 1). We therefore replace A lc with in the analysis and figures that follow (as well as Figure 1). Figure 4 (21) , the range of values of i b corresponding to libration of η around η = ηav. Also plotted are maximum, minimum and average values of e b from a series of numerical integrations for HAT-P-13-like systems with 10 o i b 50 o (red solid and dot-dashed curves).
corresponds to A * lc = 1 (compare with the discussion following equation (23) in Mardling (2007)). For i η again librates, this time around η = π + ηav. Here i K b corresponds to the minimum value of i b for which Kozai oscillations occur for a given value of γ (see Section 3.2).
While it is straightforward to determine e min b and e max b for all possible (non-relaxed) coplanar configurations and their associated fixed points (Mardling 2007), it not clear how to do this for inclined systems and their associated limit cycles for , that is, systems which η either circulates or librates around π. Here we content ourselves with trial and error expressions, obtained by trying different integer values of n1 and n2 in e (n1A lc + n2), and comparing these with numerical solutions. The results are summarized in Table 2  curves cross. Given the form of these expressions and the accuracy with which they fit the numerical data, it seems likely that they are generic. on the Love number of planet b. Note that for the case k b = 0 we have also set k * to zero (although in the k b = 0.3 case the quadrupole moment of the star contributes only 1.4% to the total value of γ, assuming the star spins with the same period as the Sun). We conclude that the analysis and conclusions  Table 2, while that for e (av) b is given by (19). Note that η librates around zero for 0 gives a fairly accurate estimate of k b for those inclinations.   Batygin, Bodenheimer & Laughlin (2009) are valid as long as the mutual inclination of planets b and c is less than around 10 o ; for higher values of i b , a measurement of e b does not unambiguously determine k b , although one can make arguments about the likelyhood of a system being near the top or bottom of the modulation cycle of e b . While a Rossiter-McLaughlin estimate of the stellar obliquity will not help to constrain i b (see discussion in Section 4.1), it should be possible to fit combined transit and spectroscopic data for the mutual inclination (Nesvorný 2009). We now argue that for realistic values of Q b , the HAT-P-13 system cannot have a mutual inclination in the range 54 o to 126 o . (all other parameters are the same as for the cases presented in Section 3). The system exhibits Kozai oscillations, detail of which is shown in panels (b) and (e) for e b and i b respectively. However, given that there is no commensurate variation of a b , these come at the price of severe tidal decay of the orbit of planet b as shown in panel (c). An upper bound for the timescale for decay of the semimajor axis of planet b's orbit due to planetary tides is τa = τcirc/e 2 b , where τcirc is given by (5). With e b at least an order of magnitude higher than for the cases considered in Section 3, here orbital decay occurs more than one hundred times faster during the high-eccentricity phase for a given value of Q b (note that for this example we have used a Q-value 10 times as large to clearly demonstrate Kozai cycles). The effect of this strong tidal damping is to quickly reduce the mutual inclination to a value for which Kozai oscillations no longer occur (around 53 o ; see discussion below) and the apsidal lines become locked with η ≃ 25π (panel (f)). The start of the slower phase corresponds to a value of planet b's semimajor axis of around 0.041 AU (the observed value is 0.0426 AU). Thereafter there is a slow decline in i b , a b and e b . 25π Figure 5. The Kozai regime. Here the initial mutual inclination and eccentricity of planet b are 60 o and zero respectively, and Q b = 100. The eccentricity of planet b is forced to high values initially (panel (b)) and there are corresponding variations in i b (panel (e)), however, this extreme behaviour comes at the price of severe tidal decay of the orbit of planet b (panel (c)) and the mutual inclination (panel (d)) until the system is no longer capable of driving Kozai oscillations. This occurs at around i b = 53 o at which point the system becomes trapped on a limit cycle for which η = 25π + ηav and e av b ≃ 0.02. This value of i Table 2) is consistent with the analysis summarized in Figure 6. Note that mc = m min c / cos 60 o .  (2007) give an expression for the maximum eccentricity, e K b , attained during a Kozai cycle in the presence of other perturbing forces such as spin, tides and relativity. Their equation (34) (which holds in the case that the initial inner eccentricity is zero) can be rearranged so that

Fabrycky & Tremaine
where x is the minimum positive root of the cubic with C1 = 5 3 cos 2 i b (0) and C2 = 2 9 γ, i b (0) being the initial inclination of planet b (which equals the initial relative inclination since ic(0) = 0). Note that this expression does not include any contribution from octopole terms in the disturbing function, nor does it include terms of order e 2 b or higher in the tidal distortion. Figure 6(a) shows the dependence of e K b on i b (0) (solid curve), taking into account the fact that the minimum mass of planet c is scaled by (cos i b ) −1 so that the apsidal advance of planet b is completely dominated by planet c when ic = 90 o (see also Figure 3 of Fabrycky & Tremaine (2007)). The circles correspond to "true" maximum values of e b for this case, calculated using the averaged code (with no dissipation) so that all tidal and octopole terms are included. While the maximum and minimum values of e b do not vary from one Kozai cycle to the next when only point-mass quadrupole terms are included in the equations of motion, this is not true when extra accelerations are included, although the variation is generally small (less than 10%). We therefore integrated for several modulation cycles and recorded the maximum value of e b during that time. The results suggest that the terms not included in (23) make a significant difference for high inclinations. Also shown in panel (a) is e K b when γ = 0 (dashed curve). Figure 6(b) shows the dependence of the orbit decay timescale on mutual inclination, for which an estimate is given by That this gives a good estimate (rather than using, say, the average value of e b during a Kozai cycle) is supported by the example in Figure 5 for which panel (c) provides an estimate over the first 10 4 yr of τa = 1.25×10 5 yr while (25) with e K b = 0.4 (panel (b)) gives τa = 2.5 × 10 5 yr, that is, (25) represents an upper bound. The circles in panel (b) of Figure 6 correspond to those in panel (a) and show that (23) may be used in (25) to estimate τa in spite of the fact that not all relevant terms are included. We may therefore conclude from panel (b) that given an estimated age of around 5 Gyr, the mutual inclination of the HAT-P-13 system cannot be between 54 and 126 deg for values of Q b less than 10 6 . Note that if Q b is as high as 10 7 , the system will not yet have relaxed to the limit-cycle state.
Finally, in order to be able to place bounds on the mutual inclination of any given observed system, it is important to know the range of inclinations for which Kozai oscillations of e b occur given a particular value of γ. Figure 6(c) shows that Kozai oscillations are completely suppressed for γ > ∼ 9 for any inclinations. This also has repercussions for the maximum possible stellar obliquity of a system, discussed in Section 6.

Limit-cycle behaviour of retrograde systems
While relaxed prograde systems are characterized by the libration of the angle η = ̟ b − ̟c, relaxed retrograde systems are characterized by libration of the angle ζ = ̟ b + ̟c − 2Ω b . This can be understood as follows.
Noting that η = 2ω b − ζ, and thatη ≃ 2ω b when ζ librates, equations (9)-(11) can be writteṅ witḣ where and with WΩ defined in (14). Parameters which are zero for retrograde coplanar systems are whilė Note that the rate of precession of the node of planet b about the invariable plane normal is negative for prograde orbits and positive for retrograde orbits (equation (A15)). Given the symmetry properties of the various inclination functions (listed in the paragraph following (A9)), we see that the relaxation behaviour of those retrograde systems for which ζ librates is identical to that of prograde systems for which η librates (albeit with a slightly different limit cycle frequency). In particular, the theory developed in Sections 3.1 and 3.5 for prograde systems caries directly over to retrograde systems if one replaces f2

The relationship between e (av) b
, the equilibrium radius of planet b, its Q-value and its core mass The above analysis implicitly assumes that the value taken for R b is its equilibrium value. In fact, for a given set of orbital parameters and planetary core mass, the equilibrium value of R b , R (eq) b , depends on Q b . This, in turn, determines e (av) b for a given value of i b . Moreover, the Love number depends on the core mass as well as the radius. We can quantify these relationships using the data provided in Table 1 of Batygin, Bodenheimer & Laughlin (2009) as follows.
For each value of the core mass of planet b, mcore, Batygin, Bodenheimer & Laughlin (2009) provide three sets of values of (R (eq) b /RJ ,Ė tide ), whereĖ tide is the tidal power required to maintain the radius at R (eq) b given the cooling rate −L b for such a structure. 4 The three values of R (eq) b /RJ correspond to the best-fit observed value ±1σ, these being (1.20,1.28,1.36). Desirable properties of the relationship between L b =Ė tide and R (eq) b /RJ are limR b →R J L b = 0 and d 2 (log 10 L b )/d(log R b ) 2 < 0 (see Figure 3 of Bodenheimer, Lin & Mardling (2001)). A function with these properties is with least-squares parameters bn and cn depending on mcore, as well as the choice for L b (R (eq) b /RJ = 1.01)/L⊙ ≡ λn (see Figure 7(a)). These are listed in Table 3. Also listed are least-squares exponential-fit parameters for the Love number for each core mass such that Fitting laws (33) and (34) are plotted in Figure 7. Following Mardling & Lin (2002) Section 4, the rate of change of the radius of planet b is given bẏ 4 Note that stellar insolation is also included in their calculations so our values for e (av) b below are slightly overestimated.
where α b is the moment of inertia coefficient of planet b 5 and E b = − 1 2 Gm * m b /a b is the orbital binding energy, while to O(e 2 b ), O(Ω 2 e /n 2 b ) and O(Ω 2 q /n 2 b ), the orbit-averaged expression for E tide is given by 6 Here Ωe = Ω b ·ê b and Ωq = Ω b ·q b , where Ω b is the spin vector of planet b andq b =ĥ b ×ê b , withê b andĥ b unit vectors in the direction of periastron and planet b's orbital angular momentum per unit mass respectively. Note that (36) assumes Ω b ·ĥ b = n b , that is, planet b is synchronously rotating with its orbit. The quantities Ωe and Ωq can remain non-zero for much longer than the tidal damping timescale if Ω * × h b = 0 and/or h b × hc = 0, where Ω * is the star's spin vector and hc is the orbital angular momentum per unit mass of planet c. The terms involving Ωe and Ωq in (36) are associated with the so-called obliquity tide (Fabrycky, Johnson & Goodman 2007); for the HAT-P-13 system they contribute less than 1% to Ė tide for all relative inclinations studied numerically, and will from hereon be ignored. The planet's radius will increase or decrease according to whetherĖ tide is greater or less than L b , at a rate enhanced by the factor (m * /m b )(R b /a b ) which is 21.5 for the HAT-P-13 system. Once the system achieves equilibrium (ie, onceĖ tide = L b ), the radius of planet b will shrink on the timescale τa, that is, the system will remain in the equilibrium state appropriate for the current value of a b . An example illustrating this behaviour is shown in panel (i) of Figure 10. For fixed Q b and mcore, (33) and (36) may be equated to give a relationship between e b and R (eq) b (using (34) to express k b in terms of R b ). Equation (19) provides a second relationship between these two variables for fixed i b , and combining the two gives one between R  Table 2. Given τe ∝ τa ∝ e −2 b (Section 2), we were able to confirm numerically that using e (av) b rather than e max b gives a good estimate for the decay timescale of e b . This is in contrast to systems with high eccentricities (see Section 3.2). Figure 8(c) presents an alternative view for low inclinations, with each set of curves corresponding to a different core mass. Given the lower bound placed on Q b by τe (Section 2), consistent with the fact that higher rather than lower values of Q b are suggested by the orbital parameters of short-period planets (Wu 2003 (v) Since 0 sin ic 1, cos i bc cos ∆Ω 1 from (iv). But since cos i bc = cos 50 o = 0.64, this contradicts (ii). Thus a mutual inclination of 50 o is inconsistent with observations. The argument against coplanar retrograde systems follows similarly, with cos i bc ≃ −1 and the librating angle ζ = ω b + ωc − ∆Ω = 0 so that ∆Ω = −48 -44 o with cos ∆Ω > 0, making sin ic < 0. Figure 8 also suggests that lower rather than higher core masses are favoured. Note that for inclined systems for which , τe can be estimated by replacing ∆0 with ∆ in (8), where ∆ is given by (12). Thus since ∆ ∆0, the lower bound for Q b given by τe for coplanar systems represents a lower bound for higher inclinations in this range.
More accurate measurements of e b and R b will allow refinement of the statements above. Note that a small value of i b does not imply that a Rossiter-McLaughlin measurement of the stellar obliquity will result in a small value; its maximum depends not only on the mutual inclination of the two orbits, but also on the stellar obliquity relative to orbit c (see Section 4). , the equilibrium radius of planet b, its Q-value and its core mass for the HAT-P-13 system (note different scales on panel (c)).

Timescale for decay of the mutual inclination in a relaxed system
Taking i b as a proxy for the mutual inclination (since we are assuming most of the angular momentum of the system resides in the outer orbit), we can use (A13) to estimate a timescale for the decay of the mutual inclination of a relaxed system. Figure 9(a) shows 10 5 years of evolution of i b for the systems shown in Figure 2, for which Q b = 10 and k b = 0.3. In order to clearly demonstrate long-term trends, we have plotted the quantity i b − i bc (0), where i b is measured relative to the invariable plane and i bc (0)  and η (lc) respectively (equation (17)), and recalling that where with ∆ defined in (12), i * b is the first root of Ψ and τi is time it takes for the system to relax to i b = i * b . For i b = 30 o , τi ≃ 80 Gyr. The function Ψ(i b ; ∆, γε 3 c ) is plotted in Figure 9(b) as a function of i b for HAT-P-13 system parameters with k b = 0.3. We see that Ψ > 0 for i b < i * b ≃ 24 o , that is, i b actually increases for this range of inclinations, while for greater inclinations i b decreases. The red dashed lines in Figure 9(a) have slopes equal to di b /dt , and are in good agreement with the general trend of the numerical solutions. Also shown are numerical solutions for i b = 40 o and 50 o . The slope of the trend for i b = 40 o is 0.4(10 5 /Q b ) o Gyr −1 , while that for i b = 50 o is a factor of three higher. Thus the timescale for the decay of the mutual inclination of HAT-P-13-like systems for which i b < 50 o is considerably longer than the age of the system for reasonable values of Q b .
A similar analysis can be done for retrograde systems, and we can conclude generally that inclined systems cannot relax to the coplanar prograde or retrograde state as long as e (av) b > 0, and that they relax to a mutual inclination given by one of the roots of Ψ(i b ) = 0 (or Ψ(π − i b ) = 0 for retrograde) as long as τi < τc and τi < τa. Note that the mutual inclination changes by the same amount as i b , and that (37) does not apply to systems with mutual inclinations greater than i (c 1 ) b for which it represents a lower bound.

A SCENARIO FOR THE ORIGIN OF HAT-P-13-LIKE SYSTEMS
The high eccentricity of HAT-P-13c suggests either a violent scattering history, or Kozai-type interactions with another planet or star, or that it was formed through gravitational collapse; it seems unlikely that such a high eccentricity could result from single planet-disk interactions alone (see, for example, Artymowicz (1993)). Kozai forcing would also affect planet b, and it is unclear without a detailed study whether or not a suitable configuration exists which would not cause the rapid decay of its orbit.
Here we focus on the scattering scenario, introducing a third planet ("planet d") which is ultimately ejected from the system. The mass ratio of planet b to planet c is too low for c to have attained its high eccentricity during a scattering event with b. Three models are presented, all of which have identical initial conditions except for the eccentricity of planet d, e d , which for model 2 differs from 1 by two parts in 10 5 while for model 3 it differs by 0.05. The outcomes are significantly different, with model 1 producing a value for ec almost identical to that for HAT-P-13c, while models 2 and 3 produce lower values at 0.46 and 0.36. Of particular interest is the relative inclination of the orbits of b and c following the escape of d, and Table 4. Data for models 1, 2 and 3 at t = 0 and around time planet b first achieves equilibrium Equilibrium has not yet been established.
the accompanying stellar obliquity relative to planet b, both of which differ significantly from model to model, as does the equilibrium value of R b . The initial configuration data are listed in Table 4, together with those at the time planet b first achieves an equilibrium radius, that is an equilibrium between the rate at which tidal energy is injected and the rate at which the planet can cool. Data of particular significance are highlighted in bold. Inclinations with superscripts (i) and (f ) are measured with respect to the original and final invariable planes respectively, the former including planet d and the latter not. The quantity θIP is the angle between the invariable planes before and after the escape of planet d. Again we use the averaged code of Mardling & Lin (2002), this time without suppressing the evolution of the radius of planet b which is evolved according to the scheme discussed here in Section 3.4 with mcore = 80M⊕. Note that this code involves averaging over the innermost orbit only; the two outer orbits are integrated directly and are therefore susceptible to instability through the overlap of mean-motion resonances (Mardling 2009a). While the mass of planet b is taken as that of HAT-P-13b, the mass of planet c is taken as the minimum mass of HAT-P-13c, that is, it is not scaled by cos i b . The mass of planet d is taken to be similar but less than that of planet c, ensuring that the probability that planet c is left with a high eccentricity is significant. Initial values of the eccentricities and inclinations are chosen to be consistent with formation and subsequent migration in a protoplanetary disk, and the period ratio (which is similar to that for Jupiter and Saturn) puts the system just inside the stability boundary of a system with these masses and eccentricities (see, for example, Figure 9(a) of Mardling (2009a)), consistent with having just emerged from the protection of the disk. 7 Moreover, the initial values of ac and the semimajor axis of planet d, a d , are such that the energy needed for escape of planet d and provided by the orbit of planet c reduces ac to a value similar to that of HAT-P-13c. The initial radius of planet b is 1.5RJ , consistent with the radius of a young planet recently arrived at its present location (see Figure 1 of Bodenheimer, Lin & Mardling (2001); disk lifetimes suggest that such a planet would probably have a higher radius). The initial mean longitudes and orientation angles of planets b and c are specified relative to the outer orbit; these are λ bd = 0 o , i bd = 5 o , ω bd = 0 o , Ω bd = 0 o and λ cd = 180 o , i cd = 10 o , ω cd = 0 o , Ω cd = 180 o respectively and are such that longitudes are measured with respect to the periastron direction of orbit d.
The time evolution of various quantities is plotted in Figures 10, 11 and 12, with detailed descriptions and discussion provided in the captions. The choice of an initial period ratio for orbits c and d of around 5:2 affects possible outcomes in the following ways (given that we wish to approximately reproduce the HAT-P-13 system). In view of the fact that we wish to place the system near the edge of the stability boundary, the choice of e d (0) is restricted by the fact that a value significantly higher would render the system far from the boundary and hence violently unstable. Our choice of m d affects the choice of ac(0); a smaller value of m d requires a smaller value of ac(0), however, it tends to produce smaller values of ec following escape of d and is less likely to produce high stellar obliquities (see model 3). On the other hand, a heavier companion is more likely to eject planet c (an outcome also possible for m d = 12MJ ; here we restrict ourselves to systems in which the outer planet is ejected). A general detailed study is required to quantify possible outcomes of such a scattering scenario, especially in the light of the recent discovery of two retrograde systems Anderson et al. 2010). For now our aim is to demonstrate its potential to produce a range of eccentricities, mutual inclinations, stellar obliquities and planetary radii.
Of particular interest is the stellar obliquity relative to the orbit of planet b, ψ * b , a quantity which can be measured directly (at least in sky projection; see Fabrycky & Winn (2009) for a discussion of the statistical properties of this quantity). In the following section we outline the mechanics of stellar obliquity in a two-planet system.

Stellar obliquity in a two-planet system
The stellar obliquity relative to each planetary orbit in a system reflects conditions at the time of its formation. In the scattering models presented in this section we assume zero stellar obliquity relative to planet b initially, while the stellar obliquity relative to planet c is 15 o . Table 4 lists the ranges for ψ * b following the escape of planet d, demonstrating that significantly different outcomes are possible from models with very similar initial conditions. The variation in ψ * b depends on two quantities: the angle between the stellar spin axis and the normal to the invariable plane, θ * , and the angle between planet b's orbit normal and the normal to the invariable plane, i b . When ψ * b = 0, the variable torque on the spin bulge of the star from planet b results in nutation of the star's spin axis, that is, a variation of θ * . This can be quite significant; in model 3 it is 7 o compared with 1 o for model 1 (see Table 4). 8 Its average value, however, depends on the history of the system; if no planets have been ejected since the system's formation, θ * is likely to be modest, while the escape of one or more planets can, depending on how much angular momentum the planets carry away, result in an invariable plane normal which points in a significantly different direction to the original (θIP in Table 4), thereby affecting θ * . Barker & Ogilvie (2009) have shown that the decay timescale for ψ * b is around τa, the timescale for the decay of the orbit, so very little reduction in the average value of θ * is expected over the lifetime of the orbit of HAT-P-13b.
Regarding the relevant angles as spherical polar angles, ψ * b is given in terms of θ * , the spin axis node angle, ϕ * , as well as i b and Ω b , by Since the variations in θ * and i b are small (but see next section), ψ * b cycles approximately between |i b + θ * | and |i b − θ * | as ϕ * − Ω b cycles between 0 and 2π over b's precession cycle, 9 or, since the orbit of planet c contains most of the angular momentum of the system so that θ * ≃ ψ * c, the variation is approximately |i b ± ψ * c|. In fact, since the nutation period of the star is equal to the precession period of planet b, extrema of ψ * b coincide with extrema of ψ * c (see panel (e) of Figures 10 to  12). A significant difference between the first two models and model 3 is that before the escape of planet d, the latter suffers a significant transfer of angular momentum between orbits c and b during a period of high eccentricity of planet c, resulting Figure 10. A chaotic origin for the HAT-P-13 system (note different timescales for each panel). At time t = 0 three planets are present, with migration leaving planet b 0.043 AU from the star with a radius of 1.5R J and a core of mass 80M ⊕ , and planets c and d with period ratio 2.46. While the presence of an outer disk has previously limited the variation of the eccentricities thereby protecting the system against instability, by t = 0 the disk surface density has reduced sufficiently to allow the system to become unstable. Panel (a) shows all three eccentricities, each remaining moderate for the first 6000 years. After 16,000 years planet d escapes the system. Following escape, the orbit of planet b precesses around the new invariable plane whose normal is approximately parallel to that of planet c (since with d removed from the system it now contains 99% of the total angular momentum; compare panels (b) and (c)). The inclinations of the three orbits to the original invariable plane are modest initially, and they remain so during and after the scattering process. Following the escape of planet d, the stellar obliquity relative to b's orbit (panel (e), black solid curve) oscillates about a mean equal to i b (blue dashed curve) and with an amplitude equal to ψ * c, the stellar obliquity relative to c's orbit (red dashed curves; see also panel (f) and Figure 13). Panel (g) shows the early behaviour of e b ; before the scatter planet b is effectively decoupled from the rest of the system and its eccentricity monotonically decreases on the tidal circularization timescale, while after the scatter it is governed by planet c. An artificially low Q-value of 40 is used for planet b, and its radius and Love number are evolved according to (33) and (34) respectively (k b (0) = 0.13 for R b (0) = 1.5R J and mcore = 80M ⊕ ). Panels (h) and (i) show the long-term evolution of e b and R b respectively. The system evolves to a limit cycle after about 0.12 Myr, with e (av) b decreasing as R b increases (given the dependence of e (see Figure 8). During the 3.5 Myr integration, the semimajor axis decreases from 0.04300 to 0.04055 AU, giving an orbital decay timescale τa = a b /ȧ b of 158(Q b /10 5 ) Gyr. Also shown in panel (i) are curves R b (t) = R b (0)exp(t/τ + ) (red dashed curve) and R b (t) = R b (t * )exp[−(t − t * )/τa] (blue dotted curve), where τ + = (m b /m * )(a b /R b )τa and t * = 3.5 Myr. These curves demonstrate that initially the rate of change of the radius is dominated by tidal heating, that is, the cooling term contributes very little (see equation (35)), and once equilibrium is reached it is continually reestablished as a b decreases on the timescale τa.
in a particularly high maximum value of ψ * b . The importance of this difference can be understood as follows. In each model, the (average) value of ψ * c depends on its value before the scatter, as well as the change in the inclination of c's orbit relative to the original invariable plane. Since the orbit of c effectively coincides with the invariable plane once d has left the system, i b ≃ ψ * c for a system in which there is no close encounter between b and c (at least for the cases considered here for which ψ * b (0) = 0), and the variation in ψ * b is approximately equal to 2ψ * c. On the other hand, i b is significantly greater than ψ * c in model 3 because of the close encounter of b and c, and as such the maximum value of ψ * b is high. Figure 13 illustrates the mechanics of stellar obliquity in a two-planet system for the cases (a) ψ * c = 20 o and i b = 70 o (a configuration ruled out for the HAT-P-13 system because 54 < i b < 126 o ; see Section 3.2) and (b) ψ * c = 70 o and i b = 20 o , Figure 11. A different outcome for the HAT-P-13 system (note different timescales for each panel). Initial conditions are the same as for model 1 but with a difference in e d of 2 × 10 −5 . Following the scatter of planet d from the system, ec = 0.45 (panel (a)), resulting in a smaller value of e (av) b (panel (g)). This time the mutual inclination reaches 20 o , as does the stellar obliquity relative to c's orbit (panels (c) and (f)). This results in a maximum stellar obliquity relative to b's orbit of i b + ψ * c = 40 o and a minimum of 3 o > i b − ψ * c (panel (e); the variation of ψ * b is not a simple sinusoid when its mimimum is near zero, similar to the behaviour of eccentricity in the same circumstances). Note that the precession timescale is longer than that of model 1 by a factor (ε (2) c /ε (1) c ) 3 , where the superscripts refers to the model numbers (see equation (14)). The long-term evolution of e b (panel (g)) is a limit cycle for which e (av) b increases slightly for the first 10 Myr or so, corresponding to the relatively rapid (but still very slow) decrease of R b (panel (h)), then decreases once equilibrium is established. The theoretical prediction according to Section 3.4 gives R  0)) and t * = 20 Myr. These curves demonstrate that the rate of change of the radius is never dominated either by heating or cooling (see equation 35), being close to equilibrium initially. After around 10 Myr equilibrium is established and R b decreases on the timescale τa.
so that 50 o ψ * b 90 o for both. Also indicated is the original invariable plane in the case that the system contained a third planet which has since escaped, and whose angular momentum was such that the stellar spin axis was parallel to the invariable plane normal. The ramifications for the origin of retrograde systems are clear (although configuration (b) would probably require a close encounter with a passing star rather than the escape of a companion planet to cause such a dramatic change in c's orbit). Given the amount of angular momentum transferred from c to b depends on all the system parameters, it is easily conceivable that a scattering scenario similar to the one described here (either bound or flyby, perhaps involving exchange of c and d) could produce retrograde systems such as HAT-P-7 and WASP-17. It is interesting to note that it is possible for the relative inclination of a system to pass through 90 o without destroying planet b through Kozai oscillations as long as γ is large enough (see Section 3.2).
We end by considering the effect on the mutual inclination of non-zero stellar obliquity with respect to the invariable plane, which can be significant if θ * is significantly non-zero. large angular momentum transfer to b Figure 12. High stellar obliquity due to a close encounter between b and c (note different timescales for each panel). Initial conditions are the same as for model 1 but with e d = 0.1. This time the mutual inclination reaches 35 o , while the maximum stellar obliquity relative to c's orbit after the scatter is 19 o (panels (c) and (f)). This results in a maximum stellar obliquity relative to b's orbit of i b + ψ * c = 53 o and a minimum of i b − ψ * c = 23 o (panel (e)). Unlike models 1 and 2 where the stellar obliquity can be attributed almost entirely to angular momentum transfer between c and d (and the initial value of ψ * c), i b and ψ * c are significantly different, a result of a close encounter between planets b and c around 17,000 yr when ec reaches 0.87 (see panel (g) in which the periastron separation, pc, reaches a minimum of 4a b ). During this high-eccentricity phase, the torque from planet c tilts the orbit of planet b through 30 o , and while some of the associated angular momentum is subsequently returned to c's orbit, the mutual inclination remains high following the escape of planet d. Since the value of ec is relatively low at 0.36 after the escape (panel (a)), e (av) b is a mere 0.001 (panel (h)), with a similarly small amplitude in spite of the relatively high mutual inclination. As a consequence, tidal heating is weak and the predicted equilibrium radius of planet b is only 1.14R J (recall Q b = 40). Panel (i) shows the evolution of

Effect on the mutual inclination of non-zero stellar obliquity with respect to the invariable plane
The numerical results presented in Section 3 assume zero stellar obliquity relative to the invariable plane, that is, θ * = 0. Consequences of this are that the relaxed state illustrated in Figure 1 (red curves) and Figure 2(a) exhibits constant amplitude variations, and that the inclination i b is effectively constant. The latter is clearly not the case in the three models discussed above (panel (c) of Figures 10, 11 and 12). Figure 14 and zero respectively. The limit-cycle amplitude varies noticeably and the modulation period is slightly longer, with the behaviour persisting during a 10 6 year integration. This can be understood as follows. The torque on the orbit of planet b due to the star's spin oblateness is given by equation (48) in Mardling & Lin (2002), and this contributes to the rate of change of i b according to equation (29) of the same paper. Expressing the stellar spin vector Ω * and the basis vectorsê b ,q b andĥ b referred to in Mardling & Lin (2002) 10 in terms of the invariable plane reference basis via θ * , the star's node angle ϕ * , and the Euler angles ω b , Ω b and i b , the contribution to di b /dt from the star's spin oblateness becomes igure 13. Illustration of high stellar obliquity in cases where (a) planet c is scattered through a small angle during its interaction with a third planet, and planet b is scattered through a large angle during its interaction with planet c, and (b) planet c is scattered through a large angle during its interaction with a third planet (or passing star), and planet b is also scattered through a large angle during its interaction with planet c, but in such a way that it ends up with a small value of i b . Here we are assuming that the stellar spin direction is parallel to the invariable plane normal before interaction with the fourth body, and that the system was originally coplanar. While both configurations have the same range of values of ψ * b = |i b ± ψ * c|, configuration (a) is ruled out for HAT-P-13 because 54 < i b < 126 o (see Section 3.2), while configuration (b) cannot be ruled out (although it would probably require a close encounter with a passing star rather than the escape of a companion planet to cause such a dramatic change in c's orbit). Figure 14. The effect of non-zero stellar obliquity relative to the invariable plane on the relaxed state for (a) e b and (b) i b . The limit-cycle amplitude (ie, the amplitude of variation of e b ) varies regularly on a timescale of twice the limit-cycle period, the latter being 2π/2ω b , and its modulation period is slightly longer. The amplitude of variation of i b is significantly enhanced by an amount proportional to sin(2θ * ). For comparison, the variation of i b for the case θ * (0) = 0 is shown (red curve in panel (b)). The small variation is due to nutation of the stellar spin axis through about 2.5 o .
Taking θ * , i b and ϕ * to be approximately constant over an orbit precession cycle (the precession rate of the stellar spin axis is 40 times slower than the orbit precession rate for this example), and putting Ω b =Ω b t + Ω b (0) whereΩ b is given by (A15) and is also approximately constant, the variation in i b due to the star's spin oblateness, ∆i b , is obtained by integrating (40) over an orbit precession cycle to give ∆i b = γ spin * ε 3 c sin(2ψ * c) where γ spin * is the contribution to γ from the star's spin quadrupole moment and is given in Table 1 for HAT-P-13 for mc = m min c . Putting ψ * c = 50 o then gives ∆i b = 1.1 o , in good agreement with Figure 13(b). Such a variation in turn produces a variation in e max,min b of the order seen in panel (a) of the same figure.

SUMMARY
The ideas and results presented in this paper can be summarized as follows: 1. Generalizing the results of Mardling (2007) to non-coplanar systems for which most of the angular momentum of the system resides in the outer orbit, we find that under the action of tidal dissipation in planet b, the system evolves to a limit cycle rather than a fixed point in e b − η space, with the average value of e b , e (av) b , decreasing and the limit cycle amplitude increasing with increasing mutual inclination, and limi b →0 e (av) b = e (eq) b . For the HAT-P-13 system, limit cycle behaviour occurs for i b < ∼ 33 o (libration of η around 2nπ) and 46 o < ∼ i b < ∼ 54 o (libration of η around (2n + 1)π). No limit cycle exists for 33 o < ∼ i b < ∼ 46 o (η circulates).
2. For systems with 54 o < i b < 90 o (90 o < i b < 126 o ), Kozai oscillations coupled with tidal dissipation in planet b act to reduce (increase) the mutual inclination until i b < 54 o (i b > 126 o ) on a timescale much less than the age of the system. We conclude that the HAT-P-13 system cannot have a mutual inclination between 54 o and 126 o for Q b < ∼ 10 6 .
3. For retrograde systems, the limit-cycle behaviour is the mirror image of that for prograde systems, apart from a slightly different limit-cycle frequency (compare (32) with (16)). (2009)  , its Q-value and its core mass (Figure 8), and conclude that for Q-values greater than the lower bound imposed by the timescale for decay of e (av) b , the orbits of planets b and c are likely to be either near prograde coplanar, or have mutual inclinations between around 130 o and 135 o . Lower rather than higher core masses are favoured. More accurate measurements of e b and R b will allow refinement of these statements.

The analysis and conclusions of Batygin, Bodenheimer & Laughlin
6. Inclined systems cannot relax to the coplanar prograde or retrograde state as long as e (av) b > 0, and instead relax to a mutual inclination given by one of the roots of Ψ(i b ) = 0 or Ψ(π − i b ) = 0 (see equation (38)). This will occur as long as τi < τc and τi < τa, both true for the HAT-P-13 system, however, in this case, τi is much greater than the age of the system. 7. A viable formation scenario for the HAT-P-13 system is that it originally contained a third planet which was scattered out of the system when the protoplanetary disk density dropped below the critical level for stability. Such a scenario is capable of producing a high eccentricity for planet c as long as the mass of planet c is sufficiently high, and may also produce significant mutual inclination, stellar obliquity and inflated planetary radii. 8. A Rossiter-McLaughlin measurement of the sky-projected stellar obliquity relative to planet b, together with a measurement of the mutual inclination, will allow us to constrain the stellar obliquity relative to planet c and hence obtain knowledge about the formation history.
In conclusion, it is likely that many more HAT-P-13-like systems will be discovered in the future including systems in which the outer body is a binary star companion. Such systems will contribute a wealth of information not only about the internal structure of the short-period planet, but also about the formation history of the system.

Note added in proof:
Since this paper was submitted more refined data for the HAT-P-13 system have become available (Winn et al. 2010), in particular, a new estimate for the eccentricity of planet b of 0.0142 +0.0052 −0.0044 . Moreover, there is now evidence for a distant third body in the system, as well as a Rossiter-McLaughlin measurement of the sky-projected stellar obliquity, the latter strongly suggesting that the stellar spin and the orbit normal of planet b are aligned. The conclusions drawn here remain valid, in particular those regarding coplanarity or otherwise when the refined value of ω b − ωc is used.

ACKNOWLEDGMENTS
The author wishes to thank Dan Fabrycky for inspirational discussions and encouragement, for carefully reading the manuscript, and especially for his generosity in allowing her to include his elegant argument ruling out the 45 -50 o and retrograde coplanar configurations for the HAT-P-13 system. Thanks also go to Eric Ford for a very valuable comment. respectively, where x represents any one of e b , a b /ac or sin ic (for example, x 2 may represent e b sin ic, and cos ic = 1 + O(x 2 )).