Are falling planets spinning up their host stars?

Mon. Not. R. Astron. Soc. 000, 1–16 (2008) Printed 21 March 2011 (MN LATEX style file v2.2) Are falling planets spinning up their host stars? D. J. A. Brown1⋆ , A. Collier Cameron1 , C. Hall1 , L. Hebb1 , B. Smalley2 , 1 SUPA, School of Physics and Astronomy, University of St Andrews, North Haugh, St Andrews, Fife KY16 9SS, UK. Group, School of Chemistry and Physics, Keele University, Staffordshire, ST5 5BG, UK. arXiv:1103.3599v1 [astro-ph.EP] 18 Mar 2011 2 Astrophysics Accepted 0000 December 00. Received 0000 December 00; in original form 0000 October 00 ABSTRACT We investigate the effects of tidal interactions on the planetary orbits and stellar spin rates of the WASP-18 and WASP-19 planetary systems using a forward integration scheme. By fitting the resulting evolutionary tracks to the observed eccentricity, semimajor axis and stellar rotation rate, and to the stellar age derived from isochronal fitting, we are able to place constraints on the stellar and planetary reduced tidal quality factors, Q′s and Q′p . We find that for WASP-18, log(Q′s ) = 8.21+0.90 −0.52 and +0.305 log(Q′p ) = 7.77+1.54 , implying a system age of 0.579 Gyr. For WASP-19 we obtain −1.25 −0.250 +2.19 +1.86 ′ ′ values of log(Qs ) = 6.47−0.95 and log(Qp ) = 6.75−1.77 , suggesting a system age of +1.1073 1.60+2.84 −0.79 Gyr and a remaining lifetime of 0.0067−0.0061 Gyr. We investigate a range of evolutionary histories consistent with these results and the observed parameters for both systems, and find that the majority imply that the stars have been spun up through tidal interactions as the planets spiral towards their Roche limits. We examine a variety of evidence for WASP-19 A’s age, both for the value above and for a younger age consistent with gyrochronology, and conclude that the older estimate is more likely to be correct. This suggests that WASP-19 b might be in the final stages of the spiral-in process, although we are unable to rule out the possibility that it has a substantial remaining lifetime. Key words: stars: rotation – planetary systems 1 INTRODUCTION There are currently 538 officially confirmed extrasolar planets known to astronomy1 , a substantial fraction of which belong to the ‘hot Jupiters’ category of Jovian-type planets that orbit within 0.1 AU of their host star. Improved observational precision has, in recent years, led to the discovery of extreme examples of this class of planet with orbital semimajor axes of the order of 0.02 AU or less(see, for example, Sasselov (2003), Hellier et al. (2009), Hebb et al. (2010)). The Jovian masses and small orbital separations of these planets imply that they interact tidally with their host stars. Tidal interactions between stars and planets were first studied in the context of the Solar system (Goldreich & Soter 1966), and have also been extensively studied as they apply to binary systems (Hut 1980), but it is only comparatively recently that such studies have been extended to consider exoplanetary systems. Tidal interactions lead to long-term changes in the orbital parameters of a planet, specifically the ⋆ E-mail:[email protected] As of 2011 March 15. A full list of known exoplanets can be found at http://exoplanet.eu 1 eccentricity and the semi-major axis (see e.g. Mardling & Lin (2002); Jackson, Greenberg & Barnes (2008); Barker & Ogilvie (2009); Jackson, Barnes & Greenberg (2009)). The end result of this process is the spiral-in of the planet towards its host until it reaches the Roche limit, where it undergoes mass transfer through the L1 point (Gu, Lin & Bodenheimer 2003), becoming disrupted in the process. Previous studies of tidal evolution have, for the most part, focused on the effect of tidal interactions on these two orbital parameters whilst neglecting the evolution of other parameters that are involved. Dobbs-Dixon, Lin, & Mardling (2004) demonstrated that the effects of tidal interactions extend to the rotation rates of the two bodies, but until recently little work had been done to investigate this aspect of tidal theory as applied to exoplanetary systems. There is, however, mounting interest in this aspect of exoplanetary tides, with several studies investigating the large scale effect of the presence of exoplanets on stellar rotation (Pont 2009; Alves, Do Nascimento & de Medeiros 2010; Lanza 2010), and some work investigating tidal interactions specifically (Leconte et al. 2010). It is well known that the rotation rate of an isolated star declines with age through the action of magnetic brak- 2 D. J. A. Brown et al ing (Weber & Davis 1967), which leads asymptotically to a power-law decay of rotation rate with the inverse square root of stellar age (Skumanich 1972). Barnes (2007) presented a method for determining the ages of stars using their rotation periods and colours, based on a Skumanich-type magnetic braking law and calibrated using the solar rotation period, as an alternative to stellar model fitting. What appears to be less widely recognized, or at least acknowledged, is that the torques that act on both star and planet as a result of tidal interactions can affect this natural rotational evolution, and may be sufficiently strong to overwhelm it, at least for short periods of time. In this paper we simulate the evolution of hot Jupiter systems to investigate the effects on the orbital eccentricity, orbital separation, planetary rotation rate and stellar rotation rate of tidal interactions between hot Jupiters and their host stars. In Section 2 we set out the mathematical basis for our models, whilst in Section 3 we describe the computational methods that were used. Section 4 and Section 5 contain discussion of the application of our simulations to the specific cases of the WASP-18 and WASP-19 systems respectively. It should be noted that these systems were not selected at random, but were chosen owing to pre-existing disagreement between age estimates calculated using different methods. These transiting hot Jupiters also have two of the shortest known orbital periods, and are therefore particularly susceptible to the influence of tidal interactions. The stellar and planetary moments of inertia are determined by their effective squared radii of gyration, αs and αp . Angular momentum is carried away from the system via a magnetically-channeled, thermally-driven stellar wind, at a rate described by a standard Weber-Davis model L̇wind = −Is κΩs Min(Ωs , Ω̃)2 (5) where the stellar moment of inertia Is = αs Ms Rs2 . The physical scaling of the braking rate is determined by the constant of proportionality κ, and Ω̃ is the ‘saturation’ rotation rate above which the stellar magnetic field strength is assumed to become independent of the stellar rotation rate. We use the expressions of Dobbs-Dixon, Lin, & Mardling (2004) to describe the tidal spin evolution of the planet to the tidal torque:   2  3   Rp n Ms Ωp 9 (6) f (e) − f (e) Ω̇p = 3 4 2 αp Q′p Mp a n and the star under the influence of both tidal and wind torques: 2  3   2   9 Mp Ωs Rs n Ω̇s = f (e) − f (e) 3 4 2 αs Q′s Ms a n −κΩs Min(Ωs , Ω̃)2 . The polynomials f3 (e) and f4 (e) have the form   5 6 15 2 45 4 e + e + e /(1 − e2 )6 f3 (e) = 1 + 2 8 16 (7) (8) and 2 TIDAL AND WIND EVOLUTION Following Eggleton, Kiseleva, & Hut (1998), Mardling & Lin (2002) and Dobbs-Dixon, Lin, & Mardling (2004) tidal energy is dissipated within a star and planet whose spin axes are aligned with the orbital axis at rates defined by the tidal quality factors Q′s = 3Qs /2k∗ and Q′p = 3Qp /2kp , where k∗ and kp are the tidal Love numbers of the two bodies. The eccentricity of the orbit evolves at a rate   5  81 n Ms Rp 11 Ωp ė = f (e) + −f 1(e) + 2 e 2 Q′p Mp a 18 n  5   81 n Mp Rs 11 Ωs + −f 1(e) + f2 (e) (1) 2 Q′s Ms a 18 n where f1 (e) =  1+  5 6 15 2 15 4 e + e + e /(1 − e2 )13/2 4 8 64 (2)  1+  3 2 1 4 e + e /(1 − e2 )5 . 2 8 (3) and f2 (e) = The star and planet have masses Ms and Mp and radii Rs and Rp respectively, and rotation rates Ωs and Ωp . The orbital frequency is defined by Kepler’s 3rd law, n2 = G(Ms + Mp )/a3 . If the system is fully aligned, then the total angular momentum of the orbit and the axial rotation of the star and the planet, perpendicular to the orbit, is s Ga(1 − e2 ) Ltot = Mp Ms + αs Ms Rs2 Ωs + αp Mp Rp2 Ωp . (4) Mp + Ms f4 (e) =  1 + 3e2 +  3 4 e /(1 − e2 )9/2 . 8 (9) The total angular momentum of the system evolves as s # " ȧ(1 − e2 ) − 2aeė G p + L̇wind = Mp Ms Mp + Ms 2 a(1 − e2 ) +αs Ms Rs2 Ω̇s + αp Mp Rp2 Ω̇p . (10) Dividing the angular momentum loss rate by the orbital angular momentum s Ga(1 − e2 ) (11) Lorb = Mp Ms Mp + Ms and rearranging, we obtain the expression for the evolution of the orbital semi-major axis:   Is κΩs Min(Ωs , Ω̃)2 ȧ ė e2 Is Ω̇s Ip Ω̇p =2 − − − .(12) a e 1 − e2 Lorb Lorb Lorb By integrating this equation we obtain an estimate of the time remaining to spiral-in for a planet orbiting a slowlyrotating star,  5 a 2Q′s Ms , (13) tremain = 117n Mp Rs giving a result almost identical to the estimate of the same quantity derived from a slightly different formulation by Levrard, Winisdoerffer, & Chabrier (2009). The quantity a/Rs is derived directly from the transit duration for a transiting planet, enabling tremain to be estimated from directlyobserved quantities for a given value of Q′s . Are falling planets spinning up their host stars? Leconte et al. (2010) recently suggested that parametrizing the tidal evolution equations in this manner, with the stellar and planetary tidal quality factors constant in time, is not equivalent to the more traditional approach of considering either a constant phase lag or a constant time lag of the tidal bulge. They further suggested that such equations, if truncated to O(e2 ), produce both qualitatively and quantitatively incorrect evolutionary histories for systems with e > 0.2, and promote the use of equations derived from Hut (1981). Despite using a constant tidal quality factor model the tidal equations that we use in this work are not truncated in this fashion, and thus do not suffer from the problems that Leconte et al. describe. 3 COMPUTATIONAL METHOD An estimate of the value of the scaling constant for magnetic braking, κ, is required to implement (5). Following the standard Weber-Davis model, the rate of change of the rotation rate of an isolated star owing to magnetic braking can be described by Ω̇s = −κΩ3s . (14) 3 Starting from a given set of initial conditions (t0 , e0 , Ωs,0 , Ωp,0 , a0 ), we integrate the four equations (1), (6), (7), and (12) using a fourth-order Runge-Kutta scheme, adapted from algorithms in Press et al. (1992). The unit of time is 1 Gyr, and the integration is allowed to run for the approximate main-sequence lifetime, tMS ≃ 7(Ms /M⊙ )−3 Gyr for a star of mass Ms , or until the planet reaches the Roche limit, as defined by Eggleton (1983). To calculate the stellar radius of gyration the metallicity and mass of the star are used to select an appropriate table from Claret (2004, 2005, 2006, 2007); at each timestep we interpolate through this data using the current time to derive a value for αs . This is not an entirely accurate method, as the series of papers by Claret provides only tables for discrete values of metallicity and stellar mass that in many cases do not coincide with the metallicity of the systems being studied. In these cases we therefore use the table that most closely corresponds to our data. We use the observed values of e, Ωs and a, and their uncertainties, to evaluate the goodness-of-fit statistic at each timestep according to χ2 (t) = (eobs − e(t))2 (Ωs,obs − Ωs (t))2 (aobs − a(t))2 + .(18) + 2 σe2 σΩ σa2 s Integrating this under the assumption that Ω2 << Ω1 and t2 >> t1 gives To this we add a Bayesian prior on the stellar age t to obtain the statistic Ω−2 2 , (15) 2t2 which leads to the standard power law for magnetic braking, C = χ2 (t) + κ= P ∝ tα , (16) with α = 0.5 (Skumanich 1972). This neglect of the initial conditions is justified by the observed strong convergence of spin rates to a narrow period-colour relation in the Hyades (Radick et al. 1987), and in other clusters of similar ages. We set Ω2 = Ωs,Hyades , calculated by scaling the Prot -(JK) colour relationship found for the Coma-Berenices cluster by Collier Cameron et al. (2009) and using J and K magnitudes taken from the SIMBAD on-line data archive, to obtain an estimate for κ at the age of the Hyades. This value is assumed to be constant throughout the evolution of the system. For simplicity, in all cases the planet’s rotation is initially assumed to be tidally locked such that Ωp,0 = n, but Ωp is permitted to evolve independently thereafter. More recent work has found that the braking law exponent, α, diverges slightly from the ideal Skumanich value in some environments (see e.g. Collier Cameron et al. (2009)). If we generalise equation 14 to raise Ω to the power β, then we derive 1 κ= Ω− α , t/α (17) a universal expression for κ with β = α1 + 1. We calibrated our implementation of the magnetic braking power law using the stars listed in table 4 of Collier Cameron et al. (2009), finding that an exponent of α = 0.495 ± 0.002 gives good agreement between the measured rotation periods and the evolved stellar rotation rates calculated using our method. We found that there was little to differentiate the evolutionary tracks produced using values of α within this range, and therefore adopted the central value of α = 0.495. (tobs − t)2 , σt2 (19) where tobs and σt are the stellar age estimated from the star’s density and effective temperature (Sozzetti et al. 2007), and its associated uncertainty. The age of the system is taken to be the time of the step at which C is a minimum. This forward integration method was built into two separate computational schemes: a grid search that carried out the integration for each node in a four-dimensional grid in a0 -e0 -log(Q′s )-log(Q′p ) parameter space, and a Markov-chain Monte-Carlo (MCMC) optimisation scheme to determine the posterior probability distributions of log(Q′s ), log(Q′p ), e0 and a0 , and to refine the estimates of the initial parameters returned by the grid search. The best-fitting evolutionary history investigated by the grid search is taken to be the solution that returns the absolute minimum value of the C statistic, whilst under the MCMC scheme the best-fitting set of parameters are taken to be the median values of the respective posterior probability distributions, with 1 − σ errors derived from the values that delineate the central 68.3 percent of the distribution. We choose this approach over the absolute minimum, as the latter strongly depends on the precise sampling of the parameter space by any given Markov chain. It was assumed that the orbit of the planet would monotonically shrink throughout the integration, thus placing a lower limit on the initial semi-major axis of a0,min = aobs . The upper limit on the same parameter was set at a0,max = 0.1 AU, a value which encompasses 95 percent of the current distribution for transiting planets. The ranges of the tidal quality factors were set to 4.0 6 log(Q′p ) 6 10.0 and 5.0 6 log(Q′s ) 6 10.0 such that the commonly accepted ranges of values (105 − 106 for Q′s , 107 − 108 for Q′p (Baraffe, 4 D. J. A. Brown et al Chabrier & Barman 2010)) were encompassed together with an additional section of parameter space. The initial eccentricity distribution for the ensemble of exoplanets is the subject of much discussion. A value of e0,max = 0.2 encompasses 93 percent of the present distribution and limits our parameter space to the region of validity for our tidal equations, as defined by Leconte et al. (2010). However the migration method by which hot Jupiters find themselves at such short orbital periods is still uncertain. One particular mechanism, Kozai scattering (Kozai 1962), can pump the eccentricity to large values, and it has been suggested that short period planets are captured at periastron from such highly eccentric orbits. If this is the case then a much higher upper limit on e0 could be justified. We must thus take into account the initial times at which we will be starting our integrations. Meibom & Mathieu (2005) carried out a survey of tidal circularisation in binary stars. They found that open clusters show a characteristic orbital period at which binaries with the most common initial eccentricity circularise (they define a circularised orbit as one with e = 0.01). This period was found to vary with the age of the cluster. However it is questionable whether exoplanets, owing to the much lower secondary mass, will tend to circularise at the same orbital period. In fact Hansen (2010) applied the formalism of Meibom & Mathieu (2005) to exoplanetary systems, finding that circularisation periods were generally much shorter than for stellar binaries. We therefore set e0,max = 0.8, a value encompassing 99 percent of the current distribution for transiting exoplanets, to allow for the possibility of eccentric orbits produced by Kozai oscillations. It is important to note that stellar and planetary obliquities greater than 0.0 are not considered by the integration methods used herein owing to the algebraic formulation that underlies our computational methods. We therefore assume that the orbital and spin angular momenta vectors are aligned, or anti-aligned. 3.1 Markov-chain Monte-Carlo simulation Under our MCMC optimisation scheme, the coordinates of the origin for the Markov chain are set to the centre of the parameter space. An exception is made if this lies within 3−σC of the coordinates of the best-fitting initial conditions, in which case the origin coordinates are set to lie outwith the 3 − σC contour in parameter space. The burn-in phase was judged to be complete when the minimum value of the test statistic from the current integration exceeded the median of all previous minimum test statistic values for the first time (Knutson et al. 2008). 4 THE WASP-18 SYSTEM The hot Jupiter WASP-18 b orbits the star HD10069, an F6 star with an effective temperature Tef f = 6400 ± 100 K and a V magnitude of 9.3, situated approximately 100 pc from the Earth. A fit to stellar isochrones using the model of Girardi et al. (2000) gives a stellar age of 0.630+0.950 −0.530 Gyr, which can be further constrained using the observed lithium abundance. Hellier et al. (2009) (hereafter H09) therefore assign the system an age of 0.5 − 1.5 Gyr, a value that the Table 1. System parameters for WASP-18, taken from H09 with the exception of the rotation period, which is derived from the observed v sin I under the assumption that the inclination of the stellar axis is 90o , the J-K colour, which was calculated using data taken from the SIMBAD on-line database, and λ (Triaud et al. 2010). Parameter Ms Rs J −K Tef f v sin I Prot,s Value Units 1.25 ± 0.13 Msun 1.216+0.067 −0.054 0.278 ± 0.032 6400 ± 100 10.77 ± 0.04 5.64 ± 0.28 Rsun K km s−1 days System age 0.630+0.950 −0.530 Gyr Mp 10.30 ± 0.69 Rp a e i Porb 1.106+0.072 −0.054 0.02026 ± 0.00068 0.0092 ± 0.0028 86.0 ± 2.5 0.94145299 ± 0.00000087 Mjup 5.0+3.1 −2.8 λ Rjup AU o days o current work attempts to improve upon. WASP-18 b was the first planet confirmed to have a period of less than a day, orbiting its host in just 0.94 days. Recent measurements of the Rossiter-McLaughlin effect by Triaud et al. (2010) have shown that the system is well aligned; the assumption of spin-orbit alignment required for our analysis is therefore justified. In the absence of a measured rotation period for WASP-18, we assume that the inclination of the stellar axis is 90o such that sin I = 1.0 and vrot = v sin I. An observed system age of 0.630 ± 0.530 Gyr was used for the calculation of the C-statistic for this system, in line with the age found from stellar isochrones by H09. WASP18 b has a radius of the order of that of Jupiter, but is significantly more massive. Lacking a means of calculating αp however, we set αp = αJupiter as a reasonable estimate. In light of the observed system age being consistent with the age of the Hyades, we compared Prot,s from Table 1 to Prot,s (t = tHyades ), as calculated for our derivation of κ for the system, finding the observed stellar rotation period to be shorter than the expected rotation period of 7.00 days at the age of the Hyades. We therefore set t0 = 150 Gyr≈ tM 35 , calculating Ωs,0 by scaling the period colour relation of Collier Cameron et al. (2009) to this age. The C-statistic values returned by the grid search were converted into probability densities according to P (C) = e− C2 2 , (20) and plotted as a set of projection maps (Fig. 1). These indicate regions of high probability, and therefore low points in the C-statistic surface, for each pair of initial parameters, and appear to show that there are a range of values for each parameter that could produce plausible evolutionary histories with a good quality of fit. This is particularly noticeable Are falling planets spinning up their host stars? 5 Figure 1. Projection maps of the test statistic probability density for WASP-18. The projections indicate distinct ranges of log(Q′s ), log(Q′p ) and a0 that give high probability, corresponding to a low value of the C-statistic. e0 provides a less clear cut high probability region, with several individual values giving good results. The best-fitting initial parameters were found to be a0 = 0.0427 AU, e0 = 0.245, log(Q′s ) = 7.00 and log(Q′p ) = 7.50, giving a system age of 0.654 Gyr. in e0 and log(Q′p )′ , indicating that these parameters have a lesser impact on the orbital and spin evolution of the system than a0 and log(Q′s ). This is, to a certain extent, unsurprising; the current orbital eccentricity for WASP-18 b is more uncertain than the orbital separation, and so will have less influence on the value of the C-statistic. It is also interesting to note the generally sharp transitions in Fig. 1 between regions with P (C) > 0.4 and those with P (C) ≈ 0. This indicates that although the values of the initial parameters are somewhat uncertain, they are generally well constrained to a smaller region of parameter space than we have allowed for. The best-fitting combination of initial parameters in these projection maps was found to be a0 = 0.0428 AU, e0 = 0.245, log(Q′s ) = 7.00 and log(Q′p ) = 7.50, indicating a system age of 0.654 Gyr in good agreement with the existing estimate derived from stellar model fitting. Using the period colour relation of Collier Cameron et al. (2009) we calculate an age of 0.387 ± 0.024 Gyr for WASP-18 using the derived rotation period, inconsistent with both the age from isochrone analysis and the result of our grid search. From the evolutionary track we derive an extremely short remaining lifetime of 0.006 Gyr, implying that the planet is currently spiralling in towards its host and is on the verge of reaching the Roche limit. This seems incredibly short, and would mean that we have managed to observe WASP-18 b in a very short window of opportunity. This result is, however, at oodds with the results that we have obtained from our second integration method, the MCMC algorithm. The posterior probability distributions produced by the MCMC integration scheme (Fig. 2) indicate significant correlation between many of the pairs of initial orbital parameters, in particular a0 and log(Q′s ) although it is also noticeable in several of the other distributions. Some correlation is to be expected; the orbital circularisation time-scale depends on energy dissipation within both the star and the planet (Miralda-Escudé 2002), and so a correlation between the eccentricity and the quality factors that govern tidal energy dissipation is unsurprising. The correlation between a0 and e0 may arise as a result of the strong coupling between eccentricity and separation noted by Jackson, Greenberg & Barnes (2008). The rapidity with which the orbital separation can decrease is intrinsically linked to the efficiency with which energy is dissipated within the system, a process that 6 D. J. A. Brown et al Figure 2. The posterior probability distributions of each pair of jump parameters produced by the MCMC analysis of the WASP-18 system. Note the clear correlations between several of the integration parameter pairs; these arise from the dependence of tidal energy dissipation on the orbital eccentricity, and from the strong coupling that exists between the orbital eccentricity and the orbital separation. The location of the best-fitting parameters in the parameter space explored are denoted by the dashed lines. is governed by the tidal quality factors. Moreover inspection of (12) indicates that ȧ ∝ Q−1 , so the form of the correlation between a0 and log(Q′s ) that is observed is also as expected. In fact it is surprising that little to no correlation is observed between a0 and log(Q′p ) The best-fitting parameters found using the MCMC code are set out in Table 2, and imply a stellar age of 0.579+0.305 −250 Gyr, in agreement with the existing age obtained through isochrone fitting, and the age derived from our grid search results. Fig. 3 displays the evolution of the stellar rotation for the set of best-fitting parameters; it is important to recall that the age estimate is not evaluated solely on the basis of this parameter, and that the observed orbital eccentricity and semi-major axis play a part as well. The value of log(Q′s ) returned by the MCMC-derived integration scheme does not agree with the best-fitting value obtained from the grid search of the defined parameter space. Nor do the values of a0 and e0 . We attribute this to the adoption of the set of median parameters from the MCMC scheme to avoid chance encounters with local minima, as well as the discrete nature of the grid search. Comparison of the solution with the lowest value of the C test statistic shows a value of log(Q′s ) that fits more closely to the grid search result, albeit still in disagreement, although we note that the minimum test statistic value for the grid search is an order of magnitude lower than that for the MCMC exploration of the available parameter space. We also note that although Figs. 1 and 2 are broadly similar in form, they differ somewhat in detail. The same correlations between parameters are visible in both figures, but the cutoff in log(Q′s ) occurs at a slightly greater value in Fig.2 than Fig. 1. Additionally, the MCMC algorithm explores a more narrow range of both a0 and e0 parameter space than the grid search, but does explore the lower end of log(Q′p ) space, a region which the grid search suggests gives poor results. We again attribute this to the differences in the manner in which the two algorithms explore the parameter space that we have delineated. Using (13) with the data from Table 2 and the bestfitting value of log(Q′s ) from the MCMC scheme, we estimate that WASP-18 b will reach its Roche limit 0.086 Gyr from now. The evolutionary track produced using the MCMC results implies a remaining lifetime of 0.076+0.742 −0.044 Gyr for the planet, consistent with this value. These are, respectively, Are falling planets spinning up their host stars? 7 Table 2. The initial orbital parameters and tidal quality factors of the best-fitting tidal evolution histories produced by the grid search and MCMC integration schemes for the WASP-18 system. The 1 − σ error bars for the MCMC derived values are estimated from the parameter values that encompass the central 68.3 percent of the final parameter distributions, and in some cases are inflated by the presence of short ‘tails’ in the distributions. Parameter Grid search MCMC Units a0 0.0427 0.0258+0.0052 −0.0044 AU e0 0.245 0.0399+0.1023 −0.0351 log(Q′s ) 7.00 8.21+0.90 −0.52 log(Q′p ) 7.50 7.77+1.54 −1.25 age 0.654 0.579+0.305 −0.250 Gyr tremain 0.006 0.076+0.790 −0.044 Gyr 2.40 and 2.12+20.73 −1.23 percent of the expected main sequence lifetime of the host star, which we estimate to be 3.58 Gyr. These remaining lifetimes strongly imply that the WASP-18 planetary system has an short life expectancy, as found by H09. Although still quite short, these values are far more reasonable than the value of 6 Myr derived from the grid search. The evolution of the stellar rotation period (Fig. 3) for the set of initial parameters returned by the MCMC integration scheme indicates that the rotational evolution of the star under the influence of tides initially differs little from our simulation in which the rotation period is governed purely by a Skumanich-type magnetic braking law. However the evolution of the rotation period rapidly diverges from this ideal case, with the rate of spin-down slowing gradually as the star ages and the planet migrates towards its host under the influence of tidal interactions. The final spiral-in of the planetary companion causes a rapid and substantial spin-up of the star; during this process the stellar rotation period is reduced by almost a factor of two from its maximum 6.26 day period to 3.16 days at the Roche limit, suggesting significant spin up of the host star. Furthermore, Fig. 3 clearly shows that, for the age estimate returned by our MCMC integration scheme, the observed rotation period of the host star is quite inconsistent with the rotation period of 7.96 days that would be expected at the same age if tidal interactions played no part in the evolution of stellar rotation. Although the best-fitting parameters that we have adopted show WASP-18 b in the process of its final spiral in, the true picture is slightly more ambiguous, as evidenced by the large upper uncertainty on the remaining lifetime. Within the 1 − σ parameter ranges there exist combinations of the initial parameters that produce evolutionary tracks of a different form to that displayed in Fig. 3. These evolutionary tracks follow the stellar rotation curve produced by a purely magnetic braking scenario for a much longer period of time, and exhibit little spin up of the host star. Furthermore, in semi-major axis space they indicate that the planet stays at approximately the same semi-major axis for Figure 3. The evolution of the rotation period of WASP-18 from the best-fitting initial conditions found by the MCMC exploration of the allowed parameter space (solid line), along with the evolution that would occur in the absence of tidal interactions (dashed line). For reference, the observed parameters are plotted with associated errors, and the best-fitting age given by the MCMC integration scheme is denoted (dotted line). The tide-governed evolution diverges from the ideal rotational evolution early on, with the effect of tidal interactions on the the rotation of the host star increasing with time until the rotation period dramatically decreases during the final spiral-in of the planetary body. the duration of the star’s main sequence lifetime. In these cases the age of the system lies towards the lower end of the adopted range, owing to the necessity of conforming to the calculated stellar rotation period, and the planet has a long remaining lifetime. We note however that these solutions form only a minority of the evolutionary tracks consistent with the parameter ranges that we adopt from the MCMC results. More prevalent were tracks of the same form as that displayed in Fig. 3, but with the spin-up taking place much later such that the stellar rotation implies an age consistent with the gyrochronological estimate and magnetic-braking only scenario. We attempted to further constrain the range of possible tidal quality factors by visually fitting to the observed parameters and their associated 1 − σ uncertainties. Starting from the best-fitting initial orbital separation and eccentricity, we investigated the evolution of the system across the range of values for both log(Q′s ) and log(Q′p ). We found that changing the value of log(Q′p ) makes no difference to the evolution of the orbital separation or rotation period, but has a strong effect on the eccentricity evolution of the planetary orbit. log(Q′s ), in contrast, strongly affects the evolution of all integration parameters. Moreover we found that changing the initial semi-major axis, within the adopted range, made a significant difference to the evolution of the stellar rotation period, whilst modifying e0 merely affected the evolution of the orbital eccentricity. We therefore conclude that the two most important parameters with respect to the stellar rotation are Q′s and a0 . We were unable to constrain the permissible ranges of the two tidal quality factors any further, owing to the range of values that a0 was able take, but found that for the majority of the possible combinations the rotational evolution of WASP-18 A gradually diverges from that expected 8 D. J. A. Brown et al Table 3. System parameters and 1−σ limits for WASP-19, taken from the free eccentricity fit of H10 with the exceptions of the J-K colour, which is derived from data taken from SIMBAD, and the rotation period (Collier Cameron, priv. comm.). Parameter Ms 0.95 ± 0.10 Msun 0.93+0.05 −0.04 0.43 5500 ± 100 4.0 ± 0.2 10.5 ± 0.2 Rsun System age Mp +9.0 5.5−4.5 1.14 ± 0.07 Rp e i Porb of a purely magnetic braking scenario, and that the star is eventually spun-up by a substantial amount during the final in-spiral of its planetary companion. Fig. 4a shows the rotational period evolution resulting from several values of log(Q′s ) within the permissible range, with log(Q′p ), a0 and e0 set to values adopted from the MCMC search. 5 THE WASP-19 SYSTEM The transiting hot Jupiter WASP-19 b orbits a late G-dwarf star with an orbital period of 0.7888399 ± 0.0000008 days. The host star has been measured to have Tef f = 5500 ± 100 K, periodic sinusoidal flux variations that indicate a Units Rs J −K Tef f v sin I Prot a Figure 4. Evolutionary tracks for WASP-18 produced using several different values of log(Q′s ) within the range consistent with the observed system parameters. The tracks were calculated using the values of e0 , a0 and log(Q′p ) adopted from the MCMC solution and, from left to right (at the end of the track): log(Q′s ) = 8.02, log(Q′s ) = 8.21, log(Q′s ) = 8.70 and log(Q′s ) = 9.11. Upper panel: The evolution of stellar rotation period with time for WASP-18. The dashed line shows the evolution expected from a Skumanichtype magnetic braking law. Changing the value of Q′p has no effect on the evolution of the rotation period with time. Lower panel: The evolution of the C test statistic with time for the tracks in the upper panel. The time at which C is a minimum is taken to be the age of the system given the set of initial parameters used for that evolutionary track. Value K km s−1 days 0.0164+0.0005 −0.0006 Gyr Mjup Rjup AU 0.02+0.02 −0.01 80.8 ± 0.8 0.7888399 ± 0.0000008 days 1.28 ± 0.07 detectable level of intrinsic variability and activity, and a metallicity slightly greater than solar at [M/H] = 0.1 ± 0.1 dex (Hebb et al. 2010) (hereafter H10). Measurement of the Rossiter-McLaughlin angle for the system indicates that it is aligned (Hellier et al. 2011). The age of the system is still somewhat uncertain, with H10 able to determine a constraint of > 1.0 Gyr. We note that several studies (e.g. Hansen (2010); Weidner & Horne (2010)), as well as the exoplanet encyclopedia, cite an age of 0.6 ± 0.5 Gyr for the star and attribute it to H10. We believe that this originates from the abstract of H10 which quotes that age as one of two possibilities. However the text of that paper favours the older age constraint. For our analysis we used an age of tage = 5.5 ± 4.5 Gyr, consistent with the isochronal fit and lower bound on the age quoted by H10, and set αp = αJupiter = 0.26401. Comparing the stellar rotation period from Table 3 to the rotation period expected at the age of the Hyades, we found that Prot,s > Prot,Hyades = 8.60 days and thus set t0 = tHyades and Ωs,0 = Ωs,Hyades . Fig. 5 shows the C-statistic probability density projection maps produced from the grid search results for the WASP-19 system. They indicate a fairly small region of high probability for log(Q′s ) and log(Q′p ), but appear to show that there are broader ranges of both e0 and a0 that give strong solutions. The absolute maximum probability, and thus the minimum value of the C statistic, was found to occur when the initial conditions were a0 = 0.0939 AU, e0 = 0.735, log(Q′s ) = 6.25 and log(Q′p ) = 10.0, giving a system age of 2.11 Gyr that is consistent with the existing constraint of > 1.0 Gyr from H10. Using the period-colour relation of Collier Cameron et al. (2009) we calculate a gyrochronological system age of 0.899±0.037 Gyr, immediately suggesting that gyrochronological analysis is inappropriate for this system. Comparing the sets of projection maps for the two systems, it is apparent that there are similarities in the forms Are falling planets spinning up their host stars? 9 Figure 5. C-statistic probability density projection maps for WASP-19 b. A small region of high probability is clearly visible in log(Q′s ) and log(Q′p ) parameter space, but in a0 and e0 the range of values that give high probability solutions are much broader. The point of maximum probability, indicating the minimum in the C-statistic surface, was found to occur at parameter space coordinates of a0 = 0.0939 AU, e0 = 0.735, log(Q′s ) = 6.25 and log(Q′p ) = 10.0, giving a system age of 2.11 Gyr. of Figs. 1 and 5, particularly in the a0 − log(Q′s ) and e0 − a0 maps. However the range of probability values covered by the greyscale is slightly narrower for WASP-19, where the maximum probability is approximately 0.8, compared to WASP-18, where the maximum probability is almost 1.0. From this we deduce that the best-fitting orbital solution for WASP-19 is less certain than our best-fitting grid search solution for WASP-18. This lower maximum probability, coupled with the form of the projection maps, also suggests that the range of tidal quality factors for which valid orbital solutions exist is greater for this system than it was for WASP-18. The posterior probability distributions produced by the MCMC algorithm exhibit a strong correlation between a0 and log(Q′s ), of a similar form to that observed for the WASP-18 system albeit extended to lower values of log(Q′s ). Also present is a correlation between e0 and a0 that is somewhat similar to that apparent in Fig. 2, but the other correlations that were present in the distributions for WASP-18 are absent for WASP-19. The region of parameter space explored by the MCMC algorithm is also substantially greater for the WASP-19 system than it was for the WASP-18 sys- tem, although this is not unexpected after the comparative appearance of Figs. 1 and 5. It is also interesting to note the relative lack of agreement between the results from the grid search and the MCMC search methods for this system; the only search parameter for which the two sets of results agree is log(Q′s ). This is also apparent from inspection of Figs 5 and 6, as the region of greatest probability in the former does not match the location of the adopted parameters in the latter. As for the WASP-18 system, we attribute this to the adoption of the set of median parameters from the MCMC scheme; comparison of the solution with the lowest value of the C test statistic shows that it more closely fits with the location of the high probability region in Fig. 5. Fortunately the stellar tidal quality factor, the most important parameter when considering the evolution of the stellar rotation, is relatively well constrained in our MCMC-derived solution. The probability distribution for the initial orbital eccentricity of the WASP-19 system also shows weaker clustering than for the WASP-18 system, and the clustering that is present encompasses a much smaller range of eccentricity values above e0 = 0. This may be symptomatic of the fact 10 D. J. A. Brown et al Figure 6. The posterior probability distributions of each pair of integration parameters produced by the MCMC analysis of the WASP19 system. The correlations between the parameters are substantially different than those for the WASP-18 system, with that between log(Q′s ) and a0 being the only particularly noticeable correlation. The location of the best-fitting parameters in the parameter space explored are denoted by the dashed lines. that the eccentricity of the system is not well known; previous analysis of the parameters of the system found little difference in the fit to the observed transit lightcurve between the cases in which the orbit was forced to be circular, and in which the eccentricity was allowed to float (H10). In the latter case, the best-fitting eccentricity value was only 0.02+0.02 −0.01 , and it is this value that we used in our C statistic calculations. The best-fitting parameters given by the MCMC integration scheme, set out in Table 4, imply a stellar age of 1.60+2.84 −0.79 Gyr, in broad agreement with the loose constraint of age > 1.0 Gyr found by H10. The gyrochronological age that we have calculated for the system also agrees with this age estimate but lies at the lower end of the range, suggesting that gyrochronology provides a possible, if unlikely estimate for the age of WASP-19 A. The data in Table 4 further imply a remaining lifetime for the planet of +1.1073 0.0067−0.0061 Gyr, a mere 0.4 percent of the estimated system age. This strongly suggests that WASP-19 b is in the final, spiral-in stage of its orbital evolution, a conclusion supported by the stellar rotation period evolutionary track that results from integration of the best-fitting parameters (Fig. 7). As with the WASP-18 system we attempted to further constrain the range of possible tidal quality factors by visually fitting to the observed parameters; in this case we were able to reduce the upper limit on log(Q′s ) and raise the lower limit on a0 . The updated values are given in Table 4. Updating the limits on these parameters have no effect on the 1 − σ limits of the age or remaining lifetime. The observed spin rate of WASP-19 b is inconsistent with evolution governed only by magnetic braking; Fig. 7 shows that, at the system age implied by our MCMC solution, the observed period of 10.5 ± 0.2 days is substantially less than the period of 15.9 days implied by our magnetic braking only model. It is possible that WASP-19 b is in the process of spiraling-in to the Roche limit, spinning up its host as it does so. Fig. 8 shows the rotational period evolution resulting from log(Q′s ) values within the 1 − σ limits returned by the MCMC search scheme. These reinforce the idea that WASP19 b is on the verge of reaching its Roche limit; for the bestfitting tidal quality factors the end of the evolutionary track Are falling planets spinning up their host stars? 11 Table 4. The initial orbital parameters and tidal quality factors of the best-fitting tidal evolution histories produced by the grid search and MCMC integration schemes for the WASP-19 system. The 1 − σ error bars for the MCMC derived values are estimated from the parameter values that encompass the central 68.3 percent of the final parameter distributions. Parameter Grid search MCMC Updated Units 0.0939 0.0317+0.0228 −0.0146 0.0317+0.0228 −0.0089 AU e0 0.735 log(Q′s ) 6.25 0.0017+0.0597 −0.0016 6.47+2.19 −0.95 6.75+1.86 −1.77 1.60+2.84 −0.79 +1.1073 0.0067−0.0061 a0 log(Q′p ) 10.0 age 2.11 tremain 0.0038 +0.67 6.47−0.95 Gyr Gyr Figure 8. The evolution of stellar rotation period with time for WASP-19 for a range of log(Q′s ) values. The dashed line shows the evolution expected from a Skumanich-type magnetic braking law, whilst the three solid lines are the evolution that results from the best-fitting e0 , a0 and log(Q′p ) given by the MCMC algorithm, and, from left to right (at the observed system age), log(Q′s ) = 6.14, log(Q′s ) = 6.47 and log(Q′s ) = 6.52. These tracks imply that WASP-19 b is very close to the Roche limit, and will reach it in less than 10 Myr. Figure 7. The evolution of the rotation period of WASP-19 from the best-fitting initial conditions found by the MCMC exploration of the allowed parameter space, along with the evolution that would occur in the absence of tidal interactions (dashed line). The tide-governed evolution follows the gradual spin-down of the isolated evolution closely until spiral-in of the planet begins, at which point the stellar rotation period rapidly decreases. The observed rotation period is plotted with associated errors, and implies that WASP-19 b is in the final spiral-in stage of its orbital evolution. Furthermore it is irreconcilable with the expected spin down from magnetic braking alone, and can only be explained by invoking spin-up during tidal interactions. is at Prot = 10.5 days compared to the observed Prot = 10.5 ± 0.2 day, placing the planet precisely on the Roche limit and suggesting that we have been fortunate to observe WASP-19 b at all. Much as with the WASP-18 system though, this is not the entire picture presented by our results. The upper limit on tremain implied by our adopted solution is 1.114 Gyr, a significant 69.6 percent of the estimated age and 13.7 percent of the MS lifetime of the host star. This strongly suggests that there exist solutions that more closely follow the stellar rotation evolution expected when only magnetic braking is acting on the star. To investigate this we characterised the orbital evolutionary tracks for each log(Q′s ) − a0 node from the grid search that was consistent with our results from the MCMC algorithm. Since our investigation into WASP-18 b’s future evolution led us to the conclusion that e0 and log(Q′p ) have little effect on the evolution of the stellar rotation, we set these parameters to the MCMC values from table 4. We found that for each value of a0 there was only a very narrow range of log(Q′s ) values that gave solutions with an age > 1.0 Gyr and a rotation period at that age consistent with the observed value. In fact only 13 percent of the nodes investigated gave such solutions. The majority of the nodes gave solutions that either did not reach an age of 1.0 Gyr, instead showing the distinctive spiral-in signature at younger ages, or that followed the magnetic braking-only track for much longer such that the spin-up induced by the final spiral-in was insufficient to reduce the stellar rotation period back to the observed value. Tracks that fell into this latter category often returned ages consistent with our gyrochronological calculation, although some did return older ages consistent with the 1.0 Gyr lower limit owing to the influence of the orbital eccentricity and semi-major axis. This raises questions concerning WASP-19 A’s true age. There is evidence from isochrone fitting that the star is older than 1 Gyr. This is supported by the investigation of H10 into the space velocity of the star, which found that 65 percent of simulated stars with similar stellar properties in a small volume around WASP-19 A were older than 1.0 Gyr. But the rotation period of WASP-19 A implies a younger age of 0.899 Gyr through gyrochronology. The lithium abundance of log(A[Li]) < 1.0 quoted by H10, while consistent with either estimate as it only constrains the age to greater than 0.6 Gyr, would appear to support a younger age, as the presence of any lithium tends to rule out a more evolved star. There therefore appear to be two main possibilities; either the star is old and has been spun up by 12 D. J. A. Brown et al the infall of the planet, or the star is younger, still following its natural spin down and more dense and/or hotter than expected for its age. Both of these are consistent with the results that we have obtained from our search methods, so which is the more plausible? 5.1 Is W19 young or old? We consider our results for WASP-19 in the context of a larger ensemble of exoplanetary systems with similar stellar and planetary properties, assuming that the value of log(Q′s ) = 6.47 adopted for WASP-19 is applicable to the other systems under consideration. We find that for 7 of the 9 planetary systems in the sample the remaining lifetime is < 1 Gyr. Furthermore we find that WASP-4 b has a remaining lifetime of 0.017 Gyr, and calculate that it lies at a distance of only 1.4 times the Roche limit from its host star; this places it in apparently similar circumstances to WASP19 b. WASP-2 b and WASP-10 b have remaining lifetimes at least one order of magnitude greater than those of the rest of the sample, but they also have significantly longer orbital periods. It is interesting to consider the range of the stellar tidal quality factor that produces certain threshold values of tremain . We find that, if we exclude the longer period WASP-2 and WASP-10 systems, log(Q′s ) > 7.10 leads to an estimated tremain > 1.0 Gyr for all planets in the sample, whilst a remaining lifetime greater than 0.1 Gyr requires log(Q′s ) > 6.12. For WASPs -2 and -10 we find that the values are 6.24 and 5.24, and 6.22 and 5.22 respectively. The tidal evolution solution that we presented previously would seem to support the hypothesis that the star is old, with stellar spin up accounting for the rotation period. But the results in Table 5 would seem to add weight to the idea that the system is actually younger than expected. The likelihood of observing one system in such a condition is low, but to observe two such systems as Table 5 implies that we have done with WASPs -4 and -19, seems incredulous. Considering the remaining lifetimes in the context of the total planetary lifetime suggests a different scenario. 13 can also be used to calculate the total planetary lifetime if the initial semi-major axis is known. We calculated the total planetary lifetimes assuming that the value adopted for WASP-19 b, a0 = 0.0317 AU applied unilaterally across our sample, and in Fig. 9 plot remaining lifetime as a function of total lifetime. It is immediately apparent that WASP19 b is a special case, being clearly separated from the rest of the sample. The remaining lifetime of WASP-19 b is only between 1 and 2 percent of the total lifetime, whilst for all of the other planets in the sample this figure is greater than 10 percent. The short remaining lifetime of WASP-4 b is thus somewhat misleading, as it represents a significant portion (≈ 12.6 percent) of the total lifetime of the planet. There are several ways in which we can analyse the remaining lifetimes calculated here, one of which is to examine the probability that we have managed to observe the system in its present configuration. We calculated this for the systems in Table 5, disregarding those with no literature age estimate. For WASP-19 we calculate the probability using both our own age adopted age estimate, and the younger age of 0.6 ± 0.5 that can be found in the literature. Fig. 10 displays the results for a range of log(Q′s ) values consistent with Figure 9. Remaining lifetime as a function of total planetary lifetime for the sample of transiting planets in Table 5, as calculated using 13 assuming that the values of log(Q′s ) and a0 adopted for the WASP-19 system from our MCMC results apply universally. We plot the results for log(Q′s ) = 5.52 (squares), log(Q′s ) = 6.47 (circles) and log(Q′s ) = 7.14 (triangles). The data for WASP-19 b are denoted by filled symbols. The lines represent tremain = 1 percent (red, dash), 5 percent (green, dot-dash) and 10 percent (dark blue, dot) of the total planetary lifetime. WASP19 b is clearly separate from the rest of the sample, indicating that it is in a unique situation. the 1 − σ limits derived from the MCMC posterior probability distributions. For the value of log(Q′s ) adopted from the MCMC solution we find that two systems, including WASP19 at our older age estimate, have a probability of observation of less than 1 percent. Using the younger age estimate for WASP-19 pushes the probability up to 2 percent, which is more plausible (see analysis in Hellier et al. (2011)). Increasing the value of log(Q′s ) to our upper limit of 7.14 increases the probability of observing WASP-19 b to approximately 3.5 percent for the older age, and approximately 9 percent for the younger age. However if we consider the lower limit of log(Q′s ) = 5.52 then four systems, including both ages of WASP-19, have observation probabilities less than 1 percent. This is an entirely implausible situation. However this analysis assumes that all of the planets in our sample started at the same distance from their respective host stars, and that they experience tidal interactions of the same strength. It seems somewhat unlikely that this accurately represents reality, as the planetary systems in our sample show significant variation in their properties. Previous studies (Matsumura, Peale & Rasio 2010; Hansen 2010) have found that different planetary systems are likely to experience different strengths of stellar and planetary tide, so describing these systems with a single value is almost certainly unphysical. This would, of course, mean that the observational probabilities for several of these hot Jupiters could be significantly greater. To further muddy the waters we add evidence for an old age by carrying out additional stellar model fits. We used the updated Padova models (Marigo et al. 2008), which imply +4.65 an age of 7.94−2.93 Gyr, and the latest version of the YonseiYale isochrones (Demarque et al. 2004), which return an age of 7+2 −3 Gyr. Both estimates agree with the earlier isochrone Are falling planets spinning up their host stars? 13 Table 5. A comparison of the WASP-19 system to a sample of transiting hot Jupiter systems with Mp > MJ in close orbits around stars similar to, or cooler than, WASP-19 A. We calculate the remaining lifetime for each planet using equation 13, assuming that the stellar tidal quality factor of log(Q′s ) = 6.47 adopted for WASP-19 A is applicable across the entire sample. Note that the remaining lifetime for WASP-19 b quoted here does not agree with the value derived from the evolutionary track displayed in Fig. 7. System Mp /MJ Rp /RJ Ms /M⊙ Rs /R⊙ a /AU P /days tremain /Gyr Age /Gyr WASP-19 1.14 1.28 0.95 0.93 0.0164 0.7888399 0.012 +2.84 1.60−0.79 WASP-2 0.847 1.079 0.84 0.834 0.03138 2.1522254 1.708 +8.1 11.9−4.3 WASP-4 WASP-5 1.12 1.63 1.416 1.171 0.93 1.00 1.365 1.15 0.023 0.02729 1.3382282 1.6284246 0.017 0.080 a 7.0+5.2 −4.5 3.0 ± 1.4 b WASP-10 CoRoT-1 CoRoT-2 OGLE-TR-113 3.06 1.03 3.31 1.24 1.08 1.49 1.465 1.11 0.71 0.95 0.97 0.78 0.783 1.11 0.902 0.77 0.0371 0.0254 0.0281 0.0229 3.0927616 1.5089557 1.7429964 1.4324772 1.818 0.093 0.159 0.224 0.8 ± 0.2 0.111 a 0.1+0.7 −0.0 TrES-3 1.91 (a)Southworth (2010) 1.305 0.92 (b)Anderson et al. (2008) 0.813 0.0226 1.30618608 (c)Christian et al. (2009) > 0.7 a c d (d)Melo et al. (2006) sity requires a change to either the stellar mass or radius, or both, either of which could have a dramatic effect on our estimates of the planetary parameters. This is notwithstanding the more likely case, which is that both the temperature and density are in need of a slight adjustment. A younger age for the star does not seem, therefore, outwith the realms of possibility. Figure 10. The probability of observing the systems in Table 5 in their present configuration, as a function of log(Q′s ) within the range encompassed by the 1 − σ errors from the MCMC results. Systems without a stellar age in the literature were disregarded. The data are offset slightly from the true value of log(Q′s ) to aid clarity for the error bars. Systems for which P > 1 are plotted at P = 1. WASP-19 is denoted by filled symbols; triangles for the probability using the age estimate from the MCMC results, squares for the probability calculated using the estimate of 0.6 ± 0.5 Gyr prevalent in the literature. At log(Q′s ) = 6.47 two systems, including WASP-19 with its older age, show a probability of observation that is less than 1 percent. Increasing the stellar tidal quality factor increases the probability that we have observed the systems in the state implied by their remaining lifetimes. fit of H10. We also investigated the values of Tef f and ρs required for a stellar age of 0.6 ± 0.5 Gyr in these two models, assuming the same metallicity; the results are given in Table 6. The correction required to the effective temperature if the stellar density is correct is not improbable, but the converse is not necessarily true. Modifying the stellar den- For completeness we also attempted to force an age of 0.6 Gyr by adjusting the stellar metallicity whilst leaving the effective temperature and stellar density at the values from H10. We found that the upper limit of Z = 0.03 imposed on the Padova isochrones prevented us from reaching such a low age, leaving us with a lowest age estimate of 3.98 ± 3.96 Gyr, which could be consistent with a star the age of the Hyades. When fitting to the Yonsei-Yale isochrones we found that a value of [F e/H] = 0.32 gave an age of 2.5+3.5 −1.9 Gyr, just consistent with the younger age that we were looking for. Fitting to the 0.6 Gyr isochrone itself would require a much greater metallicity, but given that this value of the iron abundance is already 0.05 dex greater than the maximum elemental abundance found by H10 (specifically the upper 1 − σ limit on the Calcium abundance), it is clear that such a fit would be utterly inconsistent with the current spectral analysis. It therefore seems that a combination of slightly increased elemental iron abundance, greater stellar density and higher effective temperature could produce stellar model fits more consistent with the young stellar age suggested by gyrochronology. Indeed if one is using the Padova models then such a combination lies just within the existing 1 − σ ranges for the three parameters. As noted previously, the constraint placed on the age by the lithium abundance has the potential to rule out the young stellar age. The detection of a substantial lithium abundance in a stellar spectrum strongly implies a young age for the star, even one that hosts planets, as recent work by Baumann et al. (2010) shows that there is no correlation between the presence of planets and reduced stellar lithium 14 D. J. A. Brown et al Table 6. The age estimates obtained through fitting the stellar parameters of WASP-19 A to stellar models. Adjusting the stellar density, effective temperature or iron abundance in isolation lowers the age that the model fit returns. Obtaining an age of 0.6 Gyr required the parameters to be increased beyond their 1 − σ limits when done in isolation. If adjusted as a set, the required age could be obtained with more plausible values. Model ρs /ρ⊙ Tef f /K [Fe/H] age /Gyr Padova 1.13 ± 0.12 5500 ± 100 0.02 7.94+4.65 −2.93 Padova 1.45 ± 0.10 5500 ± 100 0.02 +3.35 0.631−0.627 Padova 1.13 ± 0.12 5805 ± 100 0.02 +3.35 0.631−0.626 Padova 1.13 ± 0.12 5500 ± 100 0.2 3.98 ± 3.96 Padova 1.25 ± 0.10 5585 ± 100 0.12 0.631+3.35 −0.625 Yonsei-Yale 1.13 ± 0.12 5500 ± 100 0.02 7+2 −3 Yonsei-Yale 1.52 ± 0.10 5500 ± 100 0.02 0.60+1.90 −0.54 Yonsei-Yale 1.13 ± 0.12 5885 ± 100 0.02 0.60+1.90 −0.54 Yonsei-Yale 1.13 ± 0.12 5500 ± 100 0.32 2.5+3.5 −1.9 Yonsei-Yale 1.30 ± 0.10 5625 ± 100 0.12 0.60+1.90 −0.40 content. The original analysis in H10 utilised only 34 spectra from the CORALIE spectrograph (Queloz et al. 2000), but thanks to an ongoing radial velocity observation program we now have access to additional data on WASP-19 A. Since the publication of H10 a further 3 spectra have been taken using CORALIE, and 36 spectral measurements have been obtained using the HARPS high precision echelle spectrograph (Mayor et al. 2003). We co-added the individual spectra into a single spectrum with a higher signal-to-noise ratio, from which we were able to improve the constraint on the lithium abundance. Hints of a lithium line were present but only at the level of the spectral noise, limiting us to an upper limit of log(ǫLi ) < 0.5. This allows us to place a lower limit on the age of WASP-19 A of 2.0 Gyr (Sestito & Randich 2005), making an older age for the star more palatable. The final piece of evidence to be considered is the proximity of WASP-19,to the Roche limit. Using the ob¯ served system parameters and their associated uncertainties with the formulation of Eggleton (1983), we calculate that the Roche limit for WASP-19 b is at 0.0148 ± 0.0025 AU, in agreement with the observed orbital separation of 0.0164+0.005 −0.0006 AU. Hellier et al. (2011) place the planet slightly further from the star at 1.21 times the Roche tidal radius, but we note that they have used a different formulation for calculating the limit. We conclude that the evidence is in favour of WASP19 A being old. Whilst there is evidence for the star being young, it is more circumstantial than that which points to an older star. The upper limit on the lithium abundance and the results from stellar model fitting in particular point towards an age in excess of 1.0 Gyr. If WASP-19 A is indeed old, then the exploration of tide-governed evolution presented herein suggests that WASP-19 b has spun up its host star, and might be in the final stages of spiralling into the Roche limit. During our exploration of possible evolutionary histories we found that those which returned a stellar age of > 1.0 Gyr tended to exhibit a very short remaining lifetime. There were, however, some histories that married an older age for the star to a long remaining lifetime for the planet, so we are unable to completely rule out that scenario. 6 SUMMARY AND CONCLUSIONS We have calculated an age for the WASP-18 system of +0.950 0.579+0.305 −0.250 Gyr, in agreement with the 0.630−0.530 Gyr age found from stellar isochrones by H09. Using an MCMC algorithm we find tidal quality factors of log(Q′s ) = 8.21+0.90 −0.52 and log(Q′p ) = 7.77+1.54 −1.25 . Our results imply that WASP-18 b will reach the Roche limit in 0.076+0.790 −0.044 Gyr, and that in most cases it will cause its host to spin up as it does so. We are unable to place stronger constraints the status of the system with respect to planetary infall owing to the range of evolutionary histories that we find to be compatible with its observed parameters, but a large number of the evolutionary tracks investigated imply that the planet is in the process of spiralling in to its host star. Our results for the WASP-18 system seem roughly consistent with the work of Matsumura, Peale & Rasio (2010). For the WASP-19 system we found tidal quality factors +0.67 of log(Q′s ) = 6.47−0.95 and log(Q′p ) = 6.75+1.86 −1.77 . These values give a stellar age of 1.60+2.84 Gyr, in broad agreement −0.79 with the constraint of age > 1.0 Gyr found by H10, and imply a remaining lifetime of 0.0067+1.1073 −0.0061 Gyr. We have investigated the possibility that WASP-19 A is younger than previously estimated, in line with the predictions of gyrochronology. After considering the evidence for both the old and young stellar age possibilities, we conclude that the older age is more probable based on a reanalysis of the spectral lithium abundance and updated stellar model fits. We therefore suggest that WASP-19 b could be in the final stages of its spiral-in, and could be on the verge of reaching the Roche limit. We found that this scenario was more prevalent amongst those evolutionary histories with a stellar age > 1.0 Gyr, but that there were some instances in which the older stellar age coincided with a substantial remaining planetary lifetime. We are therefore unable to rule out the possibility that it will be some time before the planet falls into the star. Tidal interactions between these two hot Jupiters and their host stars will dramatically affect the evolution of the stellar rotation periods, counteracting and then reversing the spin-down that is expected from evolution according to a Skumanich-type magnetic braking law. The observed rotation periods are irreconcilable with such an evolution of the rotation period, and strongly suggest that falling in hot Jupiters cause their host stars to spin up during their inward, tidal interaction governed migration. It is worth noting that the our results seem to point toward a more diverse range of stellar tidal dissipation strength than is commonly considered in the literature. The ranges of log(Q′s ) that we have attributed to the two systems investigated herein are slightly disparate, for which there may be several possible explanations. We turn to the work of Pinsonneault, DePoy & Coffee (2001), which shows that the mass of the convective zone is a function of Tef f , and therefore of spectral type. From their fig. 1, we note that the effective temperature of WASP-18 A lies close to the point at which the mass of the convective zone becomes negligible, imply- Are falling planets spinning up their host stars? ing a convective zone mass of MCZ ≈ 0.001 Msun . WASP19 A, with its much lower effective temperature, would have MCZ ≈ 0.030 Msun . In our work we have assumed that the star rotates as a single body, but we have made no assumption about where the majority of tidal dissipation takes place. If this process occurs in the convective zone then the discrepancy in the masses, and therefore depths, of the convective zones of the two stars could provide an explanation for the disagreement as to the value of Q′s , with a larger convective zone allowing for more efficient dissipation and hence a smaller quality factor. This study also provides a warning against using gyrochronology to estimate the ages of hot Jupiter host stars. Owing to the tidal spin-up of its host star by the in-falling planet, the age that we have found for the WASP-19 planetary system is greater than the age found using gyrochronology alone. The situation for the WASP-18 system is less clear cut, but there is no doubt that gyrochronology will not provide an accurate estimate of the system age at all points during its evolution. We therefore suggest that care should be taken when applying gyrochronology to hot Jupiter systems. The two systems studied herein are extreme examples of this class of planet, with extremely short orbital periods and very close orbits. It is likely that there will be a critical point in semi-major axis space, beyond which tidal interactions will be sufficiently weak so as not to influence the stellar rotation period. For systems in this region gyrochronology may well work very well. But for systems such as WASP-18 and WASP-19 which have a < acrit , gyrochronology should be applied with care. Determination of this critical point will require further work, and it is likely to depend on a host of factors. Lanza (2010) also suggests that gyrochronology may not always provide accurate age estimates for exoplanetary systems, finding that plotting Prot t−α as a function of Tef f for planet hosting stars gives a poor fit to the period-colour relation of Barnes (2007). This implies that exoplanet host stars are systematically faster rotators than stars with similar ages and properties that do not appear to have any associated planets. Lanza also finds that the rotation period evolution of F- and G-type planet hosts does not differ substantially from similar stars without hot Jupiters, a conclusion with which we disagree although we note that the final spiral-in of the planet does not appear to have been considered. The substantially reduced rotation period that results from tidal spin-up may provide a means of detecting stars that have either been planet hosts in the past, or that have unseen planetary companions that are in the process of spiraling in to the Roche limit. Measurement of such an anomalous rotation period would provide a strong indication of the current or previous existence of a hot Jupiter around that star. Searching for such systems could help to pinpoint targets for extra-solar planet searches. ACKNOWLEDGMENTS DJAB would like to thank Keith Horne and Moira Jardine for useful comments and suggestions made during an early draft of this manuscript. The authors are grateful to the referee Brad Hansen, for insightful and thought- 15 provoking comments that greatly improved the quality of the manuscript. The WASP Consortium consists of representatives from the Universities of Keele, Leicester, The Open University, Queens University Belfast and St Andrews, along with the Isaac Newton Group (La Palma) and the Instituto de Astrofisca de Canarias (Tenerife). The SuperWASP and WASP-S Cameras were constructed and operated with funds made available from Consortium Universities and PPARC/STFC. This research has made use of NASA’s Astrophysics Data System Bibliographic Services, the ArXiV preprint service hosted by Cornell University, and the SIMBAD database, operated at CDS, Strasbourg, France. REFERENCES Alves S., Do Nascimento Jr. J. D., de Medeiros J. R. 2010, MNRAS, 1348 Anderson D. R. et al., 2008, MNRAS, 387, L4 Baraffe I, Chabrier G., Barman T., 2010, RPPh, 73, 016901 Barker A. J, Ogilvie G. I., 2009, MNRAS, 395, 2268 Barnes S. A., 2007, ApJ, 669, 1167 Baumann P., Ramı́rez I., Meléndez J., Asplund M., Lind K., 2010, A&A, 519, A87 Christian D. J. et al., 2009, MNRAS, 392, 1585 Claret A., 2004, A&A, 424, 919 Claret A., 2005, A&A, 440, 647 Claret A., 2006, A&A, 453, 769 Claret A., 2007, A&A, 467, 1389 Collier Cameron A. et al., 2009, MNRAS, 400, 451 Demarque P. et al., 2004, ApJS, 155, 667 Dobbs-Dixon I., Lin D. N. C., Mardling R. A., 2004, ApJ, 610, 464 Eggleton P. P., 1983, ApJ, 268, 368 Eggleton P. P., Kiseleva L. G., Hut P., 1998, ApJ, 499, 853 Girardi L. et al., 2000, AAPS, 141, 371 Goldreich P., Soter S., 1966, Icarus, 5, 375 Gu P.-G., Lin D. N. C., Bodenheimer P. H., 2003, ApJ, 588, 509 Hansen B. M. S., 2010, ApJ, 723, 285 Hebb L. et al., 2010, ApJ, 708, 224 Hellier C. et al., 2009, Nature, 460,1098 Hellier C. et al., 2011, ArXiv e-prints: 1101.3293 Hut P., 1980, A&A, 92, 167 Hut P., 1981, A&A, 99, 126 Jackson B., Greenberg R., Barnes R., 2008, ApJ, 678, 1396 Jackson B., Barnes R., Greenberg R., 2009, ApJ, 698, 1357 Knutson H. A. et al., 2008, ApJ, 673, 526 Kozai Y., 1962, AJ, 67, 591 Lanza A. F., 2010, A&A, 512, A77 Leconte J. et al., 2010, A&A, 516, A64 Levrard B., Winisdoerffer C., Chabrier G., 2009, ApJ, 692, L9 Mardling R. A., Lin D. N. C., 2002, ApJ, 573, 829 Marigo P. et al., 2008, A&A, 482, 883 Matsumura S., Peale S. J., Rasio F. A., 2010, ApJ, 725, 1995 Mayor M. et al., 2003, The Messenger, 114, 20 Meibom S., Mathieu R. D., 2005, ApJ, 620, 970 Melo C. et al., 2006, A&A, 460, 251 Miralda-Escudé J., 2002, ApJ, 564, 1019 Pinsonneault M. H., DePoy D. L., Coffee M., 2001, ApJL, 556, L59 Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in C. The art of scientific computing, 2nd ed. Cambridge University Press, Queloz D. et al., 2000, in Bergeron J., Renzini A., eds, From Extrasolar Planets to Cosmology: The VLT Opening Symposium, p. 548 16 D. J. A. Brown et al Pont F., 2009, MNRAS, 396, 1789 Radick R. R. et al., 1987, ApJ, 321, 459 Sasselov D. D., 2003, ApJ, 596, 1327 Sestito P., Randich S., 2005, A&A, 442, 615 Skumanich A., 1972, ApJ, 171, 565 Southworth J., 2010, MNRAS, 408, 1689 Sozzetti A. et al., 2007, ApJ, 664, 1190 Triaud A et al., 2010, A&A, in press Weber E. J., Davis Jr. L. 1967, ApJ, 148, 217 Weidner C., Horne K., 2010, A&A, 521, A76 This paper has been typeset from a TEX/ LATEX file prepared by the author.