Chemo-kinematic Ages of Eccentric-planet-hosting M Dwarf Stars

In the following sections, we describe a data-driven approach to estimate the kinematic prior and [Ti/Fe] likelihood. For this approach, we require an unbiased sample of stars with known age, kinematics, and [Ti/Fe]. The Veyette et al. (2017) method to measure [Ti/Fe] was calibrated to match the Brewer et al. (2016, hereafter B16) catalog of detailed abundances for 1617 FGK stars. Therefore, to ensure consistency and reduce systematic errors, we used this same catalog to develop our kinematic–[Ti/Fe]—age model. To estimate the ages of the B16 stars, we used the isochrones package (Morton 2015) with the MIST stellar evolution models (Choi et al. 2016; Dotter 2016). We describe this process in detail in the Appendix.

B16 fit for and removed from their abundance estimates any systematic trends with temperature; however, this trend was only assessed over a limited range of T_eff and log g. We found that systematic trends in [Ti/Fe] still existed for stars with T_eff > 6100 K and log g < 3.6, so we excluded those stars from further analysis. We also excluded stars with best-fit A_V values >0.1. All stars in this sample are solar-neighborhood stars, and we do not expect significant extinction. These cuts, combined with the initial requirement that stars have a parallax measurement available in the literature and convergence criteria as described in the Appendix, leave 672 FGK stars for which we have reliable [Ti/Fe] and age estimates. Figure 1 shows the general trend of increasing [Ti/Fe] with increasing age.

Zoom In Zoom Out Reset image size

Of these 672 stars, 658 have radial velocities and full five-parameter astrometric solutions in Gaia DR2 (Gaia Collaboration et al. 2018; Lindegren et al. 2018), which we used to calculate the U, V, and W peculiar velocities of each star (calculations based on code adapted from Rodriguez 2016).

Almeida-Fernandes & Rocha-Pinto (2018, hereafter A18) introduced a method for estimating the age of a star based on its peculiar velocities alone. They modeled the components of the velocity ellipsoid of field stars as Gaussian distributions with age-dependant dispersion. They used the Geneva-Copenhagen Survey (GCS; Nordström et al. 2004; Casagrande et al. 2011) to fit for the dispersion of these distributions as functions of age, as well as the V component of the solar motion, ${V}_{\odot }^{{\prime} }$, and the vertex deviation, ℓ_v. Evaluating the product of the three distributions at a given age produces the likelihood function for the measured U, V, and W velocities. We employ a prior probability distribution for a star's age based on the posterior probability distribution given by Equation (10) of A18 (method UVW),

$$\begin{eqnarray}&&p(\tau | U,V,W)\propto \displaystyle \prod _{i=1,2,3}\displaystyle \frac{1}{\sqrt{2\pi }{\sigma }_{i}(\tau )}\exp \left(-\displaystyle \frac{{v}_{i}^{2}}{2{\sigma }_{i}{(\tau )}^{2}}\right),\end{eqnarray} \tag{ 2 }$$

where v₁, v₂, and v₃ are the star's velocities in terms of the components of the velocity ellipsoid as defined by Equations (4)a–c of A18 and ${\sigma }_{1}(\tau )$, σ₂(τ), and σ₃(τ) are power laws with parameters from Table 1 of A18.

A18 made various cuts to the Casagrande et al. (2011) GCS catalog to ensure high-quality kinematic data and age estimates. The age distribution of their subsample is similar to the full, magnitude-limited GCS, which is known to be biased toward bright F-type stars (Nordström et al. 2004). Since the main-sequence lifetime of a 1.1 M_☉ star is roughly half the age of the universe (Choi et al. 2016), using the full GCS significantly biases kinematic age estimates toward younger ages. In Figure 2, we show age cumulative distributions for the A18 sample compared to a volume-limited sample of the GCS (d < 40 pc). The volume-limited sample is significantly shifted toward older ages.

Zoom In Zoom Out Reset image size

This volume-limited sample still contains a number of F dwarfs with main-sequence lifetimes much shorter than the age of the universe, which biases the sample against older stars. In an attempt to create a sample of stars that better matches the true age distribution of low-mass stars in the solar neighborhood, we further restrict the volume-limited sample to stars with 0.9 M_☉ < M_⋆ < 1 M_☉. This restricts the sample to mainly G dwarfs whose lifetimes are not significantly shorter than the age of the universe and for which the GCS is mostly complete out to 40 pc. The lower limit in mass excludes stars for which the ages are not well constrained. Figure 2 also shows the age distributions for this "unbiased" sample of the GCS and our sample of B16 stars. However, we note that it is extremely difficult to assemble a truly unbiased and complete sample of stars. The B16 sample is comprised of stars originally observed as part of the California Planet Survey (Howard et al. 2010), a radial velocity (RV) exoplanet survey. As such, it is biased against stars with excessive velocity jitter and faint stars. Our log g and T_eff cuts, along with our requirement that each star have Tycho-2 and 2MASS magnitudes as described in the Appendix, introduce additional biases. Overall, however, these biases do not result in any significant age bias for our final sample, as evidenced by the similarity between the age distributions for our B16 sample and the volume-limited "unbiased" GCS sample. This is largely a result of the fact that both samples have been limited to solar-neighborhood Sun-like stars. As these stars have lifetimes on the order of the age of the universe, their age distribution should be very similar to the age distribution of solar-neighborhood M dwarfs.

For consistency, we used our sample of B16 stars to recalibrate the A18 kinematic likelihood. We found the following best-fit relations for the vertex deviation, V component of the solar motion, and dispersion in the three components of the velocity ellipsoid as functions of time:

$$\begin{eqnarray}&&{{\ell }}_{v}=0.406{e}^{-0.163\tau },\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{V}_{\odot }^{{\prime} }=0.314\tau {}^{2}-1.28\tau +17.0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\sigma }_{1}=13.6\tau {}^{0.484},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\sigma }_{2}=7.32\tau {}^{0.493},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\sigma }_{3}=4.20\tau {}^{0.703},\end{eqnarray} \tag{ 7 }$$

where all velocities are in km s⁻¹. We also calculate the U and W components of the solar motion to be 9.33 and 7.95 km s⁻¹, respectively.

We used proper motions and parallaxes from Gaidos et al. (2014) along with radial velocities from R18 to calculate the U, V, and W peculiar velocities for each M dwarf in our exoplanet-host sample. The velocities for each star are listed in Table 2. We used these velocities and Equations (3)–(7) to calculate the posterior probability distributions of the kinematic ages given by Equation (2). We used these posteriors as the prior probability in Equations (1).

Following the approach of A18 for constructing a data-driven likelihood function, we used a sample of FGK stars with measured [Ti/Fe] and stellar ages to calculate the likelihood of a star's measured [Ti/Fe] given an assumed age.

We first divided our FGK sample into 25 age bins spaced so that each bin contains roughly an equal number of stars (∼27 per bin). Motivated by the existence of two chemically distinct populations in the solar neighborhood (e.g., Fuhrmann 1998), we modeled the [Ti/Fe] distribution within a bin as a Gaussian mixture model with two components. We fit for the means, variances, and weights of each component via expectation maximization.³ We assessed the uncertainty in the mixture model parameters by bootstrap resampling the [Ti/Fe] distribution within an age bin and refitting the mixture model, repeating this 10,000 times. In order to create a continuous likelihood as a function of age, we fit an offset power law of the form

$$\begin{eqnarray}&&{\theta }_{i}(\tau )=a+b\tau {}^{c}\end{eqnarray} \tag{ 8 }$$

to each mixture model parameter, θ_i. Here τ is the average age of the stars in the bin. We determined the best-fit parameters via χ² minimization and list them in Table 3. Figure 3 shows our mixture model parameters as a function of age, their uncertainties, and our power-law fits. Note that since the weights of the two components must sum to unity, we only fit to the weights of one component.

Zoom In Zoom Out Reset image size

Figure 4 shows our Gaussian mixture distributions with means, variances, and weights given by Equation (8) with the best-fit parameters from Table 3. In order to use these distributions as likelihood functions for our planet-hosting M dwarf sample, we incorporate the uncertainty in our M dwarf [Ti/Fe] measurements by convolving the Gaussian mixture distributions with a Gaussian kernel with a standard deviation of 0.05 dex.⁴ These distributions represent the empirical probability of measuring a given [Ti/Fe] for a given stellar age, incorporating both the intrinsic scatter in the [Ti/Fe]–age relation and the uncertainty in our M dwarf [Ti/Fe] measurements. The uncertainty in the B16 [Ti/Fe] measurements is much smaller than the intrinsic scatter within an age bin and is inherently included in this scatter. For qualitative comparison, we also show kernel density estimation (KDE) distributions from 100 bootstrap resamplings of the B16 [Ti/Fe] distribution within each age bin. We used a Gaussian kernel with a standard deviation of 0.05 dex. Evaluating the Gaussian mixture model at a measured [Ti/Fe] gives the likelihood as a function of age.

Zoom In Zoom Out Reset image size

For illustrative purposes, we used Equation (4) in Jackson et al. (2008; including the corrected numerical coefficient from Jackson et al. 2009) to calculate the simplified tidal circularization timescales for the planets in our sample. Calculating the tidal circularization timescale requires assuming a value for the tidal dissipation parameter Q, a unitless quality factor that is inversely proportional to tidal dissipation—smaller Q results in stronger dissipation. We assumed a tidal dissipation parameter for the star of Q_⋆ = 10^5.5 and for the planet of Q_p = 10^6.5, as suggested by Jackson et al. (2008), although we note that the tidal Q_p is no doubt different for each planet and likely depends on each planet's mass and structure. In our calculations, we used the planetary masses, semimajor axes, and eccentricities in Table 1. For planetary radii, we used the empirical mass–radius–incident flux relation of Weiss et al. (2013, Equations (8) and (9)). In one case, GJ 436 b, we had a measurement of the planetary radius from transit observations (Maciejewski et al. 2014). We estimated stellar masses and radii from the empirical absolute K-band magnitude relations of Benedict et al. (2016, Equation (11) with coefficients from Table 13) and Mann et al. (2015, Equation (5) with coefficients from Table 1),⁵ respectively. For both relations, we used K-band magnitudes from 2MASS (Skrutskie et al. 2006) and parallaxes from Gaia DR1 (Gaia Collaboration et al. 2016b, 2016a; Lindegren et al. 2016). We accounted for the Gaia zero-point offset; however, it has little effect, as our stars are all nearby with parallaxes of order 100 mas.

Figure 7 shows the eccentricities of the planets in our sample versus the host-star age divided by the tidal circularization timescale. Most planets in our sample have very long timescales, such that it is unlikely that they underwent recent tidal evolution. Two planets, GJ 436 b and GJ 876 d, have tidal circularization timescales shorter than the age of their host star. It has been suggested that GJ 436 b has a massive unseen companion that maintains its moderate eccentricity via Kozai interactions (Beust et al. 2012), so its short tidal circularization timescale may not be so surprising. GJ 876 d has three outer companions that are in a Laplace resonance (Rivera et al. 2010). The innermost planet in the system, GJ 876 d, is not expected to interact with the outer three planets (Trifonov et al. 2018), so its nonzero eccentricity is surprising given its short circularization timescale. We discuss GJ 436 b and GJ 876 d further in Sections 5.2.3 and 5.2.11, respectively.

Zoom In Zoom Out Reset image size

As pointed out by Jackson et al. (2008), the simplified tidal circularization timescale ignores the coupled evolution of eccentricity and semimajor axis and can underestimate the true time to circularize. An alternative approach is to numerically integrate back the tidal evolution equations for both eccentricity and semimajor axis (Equations (1) and (2) of Jackson et al. 2009) from the current age of the star to its formation. Doing so results in the initial eccentricity and semimajor axis of the planet's orbit just after the protoplanetary disk dissipated, assuming no interactions with other bodies in the system. Since we have a posterior for the age of the star, we can determine the probability distribution for the initial eccentricity and semimajor axis of a planet if we assume a Q_p and Q_⋆. To do so, we integrate back the tidal evolution equations 10,000 times, each time drawing a random age from our age posteriors. We show an example of this for GJ 876 d in Figures 8 and 9, where we assume Q_⋆ = 10^5.5 and Q_p = 10^6.5 and 10^5.5, respectively. Assuming a Q_p of 10^6.5 results in a median initial eccentricity and semimajor axis of e_i = 0.7 and a_i = 0.035 au. However, if we assume a Q_p of only 10^5.5, the median initial eccentricity exceeds 1. Therefore, we can constrain the Q_p of GJ 876 d to be greater than ∼10^5.5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Following this approach, we can use our host-star ages to estimate the minimum Q_p possible for each planet by stepping through possible Q_p values until the initial eccentricity exceeds 1. For each planet, we estimate the initial eccentricity probability distribution for 100 Q_p values spaced logarithmically from 10 to 10⁷. We take the maximum Q_p at which the median initial eccentricity is greater than 1 as the minimum possible Q_p for the planet. We hold the stellar tidal dissipation parameter fixed at 10^5.5. For planets that are susceptible to tidal effects around M dwarfs (i.e., on close-in orbits), the effect of the stellar tide is negligible. The minimum Q_p values are shown in Figure 10. Again, this assumes only tidal interactions with the host star, so the results are not necessarily meaningful for dynamically interacting multi-planet systems. We have excluded from the figure planets known to be gravitationally interacting with a companion (GJ 581 b, c, e and GJ 876 b, c, e). To assess the sensitivity of our minimum Q_p estimates to measurement uncertainties, we redid the analysis using planet radii and eccentricities that were smaller by 1σ. For radii, we used the RMSE of the mass–radius relation (1.41 R_⊕ for M < 150 M_⊕, 1.15 R_⊕ for M ≥ 150 M_⊕; Weiss et al. 2013), except for GJ 436 b, where we used the radius uncertainty quoted by Maciejewski et al. (2014). For eccentricities, we used the lower uncertainties listed in Table 1.

Zoom In Zoom Out Reset image size

The tidal dissipation parameter Q_p is poorly constrained even for the planets in the solar system. Gas giants like Jupiter and Saturn have Q_p values somewhere around 10⁵–10⁶ (Goldreich & Soter 1966; Ioannou & Lindzen 1993; Ogilvie & Lin 2004). However, Lainey et al. (2009, 2012, 2017) suggested that Jupiter's and Saturn's tidal Q_p are much lower, around 35,000 and 2500, respectively. Neptune and Uranus both have Q_p values around 10⁴ (Goldreich & Soter 1966; Zhang & Hamilton 2008). Rocky planets have low Q_p values around 10–200 (Goldreich & Soter 1966; Murray & Dermott 2000; Lainey et al. 2007; Henning et al. 2009). Because of the large separation between gas giant and terrestrial Q_p values, they can be used to differentiate between rocky and gaseous planets (Barnes 2015). Since we can only provide a lower bound on Q_p, we cannot place strict constraints on the composition of these planets, though a high minimum Q_p might suggest that the planet is more like a gas giant than a rocky planet. It is also important to note that tidal evolution is strongly dependent on the radius of the planet ($1/{\tau }_{\mathrm{circ}}\propto {R}_{{\rm{p}}}^{5}$). For all but one planet, we assume a radius based on mass and incident flux. Therefore, a high minimum Q_p does not necessarily rule out a rocky composition. Rather, a high minimum Q_p indicates that a planet could be affected by tidal interactions, if its true Q_p is not well above our lower limit and its true radius is not much less than that inferred from empirical mass–radius relations, which are know to have significant scatter (Wolfgang & Lopez 2015; Wolfgang et al. 2016). A small minimum Q_p simply means the planet is insensitive to tidal effects on timescales of the stellar lifetime, e.g., if it obits at a large semimajor axis.

It is also important to note that the tidal evolution equations are only valid under the assumption that the planet's orbital period is shorter than the star's rotational period. While this may not be valid for all systems in this study, it is valid for the most interesting systems, those with short-period planets (P < 10 days) around old stars. Field early-to-mid-M dwarfs typically have rotation periods >10 days (McQuillan et al. 2013), and kinematically old mid-M dwarfs typically rotate with periods >70 days (Irwin et al. 2011; Newton et al. 2016).

GJ 176 hosts a super-Earth in an 8.8 day orbit originally discovered by Forveille et al. (2009). Trifonov et al. (2018) published updated parameters for GJ 176 b, incorporating 23 new CARMENES observations and confirming an $M\sin i\approx 9{M}_{\oplus }$ planet in an ∼8.8 day orbit. They reported a mildly eccentric orbit with $e={0.148}_{-0.036}^{+0.249}$.

Eggen (1998) proposed that GJ 176 is a member of the moving group HR 1614 based on its kinematics. Feltzing & Holmberg (2000) estimated HR 1614 to be about 2 Gyr old. However, our analysis suggests that GJ 176 is an older star with an age of ${8.8}_{-2.8}^{+2.5}\,\mathrm{Gyr}$ and rules out ages less than 2 Gyr at the 3σ level. Furthermore, De Silva et al. (2007) spectroscopically analyzed 18 proposed members of HR 1614 and found the cluster to be metal-rich with log ε_Fe = 7.77 ± 0.033 dex. De Silva et al. (2007) also found that four out of the 18 stars they studied had lower metallicities with log ε_Fe = 7.44–7.55 dex and deviated from the cluster mean abundances in all elements except the n-capture elements, suggesting that they are not members of HR 1614. We measure an iron abundance for GJ 176 of only log ε_Fe = 7.5 ± 0.1 dex, further suggesting that it is not a member of HR 1614.

GJ 176 b has the second-highest minimum Q_p (∼10³) of all single-planet systems in our sample. At a minimum mass of ∼9 M_⊕, GJ 176 b is likely more similar to Uranus and Neptune than to a massive, rocky super-Earth. This is supported by our minimum Q_p for GJ 176 b, which is close to the Q_p of Uranus and Neptune. Orbiting at only 0.066 au, GJ 176 b has likely undergone some tidal evolution. Assuming a radius of 3.5 R_⊕ and a Q_p of 10⁴, we estimate that GJ 176 b started with a high initial eccentricity of e_i ≈ 0.7 and migrated in from an initial semimajor axis of a_i ≈ 0.1 au.

GJ 179 is one of the few M dwarfs known to host a Jupiter analog. GJ 179 b is a Jupiter-mass planet ($M\sin i=0.82{M}_{\mathrm{Jup}}$) in a slightly eccentric (e = 0.21 ± 0.08) 6.3 yr orbit (Howard et al. 2010). Despite a slight Ti enhancement, the kinematics of GJ 179 suggest a moderate age of ${4.6}_{-2.4}^{+3.5}\,\mathrm{Gyr}$ but with a long tail of probability up to older ages.

Orbiting at such a large distance from its host star, GJ 179 b is not expected to be affected by tidal interactions with the host star.

GJ 436 was the second M dwarf found to host an exoplanet (Butler et al. 2004). GJ 436 b has roughly the same radius, mass, and density as Neptune but orbits with a period of only 2.6 days (Maciejewski et al. 2014). Butler et al. (2004) originally estimated that GJ 436 is more than 3 Gyr old based on its kinematics and chromospheric activity. By combining kinematics with Ti enhancement, we constrain the age to be ${8.9}_{-2.1}^{+2.3}\,\mathrm{Gyr}$, making GJ 436 the oldest planet host in our sample (with the potential exception of GJ 625; see Section 5.2.7).

The old age of GJ 436 is surprising considering that it hosts a short-period planet in an eccentric orbit (e = 0.152 ± 0.009; Trifonov et al. 2018). Bourrier et al. (2018) recently reported that the orbit of GJ 436 b is not only eccentric but also nearly perpendicular with the spin axis of the star.

One scenario that has been proposed to explain the eccentricity of GJ 346 b is interaction with a third body (Demory et al. 2007; Maness et al. 2007; Mardling 2008; Ribas et al. 2008). Tong & Zhou (2009) investigated the possible locations of a dynamical companion in either a resonant or nonresonant orbit and argued that the eccentricity of GJ 436 b cannot be maintained by either a nearby or distant companion. Batygin et al. (2009) confirmed that, for most scenarios, the presence of a second planet does not keep GJ 436 b from rapidly circularizing. However, they found that under certain initial conditions where the eccentricities of the two planets are locked at a quasi-stationary point, the eccentricity damping of GJ 436 b can be extended to ∼8 Gyr.

Beust et al. (2012) put forward another hypothesis, suggesting that GJ 436 b originally orbited at a larger semimajor axis and migrated to its current orbit via Kozai migration induced by a distant perturber. In this scenario, the eccentricity damping of GJ 436 b can be delayed by several Gyr. Bourrier et al. (2018) estimated the age of GJ 436 to be ∼5 Gyr based on its rotation period of 44 days and found that Kozai migration could explain both the eccentricity and obliquity of GJ 436 b.

Morley et al. (2017) found that additional interior heat from tidal dissipation is required to explain the observed thermal emission of GJ 436 b. They were able to constrain the tidal Q_p of GJ 436 b to 2 × 10⁵–10⁶. This agrees well with our minimum Q_p estimate of ∼10⁵. With a Q_p  > 10⁵, the tidal dissipation is weak enough that an unseen third body is not required to explain the nonzero eccentricity. However, it does mean that the orbit of GJ 436 b has been significantly altered by tidal effects over its lifetime. Assuming a Q_p = 10⁶, GJ 436 b would have initially orbited with an eccentricity around ∼0.8 at a distance of ∼0.05 au. However, this scenario does not explain the high obliquity of the current orbit reported by Bourrier et al. (2018).

GJ 536 hosts a super-Earth planet in an 8.7 day orbit (Suárez Mascareño et al. 2017a). Using additional RV data from the CARMENES survey, Trifonov et al. (2018) refined GJ 536 b's mass estimate to $M\sin i={6.52}_{-0.40}^{+0.69}{M}_{\oplus }$ and the eccentricity of the orbit to $e={0.119}_{0.032}^{+0.125}$. We estimate the age of GJ 536 to be ${6.9}_{-2.3}^{+2.5}\,\mathrm{Gyr}$.

GJ 536 b is very similar to GJ 176 b. Both orbit at ∼0.066 au with similar eccentricities, 0.15 and 0.12. GJ 176 b has a slightly higher minimum mass than GJ 536 b, 9 M_⊕ versus 6.5 M_⊕. We estimate the minimum Q_p of both planets to be  ∼10³. In the absence of interactions with other unseen planets in the system, GJ 176 b and GJ 536 b are likely similar mini-Neptune planets with extensive gaseous atmospheres. Both are also likely to have undergone tidal circularization and migration. Assuming a Q_p of 10⁴, we estimate that GJ 536 b initially orbited at ∼0.08 au with an eccentricity of ∼0.5.

GJ 581 hosts three bona fide planets: GJ 581 b (Bonfils et al. 2005b), GJ 581 c (Udry et al. 2007), and GJ 581 e (Mayor et al. 2009). The three planets orbit in a very compact configuration with semimajor axes between 0.029 and 0.074 au. Trifonov et al. (2018) used an N-body model to show that all three planets are dynamically interacting and in a stable configuration where each planet's semimajor axis is constant but their eccentricities oscillate on timescales of 50 and 500 yr. Since all three planets are interacting, our minimum Q_p estimates are invalid. Trifonov et al. (2018) showed that this configuration is stable for at least 10 Myr. Our age estimate for GJ 581 of ${6.6}_{-2.5}^{+2.9}\,\mathrm{Gyr}$ suggests that these compact, interacting systems can be stable for several Gyr. This is further supported by the apparent ubiquity of these "compact multiples" (Muirhead et al. 2015).

HD 147379 (GJ 617 A) was the first star discovered to host a planet by the CARMENES survey (Reiners et al. 2018a). GJ 617 A b has a minimum mass of ${M}_{p}\sin i\sim 25{M}_{\oplus }$ and orbits in a nearly circular ∼0.3 au orbit. As such, GJ 617 A b is unlikely to be strongly affected by tidal interactions with its host star.

Vican (2012) estimated the age of GJ 617 A to be ∼1 Gyr based on chromospheric activity and X-ray flux. However, they used the activity–age relations of Mamajek & Hillenbrand (2008), which were only calibrated down to early-K dwarfs (B − V < 0.9 mag), and GJ 617 A is a late-K/early-M with B − V = 1.34. We estimate a slightly older chemo-kinematic age for GJ 617 A of ${5.1}_{-2.4}^{+3.2}\,\mathrm{Gyr}$.

Reiners et al. (2018a) found an additional peak in the periodogram of GJ 617 A corresponding to a period of 21 days, which they attributed to the rotation period of the star. Pepper (2018) confirmed a 22 day rotation period based on 3304 KELT observations (Pepper et al. 2007). Agüeros et al. (2018) recently measured rotation periods for 12 K and M dwarfs in the 1.34 Gyr old cluster NGC 752. They found that late-K dwarfs similar in mass to GJ 617 A rotate with a rotation period of ∼15 days. Assuming a simple Skumanich-like evolution ($\tau \propto {p}_{\mathrm{rot}}^{2}$), a rotation period of 22 days for GJ 617 A would suggest an age of ∼3 Gyr, in rough agreement with our chemo-kinematic estimate.

GJ 625 was only recently discovered to host a super-Earth orbiting at the inner edge of the habitable zone (Suárez Mascareño et al. 2017b). A rocky planet orbiting at such a close distance to its host star is expected to circularize on very short timescales. Indeed, the eccentricity of GJ 625 b is consistent with zero: $e={0.13}_{-0.09}^{+0.12}$. GJ 625 b is not likely to be currently undergoing tidal migration, although it may have in the past.

GJ 625 is a peculiar case where its kinematics strongly favor a young age; however, its low [Fe/H] and high [Ti/Fe] abundances are similar to older thick-disk members. This leads to a combined age posterior that is bimodal and has a significant probability at essentially all ages. The median and ±1σ values of the posterior are ${7.0}_{-4.1}^{+2.7}\,\mathrm{Gyr}$. Estimating the age from the kinematic prior alone yields ${3.9}_{-1.9}^{+3.3}\,\mathrm{Gyr}$, whereas ignoring the kinematic prior and assuming a flat prior results in an age estimate from [Ti/Fe] alone of ${9.5}_{-3.0}^{+2.5}\,\mathrm{Gyr}$.

We can turn to other indications of an M dwarf's age to argue in favor of either the young or old interpretation of the age of GJ 625. A common indicator for the rough age of an M dwarf is its rotation period. Suárez Mascareño et al. (2017b) estimated the rotation period of GJ 625 to be P = 77.8 ± 5.5 days. The relation between rotation period and age for M dwarfs is not well understood. However, studies of young open clusters and field M dwarfs suggest that M dwarfs, like solar-type stars, spin down over time as magnetized stellar winds carry away angular momentum (Irwin et al. 2007, 2011; McQuillan et al. 2013; Newton et al. 2016; Douglas et al. 2017; Rebull et al. 2017). Newton et al. (2016) found that field mid-M dwarfs like GJ 625 show a bimodal distribution in rotation period with peaks at ∼1 and ∼100 days. The evolutionary link between these two populations is not clear; however, Newton et al. (2016) found that M dwarfs with rotation periods >70 days are kinematically consistent with an old population with an average age of 5 Gyr. GJ 625's slow ∼80 day rotation is similar to that of the slowest rotating (and presumably oldest) stars in the field of similar mass (Newton et al. 2017) and is therefore more consistent with the older peak of our age posterior.

Montes et al. (2001) listed GJ 625 as a possible member of the young Ursa Major moving group, which would imply an age of only ∼0.5 Gyr (although Montes et al. 2001 did list it with the caveat that it fails both their peculiar velocity and radial velocity criteria). We can rule out membership based on GJ 625's long rotation period, low activity (${\mathrm{log}}_{10}({R}_{{HK}}^{{\prime} })=-5.5\pm 0.2$; Suárez Mascareño et al. 2017b), and chemical dissimilarity (Tabernero et al. 2017). We also note that the BANYAN Σ web tool⁶ lists a 0% probability that GJ 625 is a member of the Ursa Major moving group and a 99.9% probability that it is a field star (Gagné et al. 2018).

Wright et al. (2016) used archival HARPS spectra to discover three potentially rocky planets around GJ 628 (Wolf 1061), with one planet, GJ 628 c, orbiting within the habitable zone. Astudillo-Defru et al. (2017) ruled out a third planet orbiting at 67 days as originally proposed by Wright et al. (2016) but found significant evidence for a third planet orbiting at 217 days.

GJ 628 is one of a couple cases where a solar [Ti/Fe] estimate results in a very broad likelihood and the posterior is dominated by the kinematic prior that favors younger ages. Based on this, we estimate an age for GJ 628 of ${4.3}_{-2.0}^{+3.1}\,\mathrm{Gyr}$. However, Astudillo-Defru et al. (2017) claimed a rotation period of 95 days for GJ 628, which is supported by the photometric monitoring presented in Kane et al. (2017). Such a long rotation period suggests an age >5 Gyr, as discussed in the previous section. If a strict gyrochronological relation does exist for M dwarfs, this long rotation period is inconsistent with our age estimate.

We estimate a minimum Q_p ∼ 10³ for the innermost planet, GJ 628 b. Such a high minimum Q_p would suggest that the ${M}_{p}\sin i\sim 2{M}_{\oplus }$ planet is not rocky but in fact has a large gaseous envelope. However, GJ 628 b and GJ 628 c have eccentricities that are consistent with zero within the measurement error. If we assume their radii and eccentricities are overestimated by 1σ, we can no longer constrain the minimum Q_p to be greater than 10. Therefore, given that they are both likely rocky and orbit within 0.1 au of their host star, they likely were tidally circularized not long after their primordial disk dissipated.

Montes et al. (2001) listed GJ 628 as a member of the young (125 Myr; Stauffer et al. 1998) Pleiades moving group. However, the BANYAN Σ web tool lists a 0% probability that GJ 628 is a member of the Pleiades moving group and a 99.9% probability that it is a field star.

GJ 649 hosts one known planet with a minimum mass similar to Saturn, ${M}_{p}\sin i\sim 100{M}_{\oplus }$ (Johnson et al. 2010b). With a semimajor axis of 1.1 au, the orbit of GJ 649 b is not expected to be influenced by tidal effects.

GJ 649 is another case where a near-solar [Ti/Fe] does little to constrain the age of the star—other than ruling out the oldest ages—and the kinematics favor younger ages. We estimate an age of ${4.5}_{-2.0}^{+3.0}\,\mathrm{Gyr}$ for GJ 649.

GJ 849 hosts a roughly Jupiter-mass planet in a nearly circular, ∼5 yr orbit (Butler et al. 2006). Orbiting at over 2 au, GJ 849 b is not expected to undergo any tidal circularization or migration. Butler et al. (2006) noted a linear trend in their RV time-series data, suggesting a possible second planet in the system on an even longer orbit. With RV measurements spanning 17 yr, Feng et al. (2015) found strong evidence for a second planet with $M\sin i\approx {M}_{J}$ orbiting with a period of 15.1 ± 1.1 yr.

Like GJ 628 and GJ 649, GJ 849 has nearly solar [Ti/Fe] and kinematics that skew the age posterior to younger ages. Based on this, we estimate the age of GJ 849 to be ${4.9}_{-2.1}^{+3.0}\,\mathrm{Gyr}$.

The planetary system around GJ 876 is a benchmark system for studying the formation and migration of planets in compact systems. GJ 876 hosts a total of four known planets, three of which are in a Laplace 1:2:4 mean-motion resonance (Rivera et al. 2010). The resonance is chaotic but expected to be stable on timescales of at least 1 Gyr (Rivera et al. 2010; Martí et al. 2013). The old age we infer for GJ 876 (${8.4}_{-2.0}^{+2.2}$ Gyr) suggests that such chaotic resonances can be stable for several Gyr, assuming the planets migrated into this configuration soon after their formation (Batygin et al. 2015).

Even though the innermost planet, GJ 876 d, is not in the resonant chain with the other three planets, our old age for GJ 876 has some interesting implications for its past migration. Originally, Rivera et al. (2010) estimated that planet d orbited with an unusually high eccentricity (e = 0.207 ± 0.055) for a planet that orbits at only 0.02 au with a period of about 2 days. Recently, Millholland et al. (2018) and Trifonov et al. (2018) reanalyzed the system by fitting dynamical N-body models to existing and new RV data. Both analyses revised the eccentricity of planet d down to values more consistent with a circular orbit. Trifonov et al. (2018) reported a best-fit eccentricity of $e={0.082}_{-0.025}^{+0.043}$, while Millholland et al. (2018) reported a best-fit eccentricity of e = 0.057 ± 0.039. Using the Trifonov et al. (2018) estimate for the eccentricity, we constrain the minimum Q_p of GJ 876 d to be >5 × 10⁵. Such a high Q_p is surprising for an ∼7 M_⊕ planet. Our minimum Q_p for GJ 876 d is an order of magnitude higher than the Q_p of Uranus and Neptune and three orders of magnitude greater than the Q_p of terrestrial bodies.

A few scenarios could explain the nonzero eccentricity of GJ 876 d. One explanation is that the planet has a peculiar structure that is very inefficient at dissipating tidal energy. Another explanation is that GJ 876 d originally orbited with a larger semimajor axis where it gravitationally interacted with the outer three planets and, at some point in its history, fell out of a chaotic resonance with the outer planets and migrated in. A third scenario is that there is another, unseen planet in the system that is gravitationally interacting with GJ 876 d.

We note that the eccentricity of GJ 876 d has undergone downward revision recently. If additional observations in future studies lead to further downward revision and the conclusion that the orbit is essentially circular, it may not be necessary to invoke one of the above scenarios to explain the orbit of GJ 876 d. In that case, our minimum Q_p will be overestimated. However, given its proximity to its host star and our age estimate of ${8.4}_{-2.0}^{+2.2}\,\mathrm{Gyr}$ for the system, GJ 876 d likely has undergone significant tidal circularization and migration in its lifetime.

Montes et al. (2001) listed GJ 876 as a member of the young Pleiades moving group, which is inconsistent with our old age estimate for the star. The BANYAN Σ web tool gives a 0% probability that GJ 876 is a member of the Pleiades moving group and an 81.7% probability that it is a field star. Interestingly, it also gives an 18.3% probability that GJ 876 is a member of the Beta Pictoris moving group.