A link between rocky planet composition and stellar age

ABSTRACT

Interior compositions are key for our understanding of Earth-like exoplanets. The composition of the core can influence the presence of a magnetic dynamo and the strength of gravity on the planetary surface, both of which heavily impact thermal and possible biological processes and thus the habitability for life and its evolution on the planet. However, detailed measurements of the planetary interiors are extremely challenging for small exoplanets, and existing data suggest a wide diversity in planet compositions. Hitherto, only certain photospheric chemical abundances of the host stars have been considered as tracers to explain the diversity of exoplanet compositions. Here, we present a homogeneous analysis of stars hosting rocky exoplanets, with ages between 2 and 14 Gyr, revealing a correlation between rocky exoplanet compositions and the ages of the planetary systems. Denser rocky planets are found around younger stars. This suggests that the compositional diversity of rocky exoplanets can be linked to the ages of their host stars. We interpret this to be a result of chemical evolution of stars in the Milky Way, which modifies the material out of which stars and planets form. The results imply that rocky planets that form today, at similar Galactocentric radii, may have different formation conditions, and thus different properties than planets that formed several billion years ago, such as the Earth.

exoplanets, planets and satellites: formation, stars: solar-type, Galaxy: evolution

1 INTRODUCTION

The internal composition of exoplanets is a key characteristic to understand their diversity and formation. The interior of small, Earth-like planets can be approximated as having a dense metallic core, a lower and upper silicate mantle, and a thin layer of ice or liquid water (Valencia, Sasselov & O’Connell 2007). The amount of elements such as iron (Fe), silicon (Si), or magnesium (Mg) plays an important role in influencing the extent of each of these layers (Sotin, Grasset & Mocquet 2007). In the Solar system, the relative abundances of Fe, Si, and Mg are similar for the Sun and Earth, Mars, and Venus (McDonough & Yoshizaki 2021). The properties of the stars around which planets form are expected to influence planetary interiors and determine the planetary core size (Dorn et al. 2015). These properties in turn directly influence characteristics such as the presence or absence of a magnetic field, which is driven by the dynamo of the molten iron core (Taylor 1963), the importance of outgassing to create a secondary atmosphere through plate tectonics (Müller, Mather & Dutkiewicz 2022), the presence of a carbon cycle that may affect habitability (Alibert 2014), and the strength of gravity on the planet as its size changes for a given mass (Yang, Komacek & Abbot 2019).

Understanding the interior compositions of rocky exoplanets is challenging. A limited number of small planets with known atmospheres have been observed with space-based spectroscopy (e.g. Kreidberg et al. 2014; Kempton et al. 2023; Mikal-Evans et al. 2023). Several small rocky planets have been observed by JWST, but thus far, no atmospheres have been detected (e.g. Lim et al. 2023; Zieba et al. 2023).

When such atmospheric measurements are not possible or an atmosphere is not present, the determination of exoplanet compositions relies on measurements of a planet’s mass and radius, which then yields its bulk density. These measurements can then be used to derive parametric mass–radius relationships (e.g. Weiss & Marcy 2014; Wolfgang, Rogers & Ford 2016; Bashi et al. 2017; Chen & Kipping 2017; Otegi, Bouchy & Helled 2020) or infer compositions through theoretical mass–radius tracks that solve the equation of states of the varying elemental constituents (e.g. Seager et al. 2007; Sotin et al. 2007; Zeng & Seager 2008; Zeng, Sasselov & Jacobsen 2016).

The interpretation of mass–radius measurements for planets that have both a core and an atmosphere, as for planets with a significant hydrogen–helium (H–He) envelope, is often highly degenerate (Selsis et al. 2007; Lopez & Fortney 2014; Kite et al. 2020; Bean, Raymond & Owen 2021). Planets without a thick atmosphere are more suited for studies of planetary core compositions. As super-Earths are expected to be stripped cores, their interior composition models are generally more uniquely defined (Valencia et al. 2007; Rogers 2015). This makes it feasible to infer the extent of the dense Fe-rich core for such planets through an estimate of the planet mean density. One way to select these planets is by exploiting the so-called radius valley, a dearth in occurrence of planets with radii of about 1.8 |$\mathrm{ R}_{\rm{\oplus }}$| (Fulton et al. 2017). The valley is a likely consequence of atmospheric mass-loss (Owen & Wu 2013; Gupta & Schlichting 2019) and separates planets with a significant H–He atmosphere (sub-Neptunes) from those that have lost their atmosphere (super-Earths). Super-Earths sit below the valley at sizes smaller than about two Earth radii, though the exact location of the valley is dependent on orbital period (Van Eylen et al. 2018) and stellar type (e.g. Fulton & Petigura 2018; Ho & Van Eylen 2023).

To measure the mean density of a planet, both its radius and its mass need to be estimated. For the majority of known exoplanets, radii and masses are inferred by combining transit and radial velocity observations (e.g. Lam et al. 2021; Van Eylen et al. 2021). Both quantities are measured relative to the radii and masses of the host star, which are often determined heterogeneously. Individual planet discovery and characterization papers use different atmospheric parameters, stellar spectra models, stellar evolution models, and fitting codes, in order to calculate stellar parameters.

Here, we homogenize the parameters for a sample of stars with planets that have masses and radii measured by the transit and radial velocity methods, and we use these new parameters to refine the properties of the planets. We then link the properties of these planets to the properties of the host stars, and discuss an observed correlation between stellar age and rocky planet composition, which we interpret in the context of Galactic chemical evolution.

In Section 2, we outline how our sample was selected and how we homogeneously characterized the host stars. We also validate the stellar parameters. In Section 3, we describe how we use the homogenized stellar parameters to determine new planet parameters. In Section 4, we link properties of the planets to their star, and outline an observed correlation between stellar age and rocky planet composition. We discuss this finding in the context of Galactic chemical evolution in Section 5. In Section 6, we present our conclusions.

2 HOMOGENEOUS HOST STAR CHARACTERIZATION

2.1 Sample selection

We start by selecting all confirmed planets with |$R_{\rm {p}} < 4 \,\mathrm{ R}_{\rm{\oplus }}$| from the NASA Exoplanet Archive¹ (Akeson et al. 2013). We also require that the planets have both radius and mass measurements from transit and radial velocity observations. This means that values in both rv_flag and trans_flag must be 1.

The DR2 (Data Release 2) Gaia ID in the NASA Exoplanet Archive is matched to the corresponding (DR3 Data Release 3) Gaia ID using the `DR2_neighbours' table in the Gaia archive, following the algorithm in Marrese et al. (2019). Using python modules astropy (Astropy Collaboration 2013) and astroquery (Ginsburg et al. 2019), the sample is then cross-matched to the Gaia DR3 Astrophysical parameters supplementary materials (Ulla et al. 2022). Here, stars with a renormalized unit-weighted error (RUWE) of |$>1.4$|⁠, or parallax over error |$< 5$| are removed as these targets are likely to be either binaries (Belokurov et al. 2020) or have large uncertainties on their heliocentric distance, respectively. Lester et al. (2021) note that there is a decreased efficiency of planet detection around binary stars, for planets with radii less than 2 R|$_{\rm{\oplus }}$|⁠. We note that 55 Cnc is a known long-period binary system (Bourrier et al. 2018); however, its Earth-like planet, 55 Cnc e, has been observed and analysed in many planet characterization and follow-up papers over the last 15 yr (see e.g. Winn et al. 2011; Demory et al. 2016; Dai et al. 2019). Kepler-20 and Kepler-21 have potentially detected comparisons from more in-depth astrometry Mugrauer (2019) and imaging Ginski et al. (2016), respectively. However, our stellar parameters agree well with those in the literature, as these systems are extremely well-characterized with asteroseismology. We require that each target has a value reported for teff_gspspec, mh_gspspec, and their associated upper and lower confidence intervals. Further, we remove stars with teff_gspspec  |$< 5000$| K, for which model ages cannot be reliably measured (Berger et al. 2020). Some stars have reported non-zero quality flags in the first 13 flag parameters, described in Recio-Blanco et al. (2023). We confirm that no flag is greater than a value of 2 for these first 13 flags.

Following these selection criteria, we are left with a sample of 89 stars hosting 116 small planets with measured masses and radii. They all have |$T_{\mbox{eff}}$| and |$[\text{M/H}]$| from GSP-Spec, Gaia G, |$RP$|⁠, and |$BP$| photometry, and parallax, |$\varpi$|⁠, from Gaia DR3, in addition to radial velocity K amplitudes and transit-derived |$\frac{R_{\mathrm p}}{R_\star }$| values from the NASA Exoplanet Archive. The values and references of the latter two quantities are noted in Table 2.

Table 1.

BASTA output stellar parameters for host stars of the 26 rocky planets in our sample – the sample for which the composition-age trend is observed.

	Planet name	Stellar mass	Stellar radius	|$T_{\mbox{eff}}$|	|$[\text{Fe/H}]$|
		(M|$_{\odot }$|⁠)	(R|$_{\odot }$|⁠)	(K)	(dex)
1	Kepler-20 b	|$0.89^{+ 0.07}_{- 0.06}$|	|$0.86^{+ 0.03}_{- 0.03}$|	|$5700^{+ 46}_{- 95}$|	|$-0.09^{+ 0.17}_{- 0.14}$|
2	EPIC 249 893 012 b	|$1.09^{+ 0.09}_{- 0.07}$|	|$1.62^{+ 0.06}_{- 0.07}$|	|$5586^{+ 79}_{- 78}$|	|$0.20^{+ 0.15}_{- 0.12}$|
3	K2-291 b	|$0.90^{+ 0.07}_{- 0.05}$|	|$0.86^{+ 0.03}_{- 0.03}$|	|$5615^{+ 62}_{- 63}$|	|$-0.03^{+ 0.14}_{- 0.13}$|
4	WASP-47 e	|$1.06^{+ 0.09}_{- 0.08}$|	|$1.10^{+ 0.04}_{- 0.05}$|	|$5542^{+ 73}_{- 73}$|	|$0.33^{+ 0.12}_{- 0.08}$|
5	EPIC 220 674 823 b	|$0.93^{+ 0.09}_{- 0.10}$|	|$1.01^{+ 0.04}_{- 0.04}$|	|$5403^{+ 86}_{- 87}$|	|$0.16^{+ 0.16}_{- 0.16}$|
6	HD 80 653 b	|$1.17^{+ 0.07}_{- 0.05}$|	|$1.18^{+ 0.05}_{- 0.05}$|	|$5904^{+ 84}_{- 86}$|	|$0.33^{+ 0.09}_{- 0.11}$|
7	CoRoT-7 b	|$0.86^{+ 0.05}_{- 0.05}$|	|$0.81^{+ 0.03}_{- 0.03}$|	|$5374^{+ 45}_{- 72}$|	|$0.05^{+ 0.15}_{- 0.13}$|
8	Kepler-107 c	|$1.19^{+ 0.08}_{- 0.06}$|	|$1.36^{+ 0.09}_{- 0.08}$|	|$5939^{+ 166}_{175}$|	|$0.34^{+ 0.10}_{- 0.13}$|
9	Kepler-107 b	|$1.19^{+ 0.08}_{- 0.06}$|	|$1.36^{+ 0.09}_{- 0.08}$|	|$5939^{+ 166}_{175}$|	|$0.34^{+ 0.10}_{- 0.13}$|
10	K2-131 b	|$0.78^{+ 0.05}_{- 0.06}$|	|$0.76^{+ 0.03}_{- 0.03}$|	|$5055^{+ 72}_{- 87}$|	|$-0.03^{+ 0.14}_{- 0.17}$|
11	HD 3167 b	|$0.84^{+ 0.07}_{- 0.07}$|	|$0.87^{+ 0.03}_{- 0.04}$|	|$5316^{+ 81}_{- 68}$|	|$-0.01^{+ 0.14}_{- 0.16}$|
12	K2-111 b	|$0.80^{+ 0.06}_{- 0.08}$|	|$1.2^{+ 0.04}_{- 0.05}$|	|$5841^{+ 101}_{- 70 }$|	|$-0.54^{+ 0.19}_{- 0.21}$|
13	Kepler-100 b	|$1.08^{+ 0.05}_{- 0.04}$|	|$1.49^{+ 0.03}_{- 0.02}$|	|$5964^{+ 25}_{- 39}$|	|$-0.00^{+ 0.12}_{- 0.09}$|
14	TOI-561 b	|$0.77^{+ 0.04}_{- 0.06}$|	|$0.82^{+ 0.03}_{- 0.03}$|	|$5415^{+ 67}_{- 79}$|	|$-0.27^{+ 0.14}_{- 0.14}$|
15	Kepler-323 c	|$0.90^{+ 0.13}_{- 0.17}$|	|$1.18^{+ 0.05}_{- 0.06}$|	|$5805^{+ 10}_{- 86}$|	|$-0.18^{+ 0.39}_{- 0.34}$|
16	K2-38 b	|$1.04^{+ 0.09}_{- 0.09}$|	|$1.05^{+ 0.04}_{- 0.04}$|	|$5680^{+ 75}_{- 80}$|	|$0.21^{+ 0.16}_{- 0.14}$|
17	HD 20 329 b	|$0.88^{+ 0.08}_{- 0.10}$|	|$1.10^{+ 0.04}_{- 0.04}$|	|$5607^{+ 78}_{- 85}$|	|$-0.09^{+ 0.17}_{- 0.17}$|
18	Kepler-21 b	|$1.30^{+ 0.10}_{- 0.07}$|	|$1.86^{+ 0.05}_{- 0.05}$|	|$6204^{+ 77}_{- 63}$|	|$0.13^{+ 0.11}_{- 0.14}$|
19	55 Cnc e	|$0.94^{+ 0.08}_{- 0.08}$|	|$0.94^{+ 0.04}_{- 0.03}$|	|$5279^{+ 62}_{- 82}$|	|$0.26^{+ 0.15}_{- 0.12}$|
20	Kepler-93 b	|$0.90^{+ 0.04}_{- 0.05}$|	|$0.92^{+ 0.01}_{- 0.02}$|	|$5674^{+ 57}_{- 58}$|	|$-0.12^{+ 0.11}_{- 0.14}$|
21	Kepler-10 b	|$0.91^{+ 0.09}_{- 0.11 }$|	|$1.05^{+ 0.04}_{- 0.04}$|	|$5716 ^{+ 88}_{- 84}$|	|$-0.09^{+ 0.18}_{- 0.18}$|
22	KOI-94 b	|$1.21^{+ 0.09}_{- 0.07}$|	|$1.34^{+ 0.09}_{- 0.08}$|	|$6232^{ 161}_{- 203}$|	|$0.2^{+ 0.14}_{- 0.17}$|
23	Kepler-406 b	|$1.02^{+ 0.07}_{- 0.06}$|	|$0.98^{+ 0.04}_{- 0.03}$|	|$5893^{+ 57}_{- 102}$|	|$0.08^{+ 0.16}_{- 0.14}$|
24	TOI-1444 b	|$0.89^{+ 0.08}_{- 0.08}$|	|$0.91^{+ 0.03}_{- 0.03}$|	|$5417^{+ 74}_{- 61}$|	|$0.06^{+ 0.17}_{- 0.14}$|
25	HD 213 885 b	|$1.00^{+ 0.09}_{- 0.09}$|	|$1.06^{+ 0.04}_{- 0.04}$|	|$5985^{+ 82}_{- 75 }$|	|$-0.07^{+ 0.18}_{- 0.16 }$|
26	HD 137 496 b	|$1.05^{+ 0.10}_ {-0.09}$|	|$1.43^{+ 0.05}_{- 0.05}$|	|$5915^{+ 85}_{- 81 }$|	|$-0.04^{+ 0.18}_{- 0.16}$|

Table 2.

RV semi-amplitude (K), used to re-calculate planet masses, and transit-derived planet-to-star ratio (⁠|${R_{\mathrm p}}/{R_\star }$|⁠), used to re-calculate planet radius, with references for each as follows: A – Bonomo et al. (2023), B – Hidalgo et al. (2020), C – Burke et al. (2014), D – Kosiarek et al. (2019), E – Bryant & Bayliss (2022), F – Guenther et al. (2017), G – Frustagli et al. (2020), H – John, Collier Cameron & Wilson (2022), I – Adams et al. (2021), J – Vanderburg et al. (2016), K – Fridlund et al. (2017), L – Marcy et al. (2014), M – Brinkman et al. (2023), N – Burke et al. (2014), O – Toledo-Padrón et al. (2020), P – Murgas et al. (2022), Q – López-Morales et al. (2016), R – Weiss et al. (2013), S – Bourrier et al. (2018), T – Dressing et al. (2015), U – Dai et al. (2019), V – Dai et al. (2021), W – Espinoza et al. (2020), X – Azevedo Silva et al. (2022).

	Planet name	K (m s−1)	Ref.	|$R_{\mathrm p}/R_{\star }$|	Ref.
1	Kepler-20 b	|$4.23^{+0.54}_{-0.54}$|	A	|$0.01925^{+0.00018}_{-0.00018}$|	C
2	EPIC 249 893 012 b	|$3.55^{+0.43}_{-0.43}$|	B	|$0.01040^{+0.00040}_{-0.00040}$|	B
3	K2-291 b	|$3.33^{+0.59}_{-0.59}$|	D	|$0.01614^{+0.00062}_{-0.00033}$|	D
4	WASP-47 e	|$4.55^{+0.37}_{-0.37}$|	E	|$0.01458^{+0.00013}_{ -0.00013}$|	E
5	EPIC 220 674 823 b	|$6.50^{+0.52}_{ -0.52}$|	A	|$0.01601^{+0.00031}_{ -0.00029}$|	F
6	HD 80 653 b	|$3.62^{+0.21}_{-0.21}$|	A	|$0.01210^{+ 0.00040}_{-0.00040}$|	G
7	CoRoT-7 b	|$3.10^{+0.70}_{-0.70}$|	H	|$0.01870^{+0.00030}_{-0.00030}$|	H
8	Kepler-107 c	|$3.29^{+0.66}_{-0.66}$|	A	|$0.01012^{+0.00013}_{-0.00013}$|	A
9	Kepler-107 b	|$1.44^{+0.68}_{-0.68}$|	A	|$0.00973^{+0.00013}_{-0.00013}$|	A
10	K2-131 b	|$8.00^{+1.30}_{-1.30}$|	A	|$0.01720^{+0.00041}_{-0.00027}$|	I
11	HD 3167 b	|$3.56^{+0.15}_{-0.15}$|	A	|$0.01725^{+0.00023}_{ -0.00023}$|	J
12	K2-111 b	|$2.30^{+0.30}_{-0.30}$|	A	|$0.01255^{+0.00050}_{ -0.00048}$|	K
13	Kepler-100 b	|$1.90^{+0.80}_{-0.80}$|	L	0.00809	L
14	TOI-561 b	|$2.18^{+0.20}_{-0.20}$|	M	|$0.01550^{+0.001}_{-0.001}$|	M
15	Kepler-323 c	|$2.80^{+1.40}_{-1.40}$|	A	|$0.01337^{+0.00023}_{-0.00023}$|	N
16	K2-38 b	|$3.02^{+0.43}_{-0.43}$|	A	|$0.01330^{+0.00071}_{-0.00076}$|	O
17	HD 20 329 b	|$5.07^{+0.41}_{-0.42}$|	P	|$0.01390^{+0.00050}_{-0.00050}$|	P
18	Kepler-21 b	|$2.70^{+0.46}_{-0.46}$|	A	|$0.00788^{+0.00005}_{-0.00005}$|	Q
19	55 Cnc e	|$6.02^{+0.24}_{-0.23}$|	S	|$0.01820^{+0.00020}_{-0.00020}$|	S
20	Kepler-93 b	|$1.89^{+0.21}_{-0.21}$|	A	|$0.01475^{+0.00006}_{-0.00006}$|	T
21	Kepler-10 b	|$2.34^{+0.21}_{-0.21}$|	A	|$0.01268^{+0.00004}_{-0.00004}$|	U
22	KOI-94 b	|$3.30^{+1.40}_{-1.40 }$|	R	|$0.01031^{+0.00014}_{-0.00014}$|	R
23	Kepler-406 b	|$2.89^{+0.60}_{-0.60}$|	L	0.01229	L
24	TOI-1444 b	|$3.30^{+0.58}_{-0.59}$|	V	|$0.01410^{+0.00060}_{-0.00058}$|	V
25	HD 213 885 b	|$5.30^{+0.39}_{-0.39}$|	W	|$0.01453^{+0.00041}_{-0.00042}$|	W
26	HD 137 496 b	|$2.14^{+0.29}_{-0.29}$|	X	|$0.0076^{+0.00030}_{-0.00029}$|	X

2.2 Sample preparation

We apply calibration steps to several Gaia DR3 data products. The effective temperature selected for use in this study is the ‘GSP-Spec’ temperature, which is determined using the Matisse–Gauguin algorithm detailed in Creevey et al. (2023) and Recio-Blanco et al. (2016), to derive atmospheric parameters from R  |$\sim$| 11 000 spectra observed with the RVS instrument on Gaia (Cropper et al. 2018). The lower and upper uncertainties are determined as the 16th and 84th percentile using Matisse–Gauguin and Monte Carlo realization. As the distributions of the atmospheric parameters for all targets are overall symmetric and Gaussian, we simplify the description of them and symmetrize the uncertainty by using the mean of the two uncertainties. The reported uncertainties are small (⁠|$< 30$|  |$\mathrm{K}$| in some cases) and do not take into account systematic errors. Following Vallenari et al. (2022), we add a systematic uncertainty of 90 K in quadrature to the reported uncertainty, consistent with the expected uncertainties in effective temperature scales (Tayar et al. 2022).

A similar approach was taken with the errors for the metallicity, with a median scatter reported by Recio-Blanco et al. (2023) of 0.13 dex added in quadrature to the symmetrized error. We also corrected the metallicity values using the polynomial calibration over |$\log g$| as described in Recio-Blanco et al. (2023).

Finally, we corrected the measured parallax values for the known parallax zero-point offset in Gaia DR3 using the approach outlined in Lindegren et al. (2018, 2021).

2.3 Determination of stellar parameters: BASTA

We make use of the BAyesian STellar Algorithm code (BASTA; Aguirre Børsen-Koch et al. 2022) to perform a homogeneous analysis of exoplanet host stars, alongside a custom computed set of stellar evolution models. Given a pre-computed grid of stellar models, BASTA uses a Bayesian approach to compute the posterior distribution of a given stellar parameter using a set of observational constraints. It is fast and flexible and allows for a homogeneous characterization even when the same quality or quantity of data is not available for all targets in a given sample.

The Bayesian framework used is, briefly, as follows. If we treat |$\Theta$| as the desired parameters, mass, radius, and age, and D as the input data, e.g. |$T_{\mbox{eff}}$|⁠, |$\log g$|⁠, and |$[\text{Fe/H}]$|⁠, then the likelihood of finding |$\Theta$| given D is

$$\begin{eqnarray} P(\boldsymbol{\Theta} | \boldsymbol{D}) = \frac{P(\boldsymbol{D}| \boldsymbol{\Theta})P(\boldsymbol{\Theta})}{P(\boldsymbol{D})}. \end{eqnarray}$$

(1)

The full likelihood is the product of the likelihood of groups of observables |$D_i$|⁠:

$$\begin{eqnarray} P(\boldsymbol{D}| \boldsymbol{\Theta}) = \prod _i P (\boldsymbol{D_i}| \boldsymbol{\Theta}), \end{eqnarray}$$

(2)

where

$$\begin{eqnarray} P(\boldsymbol{D_i}| \boldsymbol{\Theta}) = \frac{1}{\sqrt{2\pi |{\boldsymbol{ C}_i}|}} \exp (-\chi ^2/2). \end{eqnarray}$$

(3)

Here, |$|{\boldsymbol{ C}_i}|$| is the determinant of the covariance matrix, which is used to account for known correlations in parameters in the fitting procedure. The computed posterior is used to calculate the marginalized posterior for the stellar model parameter |$\theta$|⁠, where |$\Theta ^{\prime }$| is all model parameters excluding |$\theta$|⁠, and |$x_{\boldsymbol{\varTheta }}$| is a bayesian weight used to account for the volume of parameter space a model characterized by |$\Theta ^{\prime }$| encompasses (Aguirre Børsen-Koch et al. 2022).

$$\begin{eqnarray} P(\theta | {\boldsymbol{ D}}) = \int P(\theta , \boldsymbol{\varTheta ^{\prime }}| \boldsymbol{D})x_{\boldsymbol{\varTheta }} \mathrm{ d}\boldsymbol{\varTheta ^{\prime }}. \end{eqnarray}$$

(4)

Quantities that are not directly observable like stellar age can then be inferred given the constraints from the observables.

We determine the stellar properties using grid-based stellar modelling, comparing the observed parameters to the predicted theoretical quantities. Specifically, we use the astrometric coordinates (right ascension, declination), parallax (⁠|$\varpi$|⁠), metallicity (⁠|$[\text{M/H}]$|⁠), effective temperature (⁠|$T_{\mbox{eff}}$|⁠), and Gaia magnitudes: G, BP, and RP, all from Gaia DR3, as inputs.

For stars with literature stellar masses |$M_{\star } < 1.15$|  |${\mathrm M}_\odot $|⁠, we build a grid with |$\sim 6000$| evolutionary tracks of stellar models. Our grids of stellar models were computed with the Garching Stellar Evolution Code (garstec; Weiss & Schlattl 2008). The code utilizes a combination of the equation of state by the opal group (Rogers, Swenson & Iglesias 1996; Rogers & Nayfonov 2002), and the Mihalas–Hummer–Däppen equation of state (Daeppen et al. 1988; Hummer & Mihalas 1988; Mihalas, Dappen & Hummer 1988; Mihalas et al. 1990). We use the OPAL opacities (Rogers & Iglesias 1992; Iglesias & Rogers 1996) at high temperatures supplemented by the opacities of (Ferguson et al. 2005) at low temperatures. garstec uses the NACRE nuclear reaction rates (Angulo et al. 1999) except for |$^{14}$|N(⁠|$p,\gamma)^{15}$|O and |$^{12}$|C(⁠|$\alpha ,\gamma)^{16}$|O for which the rates from Formicola et al. (2004) and Hammer et al. (2005) were used. The stellar models are computed using the Asplund et al. (2009) solar mixture.

Convection in the models is parametrized using mixing-length theory (Böhm-Vitense 1958; Kippenhahn, Weigert & Weiss 2012), where the mixing length parameter is allowed to vary in the range 1.6–2.0. We sample the parameter space by utilizing the Sobol quasi-random, low-discrepancy sequences to uniformly populate the parameter space (Sobol & Levithan 1976; Sobol 1977; Antonov & Saleev 1980; Fox 1986; Bratley & Fox 1988; Joe & Kuo 2003). The stellar grid samples initial iron abundances from |$-1.0$| to 0.6 dex, stellar masses from 0.7 to 1.2 |${\mathrm M}_\odot $|⁠, mixing-length parameter from 1.6 to 2.0, and initial helium fraction from 0.24 to 0.32. The stellar models include the effects of atomic diffusion. We compute the theoretical asteroseismic mode frequencies for the models using the Aarhus adiabatic oscillation package (adipls; Christensen-Dalsgaard 2008). For the computation of synthetic magnitudes, we use the bolometric corrections of Hidalgo et al. (2018) and we use the three-dimensional dust map from Green et al. (2019) for computing the extinction along the line of sight for each target. BASTA allows for prior probability distributions to be taken into account when computing the posterior distributions. We used the Salpeter initial mass function (Salpeter 1955) as a prior to quantify the expected mass distribution of stars favouring low-mass stars as the most abundant.

garstec models are calculated along fixed iron abundances |$[\text{Fe/H}]$| with alpha-abundance element mixtures ranging from |$[\rm{\alpha /Fe}] =-0.2$| to |$[\rm{\alpha /Fe}] =0.6$|⁠. BASTA adopts the following conversion to the bulk metallicity |$[\text{M/H}]$|⁠:

$$\begin{eqnarray} [\text{Fe/H}] = [\text{M/H}]-\mbox{corr}([\rm{\alpha /Fe}]), \end{eqnarray}$$

(5)

where corr(⁠|$[\rm{\alpha /Fe}]$|⁠) is a correction factor determined as a function of the adopted solar mixture for the selected model, as in Salaris, Chieffi & Straniero (1993). We supply |$[\text{M/H}]$| input values from the Gaia Radial Velocity Spectrometer (RVS) spectra , which are then converted by BASTA to |$[\text{Fe/H}]$| using equation (5) for each |$[\rm{\alpha /Fe}]$| grid, therefore ensuring consistency in the atmospheric input parameters.

For stars with literature |$M_{\star } >1.15$|  |${\mathrm M}_\odot $|⁠, we build a similar grid. However, for stars with a mass limit of |$1.1\,\mathrm{ M}_\odot < M_{\star } < 1.6$|  |$\mathrm{ M}_\odot $| we do not include atomic diffusion but instead include the effects of core overshooting in the models. With this more appropriate input physics, we avoid underestimating masses and metallicities when fitting for more massive stars. Some stars in our sample have been precisely characterized using asteroseismology (Huber et al. 2013; Silva Aguirre et al. 2015). In order to ensure consistency, we do not use literature analyses of these stars but instead include the measured global asteroseismic observables, the large frequency separation |$\Delta \nu$| and the frequency of maximum power |$\nu _\text{max}$|⁠, in our inference of stellar properties for these targets. These additional constraints improve the precision on the stellar age for this subset of our sample. (Soderblom 2010). We adopted asteroseismic parameters for the following stars: Kepler-21, Kepler-25, Kepler-68, Kepler-93, Kepler-100, Kepler-103, Kepler-129 – three of which are in the final rocky planet sample, and four of which are sub-Neptunes (see Section 3). Input asteroseismic values were adopted from Huber et al. (2013). For Kepler-10, previous stellar characterization using asteroseismology has been carried out (Fogtmann-Schulz et al. 2014), but no |$\nu _\text{max}$| has been reported. Hence, we did not include asteroseismic constraints of Kepler-10.

Stellar parameters for our final sample are listed in Table 1. In Fig. 1, the stars are shown in a Hertzsprung–Russell diagram, demonstrating how comparing homogeneous constraints on effective temperature, radius, and metallicity from Gaia to stellar models enables precise relative measurements of ages for solar-type stars, which span from about 2 to 14 Gyr for our sample. These tests demonstrate the accuracy and precision of |$\approx$|2 Gyr (⁠|$\approx$|30 per cent for the median age in our sample) of ages in our sample. The precision is also consistent with our median estimated age uncertainties.

Figure 1.

Homogeneous data from Gaia enable the precise relative measurement of ages for planet-hosting stars. Stellar radius versus effective temperature for the 25 host stars in our sample that host rocky exoplanets. Points are colour-coded by the inferred stellar ages from BASTA. Grey points connected via thin lines indicate previous, heterogeneously determined parameters for these host stars from the literature.

We have tested our stellar mass, radius, and age estimates in a various ways. We summarize our validation tests in the appendices. In Appendix A1, we validate stellar mass and radius against literature values, and by using alternative spectroscopic inputs. Our stellar ages are validated in Appendix A2, where we compare results from our method to those of asteroseismology, gyrochronology, kinematics, and detailed abundances. We also apply this validation to a gold-standard sample of Kepler stars studied with asteroseismology.

3 PLANET SAMPLE AND PARAMETERS

Once the new, homogeneous set of stellar parameters were obtained, and validated, we recalculated a new value of radius and mass for each planet. Planet radius is measured indirectly, by using the depth of the transit and the stellar radius as follows:

$$\begin{eqnarray} \Delta F \approx \delta \approx \left(\frac{R_{\mathrm{ p}}}{R_{\star }}\right)^2, \end{eqnarray}$$

(6)

where |$R_\mathrm{ p}$| and |$R_{\star }$| are planet and star radius, respectively, and |$\Delta F$| is the observed reduction in relative flux, which is equal to transit depth |$\delta$|⁠. We use the reported values of |$\frac{R_{\mathrm p}}{R_{\star }}$|⁠, which have taken into account limb darkening effects, and our updated value of |$R_{\star }$| in order to calculate a new |$R_{\mathrm p}$|⁠.

Planet mass is also recalculated using the new homogeneous stellar mass. We follow the same method as described in Fulton et al. (2018). In short, we substitute our new value for |$M_{\star }$| into

$$\begin{eqnarray} M_\mathrm{ p}\sin i = \frac{K}{K_0} \sqrt{1-e^2} \frac{M_{\star }}{{\rm M}_{\odot }} \frac{2}{3} P^{\frac{1}{3}}, \end{eqnarray}$$

(7)

where |$M_\mathrm{ p}$| denotes planet mass, K the radial velocity semi-amplitude, i the inclination, e the eccentricity of the orbit, |$M_{\star }$| the stellar mass, |${\rm M}_{\odot }$| the solar mass, and P the orbital period. These values, except for |$M_{\star }$|⁠, are extracted from the NASA Exoplanet Archive.

In Table 2, we list the adopted transit and RV observables used to update our planet parameters, in addition to the sources for these values. In 15/26 (⁠|$\sim 57~{{\ \rm per\ cent}}$|⁠) of cases, K and |$\frac{R_{\mathrm p}}{R_{\star }}$| are determined in the same paper.

In the case of Corot-7 b, several literature studies have arrived at significantly differing values for the RV amplitude. We opt to use the values reported by John et al. (2022), as this is the most recent study and it applies state-of-the-art stellar activity mitigation in its RV analyses. In all plots, we use values from this paper as the literature values.

At this stage, we use the inferred planet radius measurements to determine which planets are super-Earths located below the radius valley, and limit all further analysis to these planets only. To do so, we use the location of the radius valley as determined by Ho & Van Eylen (2023). Specifically, the location of the radius valley depends on orbital period and stellar mass, as |$\log (R_{\mathrm p}/{\mathrm R}_{\oplus }) = -0.09 \log (P/\mathrm{d}) + 0.21 \log (M_{\star }/{\rm M}_{\odot }) + 0.35$| (equation 11, Ho & Van Eylen 2023). For each planet in our sample, we calculate the location of the valley at its orbital period and stellar mass, and keep those planets whose radius is smaller than the valley radius. Finally, we exclude three planets that have only reported mass upper limits, and thus no or uninformative mass errors. We exclude Kepler-131c, as its mass is non-physically large and poorly constrained as previous studies have noted (Weiss & Marcy 2014). Kepler-406c is also excluded, due to having an unphysical density, calculated based on a tentative RV curve in Weiss & Marcy (2014). The K2-32 system is a close-in binary (Mugrauer 2019) whose Gaia parameters, in particular |$[\text{M/H}]$| are clearly very discrepant with the literature values, so this star is removed from the analysis.

4 A LINK BETWEEN AGE AND COMPOSITION

In Table 3, we list our adopted planet parameters. The planets in our sample are shown in a mass–radius diagram in Fig. 2. Models for rocky planets with varying compositions (taken from Zeng et al. 2016) are also shown. These tracks are created from semi-empirical relations derived under isothermal formation scenarios for rocky exoplanets with gaseous envelopes of negligible mass fraction. They represent rocky planets at different percentages of mantle materials such as MgSiO|$_3$| and MgSiO|$_4$| with respect to iron, which is assumed to be the main constituent of a core (Seager et al. 2007). In Fig. 2, we furthermore colour-code each planet by its age, where we have assumed that planet ages are identical to the age of their host stars, since studies of young stellar objects indicate that the protoplanetary disk phase lasts only a few million years following the collapse of their parent molecular cloud (Hartmann et al. 1998; Haisch, Lada & Lada 2001). Visually, this seems to reveal a correlation between planet composition and stellar age. In the age range of the sample, younger stars appear more likely to host denser, more iron-rich rocky planets. Planets above the ‘Earth like’ composition track in Fig. 2 seem to orbit about older stars than those below it. We observe a hitherto unknown trend that rocky planets that are denser and have a higher inferred Fe content appear to orbit stars that are younger, suggesting a link between the composition of rocky planets and their age. We now investigate this trend in more detail.

Figure 2.

Rocky planet density and composition depends on the age of the host star. The mass and radius of small planets, colour-coded by age, alongside composition tracks from Zeng et al. (2016). The upper line corresponds to a planet with composition of 100  per cent MgSiO|$_3$|⁠, i.e. 100  per cent rocky mantle. The middle, line corresponds to a planet with 67.5  per cent MgSiO|$_3$|⁠, and 32.5  per cent Fe, providing an analogue for Earth-like composition. The lower line represents planets that would be composed of 100  per cent Fe, essentially pure iron cores. Faint, thin lines connected to the points indicate movement from previously reported, heterogeneous values from the NASA Exoplanet Archive.

Table 3.

Parameters of small, rocky exoplanets used in Figs 2 and 3. From left to right the table lists planet name, radius, mass, and age, and their associated errors. The quoted lower (Age|$_{\rm {err,m}}$|⁠) and upper (Age|$_{\rm {err,p}}$|⁠) age uncertainties are the 16th and 84th percentiles directly extracted from the Bayesian posterior distribution function (PDF) by BASTA. For each planet, the references for our adopted values of |$\frac{R_p}{R_\star }$| and radial velocity K-amplitude are given in Table 2.

	Planet name	R	R|$_{\rm {err}}$|	M	M|$_{\rm {err}}$|	Age	Age|$_{\rm {err,m}}$|	Age|$_{\rm {err,p}}$|
		(R|$_{{\oplus }}$|⁠)	(R|$_{{\oplus }}$|⁠)	(M|$_{{\oplus }}$|⁠)	(M|$_{{\oplus }}$|⁠)	(Gyr)	(Gyr)	(Gyr)
1	Kepler-20 b	1.80	0.08	9.46	1.32	3.04	2.33	3.76
2	EPIC 249 893 012 b	1.83	0.12	9.01	1.21	9.20	1.07	1.05
3	K2-291 b	1.52	0.11	6.35	1.18	3.46	2.37	4.65
4	WASP-47 e	1.74	0.07	6.86	0.71	6.43	3.52	4.23
5	EPIC 220 674 823 b	1.77	0.09	8.03	0.93	11.07	5.56	5.70
6	HD 80 653 b	1.56	0.09	5.70	0.41	3.48	2.26	1.97
7	CoRoT-7 b	1.53	0.14	5.83	0.68	3.42	2.57	3.63
1	Kepler-107 c	1.50	0.10	9.81	2.03	5.06	2.36	2.54
2	Kepler-107 b	1.44	0.10	3.71	1.76	5.06	2.36	2.54
10	K2-131 b	1.43	0.09	7.63	1.31	6.52	4.34	7.22
11	HD 3167 b	1.64	0.15	4.92	0.39	10.62	6.21	5.22
12	K2-111 b	1.76	0.12	5.39	0.80	14.03	2.43	2.53
13	Kepler-100 b	1.31	0.06	5.92	2.50	6.72	0.93	0.83
14	TOI-561 b	1.35	0.07	2.19	0.23	12.19	6.43	4.30
15	Kepler-323 c	1.73	0.09	6.19	3.21	11.25	5.36	5.00
16	K2-38 b	1.52	0.10	7.68	1.21	4.59	3.06	3.92
17	HD 20 329 b	1.67	0.10	7.31	0.85	13.50	4.62	3.39
18	Kepler-21 b	1.60	0.05	7.14	1.27	4.28	0.61	1.2
19	55 Cnc e	1.86	0.08	8.20	0.63	7.26	4.63	6.26
20	Kepler-93 b	1.48	0.03	4.64	0.55	6.21	1.41	1.83
21	Kepler-10 b	1.45	0.06	3.25	0.41	9.84	4.68	4.57
22	KOI-94 b	1.50	0.22	8.80	3.76	3.32	1.96	1.90
23	Kepler-406 b	1.32	0.06	6.22	1.33	1.76	1.21	2.81
24	TOI-1444 b	1.41	0.09	3.74	0.72	8.52	5.17	6.55
25	HD 213 885 b	1.68	0.08	8.44	0.88	4.56	2.73	3.22
26	HD 137 496 b	1.19	0.07	4.08	0.63	7.57	1.81	1.82

4.1 Planet density and inferred iron fraction

In Fig. 2, we showed a mass–radius diagram with planets colour-coded by age. From this diagram, a visual correlation between planet composition and age appears to be apparent. An alternative way of showing this effect is to plot the stellar density as a function of age, which is what we do in Fig. 3, which similarly reveals a visual correlation between planet density and age.

Figure 3.

Planet density, in log scaling on the y-axis, against stellar age on the x-axis, for planets and stars in our rocky planet sample. A clear trend can be seen, with denser planets orbiting around younger stars.

To quantify this trend, we first note that planet density is an incomplete metric for measuring planet composition. This is because planets with the same density do not necessarily have the same composition, as planet composition tracks do not follow equal-density tracks. To address this, we can infer the planet composition using a grid of planet interior models (Zeng et al. 2019, Table 2 therein), which gives values for radii as a function of mass, at incremental values of iron per cent. This provides a sparse grid of points for which values of mass, radius, and iron mass fraction are available. We take this grid and linearly interpolate in two dimensions, namely in mass–radius space, using `scipy.interpolate.LinearNDInterpolator'. The newly interpolated grid that we create is shown in Fig. 4. In order to assign a value of iron fraction to each planet in our sample, we draw 100 random samples from a Gaussian distribution with |$\mu$|⁠, |$\sigma$| treated as the measurement and uncertainty respectively, for planet mass and radius. For each random tuple, a value of |$\rm{\%Fe}$| is calculated from the interpolation, and these 100 values are used to create a posterior distribution function for the |$\rm{\%Fe}$|⁠. We take the median of |$\rm{\%Fe}$| as the final value, and lower and upper uncertainties as the 16th and 84th percentiles, respectively. Planets with best-fit masses and radii that fall above the ‘pure-rock’ line in are assigned a value of |$\rm{\%Fe} = 0$|⁠. While this may not be a realistic physical description of the planet compositions, we note that all such planets are consistent with pure-rock models within their uncertainties, and that some of this discrepancy may be due to systematic errors in the models.

Figure 4.

Interpolated grid of iron mass fraction values for planets at a given mass and radius, measured in Earth-scaled values. The colour bar represents the per cent mass fraction of iron, as interpolated from table 2 of Zeng et al. (2016). White overlayed points highlight our sample of rocky planets.

4.2 Significance of planet composition-age trend

In Fig. 5, we show the inferred |$\rm{\%Fe}$| of the planets (see Section 4.1) as a function of stellar age (⁠|$\tau$|⁠). To quantify this trend, we first calculate correlation coefficients. We find Pearson, Spearman, and linear least squares regression coefficients that are |$-0.62$|⁠, |$-0.63$|⁠, |$-0.62$|⁠, suggesting a negative correlation where the planet iron percentage decreases with stellar age.

Figure 5.

Younger, denser rocky planets have higher iron content. Calculated iron mass percentages of planets, based on their measured masses and radii and models interpolated from a grid by Zeng et al. (2019), shown as a function of age. A Bayesian linear regression yields a slope and intercept of |$\alpha = -8.0^{+4.6}_{-5.0}$|⁠, |$\beta = 100^{+34}_{-32}$|⁠, and a median correlation coefficient of |$-0.62^{+0.23}_{-0.17}$|⁠. The solid line shows the best-fitting model, and the shaded area is the 1 − |$\sigma$| confidence interval. The residuals of the model are shown in the bottom panel.

In order to further quantify the observed trend, we fitted the observed age–planet iron trend using a linear model |$\tau = \alpha \rm{\%Fe} + \beta$|⁠. We did so using linmix (Kelly 2007), which performs a Bayesian linear regression, taking into account the errors on both x and y (i.e. age and iron content). We find |$\alpha = -8.0^{+4.6}_{-5.0}$|⁠, |$\beta = 100^{+34}_{-32}$|⁠, and a median correlation coefficient of |$-0.62^{+0.23}_{-0.17}$|⁠. Fig. 6 shows the posterior distributions for |$\alpha$|⁠, |$\beta$|⁠, and the correlation coefficient. The correlation is highly significant.

Figure 6.

The two leftmost panels show posteriors of the intercept and slope – defined as |$\alpha$| and |$\beta$| in our model, of the fit between age and planet iron mass fraction |$\rm{\%Fe}$| using linmix (Kelly 2007). In the rightmost panel, a PDF of the Pearson correlation coefficient is plotted, showing a clear negative linear correlation between age and |$\rm{\%Fe}$|⁠. These values were calculated using a Bayesian linear regression fit.

To explore alternative methods, we fit our data using an orthogonal distance regression method (Boggs & Rogers 1989), using scipy.odr. This method takes into account errors on age and %Fe, and can fit to non-linear functional forms. We heuristically choose two models, shown respectively in Fig. 7. First, we fit a linear model of the form |$\alpha x + \beta$|⁠, which is plotted in dark blue, atop the shaded |$1-\sigma$| confidence region in light blue. With this method, |$\alpha = -5.52 \pm 1.49$| and |$\beta = 81.0 \pm 12.2$|⁠, which agree, within errors, to our values from the Bayesian linear regression discussed above, and also suggest a non-zero slope at a significance of more than |$2 \sigma$|⁠. We repeat this method using an inverse linear model of the form |$(\gamma /x) + \eta$|⁠, giving values of |$\gamma = 174 \pm 55$| and |$\eta = 9.35 \pm 10.5$|⁠. Both methods show correlation between the two parameters, confirming the results derived using linmix in Fig. 5.

Figure 7.

Calculated iron mass percentages of planets, based on their measured masses and radii and models interpolated from a grid by Zeng et al. (2019), shown as a function of age. Upper: An orthogonal distance regression using a linear model of the form |$\alpha x + \beta$| is plotted in a dark solid line, atop the shaded |$1-\sigma$| confidence region in light shading. Lower: An ODR using an inverse linear model of the form |$(\gamma /x) + \eta$| is plotted in a dark solid line, atop the shaded |$1-\sigma$| confidence region in shading. The Solar system planets are plotted, with Earth, Venus, and Mars having similar |$\rm{\%Fe}$|⁠, between 25 and 40. Mercury has 70 per cent Fe. Solar system values used are from McDonough & Yoshizaki 2021), and the residuals of both models are shown in the respective bottom panels.

4.3 Investigating alternative correlations

At a given stellar radius, temperature, and metallicity, a higher-mass star will have a lower age than a lower-mass star. Therefore, we investigate whether the observed trend with stellar age could be caused by a physical effect that is related to stellar mass.

Fig. 8 shows stellar age versus mass for the stars in our sample, colour-coded by their inferred planet |$\rm{\%Fe}$|⁠. This plot shows that higher-mass stars span a narrower range of ages and are mostly younger relative to their lower-mass counterparts, as expected. Visually, a trend of |$\rm{\%Fe}$| with age can be seen more clearly than a trend with stellar mass. An age trend can be seen most clearly at |$M_\mathrm{ s} < 1\, {\rm M}_{\odot }$|⁠, while the data are consistent with such a trend at higher masses too.

Figure 8.

Stellar mass versus age for stars in our rocky planet sample. The spread of ages is larger for lower-mass stars, which is expected as such stars live longer. Data points are colour-coded by |$\rm{\%Fe}$|⁠, demonstrating a visual gradient that is stronger with age (from left to right) than with mass (from bottom to top).

To further explore the effect of stellar mass, in Fig. 9 we show the planet |$\rm{\%Fe}$| as a function of age, in two bins of stellar mass, split at |$M_s=0.95$|  |${\mathrm M}_\odot $| to contain a roughly equal part of our sample. The trend between stellar age and |$\rm{\%Fe}$| can be clearly seen within the lower mass bin, and less clearly in the higher mass bin, which has a narrower age dispersion. The points are color-coded by metallicity, again showing a lack of correlation between stellar metallicity and |$\rm{\%Fe}$|⁠. We note that the dearth of high density (high |$\rm{\%Fe}$|⁠) planets around old, high-mass stars may be caused by detection bias, since some of these stars will have appreciably evolved off the main-sequence.

Figure 9.

|$\rm{\%Fe}$| plotted as a function of stellar age, for two subsets of our sample: stars with |$M_{\star } < 0.95$|  |${\mathrm M}_\odot $| (left), and |$M_{\star } >0.95$|  |${\mathrm M}_\odot $| (right). The age trend can be seen clearly in the left panel. The trend is weaker in the right-hand panel, which contains fewer stars of older ages. The points are color-coded by stellar metallicity. While on average less massive stars are older and slightly metal poorer than more massive stars, we do not observe a clear trend between planet |$\rm{\%Fe}$| and stellar metallicity.

We finally plot planet density as a function of stellar mass in Fig. 10. A weak correlation is observable, with a Pearson correlation coefficient |$r=0.35$|⁠, with P value of .08, as expected due to the inverse correlation between stellar mass and age. This is weaker than the trend observed between age and density. This may be expected as on average, more massive stars are younger (see Fig. 8). However, as not all old stars are less massive, this may explain why the observed correlation with stellar density appears weaker than that observed with stellar age. The latter is more strongly correlated with alpha-abundances |$[\rm{\alpha /Fe}]$|⁠, which change with Galactic evolution, and which we expect to have a more measurable impact on planet composition (i.e. Santos et al. 2017).

Figure 10.

Stellar mass as a function of planet density. A weak correlation is seen, as expected due to the inverse correlation between stellar mass and age. This is weaker than the trend observed between age and density.

We also plot stellar age as a function of |$[\text{Fe/H}]$|⁠. Studies of Galactic stellar populations have shown that the age–metallicity relation is mostly flat until about 10 Gyr (e.g. Haisch et al. 2001; Snaith et al. 2015). This is also what we observe, as we see no significant correlation between metallicity and age except at the highest ages, above approximately 10 Gyr, as shown in Fig. 11. This explains why we do not see a strong correlation between planet composition and stellar |$[\text{Fe/H}]$|⁠. The observed planet density–age trend seems to instead show the greatest rate of change at around 6–8 Gyr.

Figure 11.

Stellar metallicity |$[\text{Fe/H}]$| as a function of stellar age. No significant correlation between age and metallicity is observed, except for the highest ages, in agreement with literature.

This is the first time that a link between the composition of rocky planets and their age is observed. Despite large uncertainties in age, this is the first homogeneous set of ages for a sample of such planets, using the same stellar models, input parameters, and astrometry. Therefore, we interpret age as an indirect link to terrestrial planet density as a result of its direct link to the chemical composition of the host star.

4.4 Potential sample biases

We consider whether the observed planet density versus age correlation could be affected by observational biases. In Fig. 12, we plot a range of parameters as a function of stellar age, such as planet radius, planet mass, orbital period, stellar mass, stellar radius, and stellar magnitude. Although the stellar parameters are calculated homogeneously in this work, a homogeneous re-calculation of planet RV semi-amplitudes and transit depths is beyond its scope. Such heterogeneity may create biases in planet parameter space. Many different facilities and missions are used for these measurements; thus, we expect that such biases, should they exist, to be stochastic, and unlikely to systematically affect the observed trend.

Figure 12.

Planetary and stellar parameters as a function of age for the 25 stars and 26 planets in this sample. Pearson correlation coefficients r are quoted in each panel to indicate the calculated monotonic relationship between the variables.

4.4.1 Sample selection: transit detection

First, we consider whether the sample selection could be biased. As most planets in our sample were first discovered by transit, we first consider whether transit detectability could influence our results. From the top left panel of Fig. 12, it can be observed that the planet size distribution does not appear to depend on age, suggesting planet transits can be detected across the entire age range within our sample. While there may be reasons to expect that small planets are more often missed around young stars because such stars are more active, such an effect is expected to be unimportant at main sequence (⁠|$\ge$|1 Gyr) ages.

We conclude that it is unlikely that transit detection has a significant impact on our conclusions.

4.4.2 Sample selection: RV observations

To be included in our sample, transiting planets require mass measurements, obtained using RV observations. Similar to with transit observations, it is a priori plausible that lower-mass planets are missed around the youngest stars, as these stars may be more active. Such an effect is not observed in the top right panel of Fig. 12. While stellar rotation or activity are a function of age, they primarily affect achievable RV precision for very young ages (⁠|$< 500$| Myr), which are not present in our sample. The effect on RV precision across main sequence ages is small. For example, Isaacson et al. (2024) suggest that stars with values for the chromospheric activity indicator |$R^{\prime }_{HK} < 4.8$| are considered inactive. As shown in Huber et al. (2022), |$R^{\prime }_{HK}$| is generally below such values for stars above about 2 Gyr. With the exception of Kepler-406 (which hosts a low-mass planet), all stars in our sample are older than this. Nevertheless, there may be a decreasing sensitivity to lower density planets towards younger ages (e.g. the lower-left region in Fig. 5), as higher stellar activity levels affect the mass determination more than the transit detection in such systems. This may affect the shape or strength of the observed trend.

4.4.3 Other effects

We consider whether other effects could influence our results. In the bottom right panel of Fig. 12, we show the observed sample as a function of age and stellar magnitude. No clear correlation is observed, suggesting again that our sample of planets with both transits and RV observations is primarily a magnitude-limited sample that is not selected as a function of age (which would have largely been unknown to observers at the time of detection or follow-up). A correlation between stellar mass and age can be observed as expected (the middle right panel of Fig. 12), since more massive stars evolve faster than less massive stars. We find no evidence that the observed correlation with planet properties is an effect of stellar mass rather than age (see previous section, see also the top panels and the middle-left panel of Fig. 12). We further look at age as a function of orbital period (the middle-left panel of Fig. 12). While many planets in our sample orbit at short periods, others have periods of several days and no particular correlation as a function of age can be observed. Overall we conclude that no particular sample bias can be identified that would lead to the observed effect. The most plausible bias would be an absence of low-mass planets orbiting young stars due to the challenge of RV observations. We find no evidence of this occurring, but also note that such an effect could only explain the absence of low-density planets around young stars, and not the similar absence of high-density planets around older stars.

5 THE GALACTIC CHEMICAL CONTEXT

Galactic chemical evolution predicts a link between stellar chemical abundances and ages. Older stars are expected to be enhanced in alpha-elements like O, Mg, Si. Younger stars are more metal-rich and less alpha-enhanced, as a result of enrichment of the interstellar medium by type Ia supernovae. Observations have qualitatively confirmed this picture (Tucci Maia et al. 2016; Ciucă et al. 2021; Viscasillas Vázquez et al. 2022).

Previous studies have predicted the link between the composition of rocky planets to the chemical composition of their host stars (e.g. Santos et al. 2017). This was observed, as a mass–period–metallicity relation (Sousa et al. 2019), and later when the star iron mass fraction, defined as [Fe/Fe + Si  + Mg], was shown to correlate to a modelled planet iron mass fraction, and to the planet bulk density (Adibekyan et al. 2021). Abundance trends in refractory and volatile elements of solar twins have also been suggested to be a tracer of the formation of Earth-like planets (Meléndez et al. 2009; Ramírez, Meléndez & Asplund 2009). Further, it is well understood that abundance ratios themselves, such as C/O and Mg/Si, directly influence planet properties, such as the size of the planetary core (Unterborn et al. 2014; Dorn et al. 2015; Unterborn & Panero 2017). On the other hand, recent re-analysis of planet radius and mass measurements suggest that the link between stellar Mg and Si abundances with rocky planet composition is significantly weaker than previously thought (Brinkman et al. 2024).

A potential link between stellar composition and rocky planet composition has separately been suggested, due to recent work on polluted white dwarfs. Hollands et al. (2021) found that old, lithium-rich white dwarf spectra showed signs of crust-dominated planet debris accretion. Elms et al. (2022) presented analysis of the oldest-known polluted white dwarf. With the use of future Gaia releases, a larger sample of polluted white dwarfs can be used to verify the trend we see in this work, by comparing the composition of pollutants from the spectra, and their ages, which are related to cooling time-scales.

Our results demonstrate for the first time, via homogeneous host star characterization, that these stellar compositional links with rocky exoplanet compositions are also linked to age. It is well known that the elements considered to be influential on rocky planet formation in this way, are strong tracers of stellar age (see e.g. the Si/age relation in Snaith et al. 2014). Our results therefore support robust theoretical predictions and suggestions that Galactic chemical evolution is the driving factor for the compositional diversity of rocky exoplanets (O’Neill, Lowman & Wasiliev 2020; Cabral et al. 2023).

We also note that the [Mg/Si] ratio itself varies with age (Bedell et al. 2018), at a different rate than the alpha elements individually, due to the enhanced production of Si in comparison to Mg, by type Ia supernovae. In simulations, this ratio has been shown to directly impact the mass of a planet at a given radius (Unterborn & Panero 2017), in agreement with our results.

Fig. 5 shows the correlation between stellar age and rocky planet composition. It may further be interpreted that younger, more iron-rich stars, are able to form planets of a greater compositional diversity. Older (⁠|$>8$| Gyr), less iron-rich stars may not have been able to produce planets with e.g. |$>50 \rm{\%Fe}$|⁠. This is also seen in Fig. 9, where stars with |$M_{\mathrm{ s}} > 0.95$|  |$M_{\star }$| do not appear to host any alpha-rich, less dense rocky planets with |$\rm{\%Fe} < 20$|⁠. This may be due to the fact that higher mass stars typically belong to a sample of stars that were born more recently. Due to chemical evolution, they are more iron-rich, and less alpha-enhanced. Thus, they may be more likely to be able to form a dense, iron-rich rocky planet due to this chemical abundance.

Fig. 9 also demonstrates that [Fe/H] does not play a significant role in the observed trend, as for each mass bin [Fe/H] is relatively scattered. This is discussed further in section 5.1 of Brinkman et al. (2024).

In comparing the left and right panels of Fig. 9, we tentatively observe that lower-mass stars appear to have fewer of the highest-density planets. However, we caution this effect could be caused by the relatively small sample size and more data are required to investigate this further. Finally, in comparing the left and right panels of Fig. 9, we also see that there are fewer planets orbiting higher mass stars at higher ages in our sample. We expect this is a natural consequence of the faster evolution of higher mass stars, which are on average younger. While not unexpected, this complicates disentangling the effect of stellar mass and stellar age on planet properties, and a larger sample is required to fully investigate these effects.

6 CONCLUSIONS

We find that planets that formed more recently (about |$2\!-\!6$| Gyr old), around younger stars, have a higher Fe content than planets that formed earlier (⁠|$6\!-\!14$| Gyr). Our findings imply that the composition of rocky planets in our solar neighbourhood is dependent on when, in the history of our Galaxy, these planets formed. We interpret this as an effect of Galactic chemical evolution. More metal-rich stars born more recently, may be more likely to host denser rocky planets, due to the higher abundance of iron in their discs.

Although beyond the scope of this work, the relationship between rocky planet density and age for systems outside of the limitations of this sample could be considered. This may include planets orbiting M and early-K dwarfs (with |$T_{\mbox{eff}} < 5000$| K), which could be aged using alternative techniques such as gyrochronology (e.g. Engle & Guinan 2023) or kinematics (e.g. Sagear et al. 2024). This could also be extended to younger (⁠|$< 2$| Gyr) systems, when mass measurements of the planets become available. Probing such data would allow for a more thorough analysis of the source(s) of the observed correlation.

Furthermore, the PLATO mission (Rauer et al. 2014) is predicted to significantly increase the sample size of known small, rocky exoplanets. Alongside planet discovery, its goals include homogeneous stellar characterization, including the use of asteroseismology to provide precise ages for planet host stars, which should allow for more data to be added to plots such as Fig. 5, potentially allowing for more in-depth investigation to the physical causes of the correlation. These denser, more Fe-rich planetary cores are important for planetary processes such as plate tectonics (Shah et al. 2022; Unterborn et al. 2022), life-protecting magnetic dynamos (Lundin, Lammer & Ribas 2007) and affect the planetary surface gravity. The results of this work suggest that age, in addition to stellar composition and mass, is an important parameter to consider when interpreting the composition of rocky planets, their formation, and their potential habitability.

ACKNOWLEDGEMENTS

AW would like to thank the Science and Technology Facilities Council (STFC) for funding support through a PhD studentship. VVE and DK have been supported by UK’s Science & Technology Facilities Council through the STFC grants ST/W001136/1 and ST/S000216/1. DH acknowledges support from the Alfred P. Sloan Foundation, the National Aeronautics and Space Administration (80NSSC19K0597, 80NSSC21K0652), and the Australian Research Council (FT200100871). DK was supported by MWGaiaDN, a Horizon Europe Marie Skłodowska-Curie Actions Doctoral Network funded under grant agreement no. 101072454 and also funded by UK Research and Innovation (EP/X031756/1). AS has been supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (CartographY; grant agreement ID 804752). Funding for the Stellar Astrophysics Centre was provided by The Danish National Research Foundation (grant agreement no. DNRF106). PP acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon Europe research and innovation programme (grant agreement no. 101076489).

DATA AVAILABILITY

This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/Gaia. The Gaia archive website is https://archives.esac.esa.int/Gaia. All data used for analyses in the paper are available in Tables 1–3.

This study makes use of numpy (https://numpy.org), scipy (https://scipy.org), matplotlib (Hunter 2007), astropy (Astropy Collaboration 2013), astroquery (Ginsburg et al. 2019), BASTA (https://basta.readthedocs.io), and linmix (https://linmix.readthedocs.io).

Footnotes

1

https://exoplanetarchive.ipac.caltech.edu/index.html, last accessed 15/11/2023.

2

https://docs.galpy.org/en/v1.9.2/#acknowledging-galpy

APPENDIX A

A1 Validation of stellar masses and radii

We validate our stellar masses and radii by using input atmospheric parameters from the SWEET-Cat catalogue (Santos et al. 2013; Sousa et al. 2021). SWEET-Cat is a catalogue of spectroscopic parameters for exoplanet host stars, extracted from high-precision spectra and fitted to model atmospheres. These input parameters are not fully homogeneous: about 25 per cent of the SWEET-Cat sample has spectroscopic parameters extracted completely homogeneously. We use the same grid and formalism as described earlier in Section 2.3, but replace the input spectroscopic parameters metallicity [M/H] and effective temperature |$T_{\mbox{eff}}$| from GSP-Spec pipeline with |$[\text{Fe/H}]$| and |$T_{\mbox{eff}}$| from SWEET-Cat. Fig. A1 shows radius and mass as determined by our original method and inputs as well as the results found using the same method but input parameters from SWEET-Cat instead. We see excellent agreement with a median residual scatter of 4.2 per cent and 4.1 per cent, respectively, validating that input parameters from Gaia DR3 GSP-Spec measurements are appropriate for this method of stellar characterization. Due to their homogeneity, we adopt the Gaia DR3 GSP-Spec measurements for our final stellar characterization.

Figure A1.

Masses and radii for the 25 stars in our sample, as determined by using our BASTA stellar properties inference method using to different sets of spectroscopic inputs: Gaia DR3 GSP-Spec (y-axis) and SWEET-Cat (x-axis). We see excellent agreement, with median residual scatters of 4.2 per cent and 4.1 per cent.

We further compare the stellar masses and radii derived here to literature values. In Fig. A2, we show the stellar masses and radii derived here, relative to the literature values currently adopted for those stars in the NASA Exoplanet Archive. Despite the heterogeneity of the comparison values, we again find excellent agreement, with a median residual scatter of 4.7  per cent and 3.8  per cent.

Figure A2.

Masses and radii for the 25 stars in our sample, as determined by using our BASTA fitting method using Gaia DR3 GSP-Spec spectroscopic parameters (y-axis), plotted as a function of literature parameters. We see excellent agreement, with median residual scatter of 4.7  per cent and 3.8  per cent.

A2 Validation of stellar ages

We validate our stellar ages by comparing our approach to a range of other methods. In the next subsections, we discuss our ages in the context of asteroseismology, gyrochronology and chemical abundances, and kinematics.

A2.1 Asteroseismology

First, we compared BASTA output ages with and without including global asteroseismic parameters (⁠|$\Delta \nu$|⁠, |$\nu _\text{max}$|⁠). The results are shown in Fig. A3 (left). We find a standard deviation for the absolute residuals of 0.9 Gyr. We observe a median residual offset and standard deviation of |$15 \pm 15~{{\ \rm per\ cent}}$| and 40  per cent, showing good agreement with our estimated uncertainties.

Figure A3.

Left: Ages determined from BASTA, for stars that exhibit asteroseismic oscillations visible in Kepler data, and reported in Silva Aguirre et al. (2015). The y-axis shows ages determined without asteroseismic input to BASTA, using the inputs described in the text. The x-axis shows the ages determined when the same input parameters are used, in addition to asteroseismic parameters. A 1:1 line is shown in grey. Right: The x-axis shows asteroseismic ages determined by Silva Aguirre et al. (2015), for stars that overlap with our sample, plotted as a function of asteroseismic ages of the same stars in our analysis, on the y-axis.

The gold standard for determining ages of single stars is asteroseismic modelling of individual oscillation frequencies. Fig. A3 (right) shows our BASTA output ages (including the asteroseismic |$\Delta \nu$| and |$\nu _\text{max}$| as inputs) compared to ages determined from individual frequency modelling in Silva Aguirre et al. (2015). We find a median absolute offset of 1.0 Gyr and the standard deviation for the absolute residuals is 1.1 Gyr. For fractional residuals, we find an offset of |$6 \pm 2~{{\ \rm per\ cent}}$| and a residual scatter of 16  per cent, again showing good agreement with our estimated uncertainties.

To conduct a more extensive validation on a larger sample of stars, we apply our methods to a separate sample of stars. For this, we use the Kepler LEGACY sample, for which ages are precisely determined using asteroseismology (Lund et al. 2017; Silva Aguirre et al. 2017). Specifically, stellar parameters in this study were determined from frequency analysis of high quality, long baseline (minimum 12 months) photometric time series data from the Kepler mission. We use the stars in this sample which are within the stellar parameter range of our sample (5000 K < |$T_{\mbox{eff}}$| < 6500 |$\mathrm{K}$|⁠, |$0.7\,{\mathrm M}_\odot < M_{\star } < 1.4$|  |${\mathrm M}_\odot $|⁠). We extract the same Gaia parameters (⁠|$\varpi$|⁠, G, |$G_\mathrm{BP}$|⁠, |$G_\mathrm{RP}$|⁠), following the same data preparation steps, and perform fitting with BASTA and identical grids of stellar models for this sample as was done with the main sample in this paper. Fig. A4 shows the ages derived by our method relative to the asteroseismic ages from Silva Aguirre et al. (2017). We find a median absolute offset of 1.28 Gyr, a median residual scatter of 47  per cent, and the standard deviation for the absolute residuals is 2 Gyr. These results are consistent with our median uncertainties of 3.2 Gyr and 46  per cent.

Figure A4.

36 stars from the Kepler LEGACY sample with similar parameters to the stars in our sample. Ages determined with the method of this paper are plotted on the x-axis, against asteroseismic ages from Silva Aguirre et al. (2017) on the y-axis. The line of unity is shown by the upper line with smaller dashes. The lower line with larger dashes shows this line of unity shifted to the median offset of 1.28 Gyr. Four stars with a difference in metallicity of 0.15 dex or greater between Silva Aguirre et al. (2017) and Gaia are in faint, non filled points, with the remaining 36 in filled points.

A2.2 Gyrochronology and chemical abundances

As a further validation, we compare our ages to literature ages for these systems. Some were determined before the availability of Gaia DR3 parallaxes, and many use different stellar model tracks. Nevertheless, we compare the ages derived here to literature stellar ages, as shown in Fig. A5. Stars with ages calculated based on detailed chemical abundances, e.g. known low-iron hosts and stars with specific alpha element abundances reported, are highlighted in the blue squares, i.e. TOI-561, EPIC 249893012, K2-111, HD 20329, K2-291 (Kosiarek et al. 2019; Hidalgo et al. 2020; Bourrier et al. 2022; Murgas et al. 2022; Brinkman et al. 2023). Stars with ages derived from age–rotation relations and measured rotation periods are shown in green, i.e. WASP-47 and CoRoT-7-b (Almenara et al. 2016; John et al. 2022). Stars with asteroseismic ages are shown in red, i.e. Kepler-21, Kepler-10, Kepler-93, Kepler-100, from Silva Aguirre et al. (2015). All other stellar ages for this sample are shown in grey – these are primarily from isochrone ages from grids built with different models and differing input values on a star to star bases, making them difficult to compare. We see a good agreement for those stars where literature ages are available from asteroseismology, gyrochronology, or detailed stellar abundances. For all stars in this sample, including the grey points that are determined heterogeneously in the literature by isochrones, we find a median residual offset of |$10.8~{{\ \rm per\ cent}}$| and a scatter of |$29.1~{{\ \rm per\ cent}}$|⁠. For the three categories of stars with precise ages, we find a median scatter of 1 Gyr. After correcting for this systematic effect, we find a median scatter of 12.9  per cent, and the standard deviation for the absolute residuals is 0.97 Gyr. These results are consistent with our median uncertainties of 40.1  per cent and 3.2 Gyr.

Figure A5.

Ages from this work, which were used in the final analysis, plotted as a function of the most recent literature ages, for host stars in our initial sample. The 1:1 fit is shown with the dashed line. Ages determined in the literature from gyrochronology are shown with larger circles, those with well-constrained ages from measured |$[\rm{\alpha /Fe}]$| abundances are shown with squares, and those with asteroseismic literature values are shown with triangles.

Finally, we compare our ages to the star iron mass fraction as estimated by Adibekyan et al. (2021), for those stars common to both studies. The result is shown in Fig. A6. We observe that stars with enhanced value of |$F_{\text{iron}}^{\text{star}}$| are younger, whereas stars with enhanced alpha elements (and thus lower iron mass fraction) are older. This qualitatively confirms the expected picture from Galactic chemical evolution. The lowest iron-fraction star in the sample is TOI-561, which was independently confirmed as an old star that is a member of the chemical Galactic thick disc (Lacedelli et al. 2021; Weiss et al. 2021; Brinkman et al. 2023).

Figure A6.

Stellar iron mass fractions |$F_{\text{iron}}^{\text{star}}$| (Adibekyan et al. 2021) versus ages determined in this work. Stars with higher |$F_{\text{iron}}^{\text{star}}$|⁠, thus more iron compared to alpha elements, are younger. This is expected from Galactic chemical evolution, providing independent validation for our age determination.

A2.3 Kinematics

We further compare our ages to dynamical radial actions (⁠|$J_R$|⁠), which provide a measure for the radial excursions of a star’s orbit around the centre of the Milky Way. We calculate these using the `galpy'²  python package, using the MWPotential2014 potential (Bovy 2015). The radial actions are calculated according to

$$\begin{eqnarray} J_i = \frac{1}{2\pi } \oint _{\mathrm{orbit}} P_i \, \mathrm{d}x_i , \end{eqnarray}$$

(A1)

where i  = |$(r, \theta , z)$|⁠. Dynamical actions are a weak age indicator, due to the fact that orbits of stars born in the plane kinematically heat up over time due to interactions with substructure in the Galaxy (Beane, Ness & Bedell 2018). We plot age and normalized |$\sqrt{J_R}$| of host stars in the rocky planet sample, in Fig. A7. We observe a weak positive trend, in broad agreement with expected chrono-kinematic relations.

Figure A7.

Dynamical radial actions (⁠|$\sqrt{J_R}$|⁠) versus ages from this work. We broadly observe that stars with older ages in our sample correspond to stars that have undergone a higher degree of secular orbital evolution.

A2.4 Summary of age validation

The tests above demonstrate that our ages derived using Gaia input parameters have relative accuracies and precision of |$\approx$|2 Gyr (⁠|$\approx$|30  per cent for the median age in our sample). Since our sample is homogeneous, our conclusions are very minimally impacted by relative age accuracy, and are primarily based on relative age precision. The precision achieved with our method is consistent with independent studies who have used Gaia parallax to derive ages (Berger et al. 2018; Berger, Schlieder & Huber 2023a) and have identified demographic trends of exoplanet properties with age (e.g. Berger et al. 2020; Chen et al. 2021, 2022; Berger et al. 2023b; Yang et al. 2023).