Isotopic Enrichment of Planetary Systems from Asymptotic Giant Branch Stars

²⁶Al and ⁶⁰Fe are short-lived radioactive isotopes (SLRs), with half-lives of 0.7 and 2.6 Myr, respectively (Castillo-Rogez et al. 2009; Wallner et al. 2015). Their decay isotopes are often found in chondritic meteorites, some of the oldest objects in our solar system, which suggests that these isotopes were present at the earliest epoch of planet formation around the Sun. Furthermore, ²⁶Al, and to a lesser extent, ⁶⁰Fe, are much more abundant than in the interstellar medium (ISM; Kita et al. 2013; Cook et al. 2021), indicating that the giant molecular cloud (GMC) that formed the Sun was either already enhanced in these SLRs (Young 2014), or some mechanism delivered them to the solar system as it was forming (Ouellette et al. 2010; Fatuzzo & Adams 2022).

In the scenario where the SLRs are inherited from the GMC, in order to obtain the required ²⁶Al/⁶⁰Fe ratio, the star-forming event that formed the Sun must have been sequential, with supernovae that produced ⁶⁰Fe triggering subsequent generations of stars which deliver the majority of the ²⁶Al later, via the winds of one or more Wolf–Rayet stars (Gounelle & Meynet 2012; Gounelle 2015). While apparently corroborated by circumstantial observational evidence that appears to show sequential triggered star formation (e.g., the Upper Sco complex; Preibisch & Zinnecker 1999), simulations of star formation do not predict such an efficient triggering process (Dale et al. 2015). Furthermore, it has been shown that a triggered star-forming region would dynamically merge into the region that triggered it (Parker & Dale 2016), resulting in age spreads (or even age dichotomies) exceeding 10 Myr, which are not observed (Jeffries et al. 2011).

The scenario in which SLRs are directly delivered to the solar system usually assumes the protosolar disk is enriched by the explosion of a nearby supernova (Ouellette et al. 2010; Lichtenberg et al. 2016). While this mechanism merely requires the Sun to form in a single populous star cluster, massive enough to contain stars that explode as supernovae (Nicholson & Parker 2017), the main issue is that any supernovae are unlikely to explode until 4 Myr (or much later, if stars >25 M_⊙ directly collapse to black holes rather than exploding as supernovae, e.g., Limongi & Chieffi 2018), by which time the protosolar disk will have evolved to form planets (or may have been destroyed by the ionizing radiation from the same massive stars that enrich the disk; Nicholson et al. 2019; Concha-Ramírez et al. 2021). This tension in timescales can be slightly mitigated if the low-mass disk-hosting star receiving the ejecta forms after the massive star(s), or is enriched by the winds of the massive stars (Portegies Zwart 2019; Parker et al. 2023). However, significant amounts of ⁶⁰Fe are only produced by the supernovae of massive stars, which do not explode until ∼10 Myr in the latest stellar evolution models (Limongi & Chieffi 2018).

An alternative production channel for ²⁶Al and ⁶⁰Fe is in the cores of asymptotic giant branch (AGB) stars (Karakas & Lugaro 2016). An AGB star is a post-main-sequence evolutionary phase undertaken by stars with initial masses 1–8 M_⊙ (Herwig 2005; Ventura et al. 2018). AGB stars dredge up their interiors and drive powerful winds, making the delivery of SLRs relatively straightforward (Trigo-Rodríguez et al. 2009). The issue with AGBs is that it was thought unlikely that an old, evolved star would have a chance encounter with the young Sun as it was forming planets (Kastner & Myers 1994).

However, the revolution in positional astronomy thanks to the Gaia mission has enabled researchers to accurately determine membership of star-forming regions, as well as being able to trace fast-moving "runaway" (RW; stars moving >30 km s⁻¹) and slower-moving "walkaway" (WW) stars (>5 km s⁻¹). In these analyses, interloping, or visiting stars can be disentangled from the host population of the star-forming region, demonstrating that older stars could encounter younger stars, and vice versa (Schoettler & Parker 2021).

In this Letter, we report the serendipitous discovery in Gaia DR3 of an interloping AGB star that has recently passed through a young star-forming region. We present the observational evidence in Section 2. We then model the enrichment of young stars and their protoplanetary disks using N-body simulations, and calculate the distribution of the yields of ²⁶Al and ⁶⁰Fe from an interloping AGB star in Section 3. We discuss caveats and conclude in Section 4.

While performing a search in Gaia Data Release 3 (DR3; Gaia Collaboration et al. 2021, 2023) for RW and WW stars in the vicinity of nearby star-forming regions, we discovered a WW star on the giant branch of the color–absolute-magnitude diagram (Figure 1) of NGC 2264, a young (∼3 Myr, e.g., Dahm 2008), relatively nearby (∼723 pc, Cantat-Gaudin & Anders 2020) star-forming region, postulated to have formed stars in a very dense configuration (Parker & Schoettler 2022).

Zoom In Zoom Out Reset image size

To identify interloping stars from Gaia DR3 data, we follow the method described in Schoettler et al. (2022) for Gaia DR2 data but update the information on the position and velocity of the center of NGC 2264 from Cantat-Gaudin & Anders (2020) and Carrera et al. (2019). These central values have lower associated uncertainties compared to those used in Schoettler et al. (2022), who used a different source for consistency with previous work.

We then collect the position and velocity information for all stars within 100 pc of this center from the Gaia archive. Instead of parallax, we use the (photogeometric) distance estimates from Bailer-Jones et al. (2021) for all stars. We apply a rest frame centered on the above values to all data and convert the position and velocity into a Cartesian reference frame. We then use a straight-line trace-back and identify any star as an interloper candidate if it can be traced back to a region within 2 pc in the x- and y-direction of the center (on the sky search radius as derived and used in Schoettler et al. 2022). We do not use the search radius definition in the z-direction (radial distance) from that analysis due to the large center distance errors. The NGC 2264 center distance from Cantat-Gaudin & Anders (2020) has a much smaller associated uncertainty, which translates to a search radius of 4 pc (2 pc as the on the sky radius + 2 pc distance uncertainty) in the z-direction.

We apply a maximum trace-back time of 5 Myr (as in Schoettler et al. 2022) and exclude any stars that were within our search region before that time. We then plot all the traced-back RW/WW candidates on a extinction/reddening-uncorrected color–absolute-magnitude diagram (Figure 1) and identify a star at a location inconsistent with the young age of NGC 2264.

The WW star (Gaia DR3 3131012157848982272) was traced back to the northern region of NGC 2264 using its proper motion and radial velocity. It has a velocity in the reference frame of NGC 2264 of 22.6 ± 1.8 km s⁻¹ and likely flew through this region ∼3 Myr ago. Gaia DR3 identifies this star as a long-period variable (LPV), a classification which encompasses AGB stars. We use the method to identify the subclass of this potential AGB star using Gaia DR3 and 2MASS photometry developed by Lebzelter et al. (2018). Our identified WW-LPV star is located firmly within the region of O-rich low-mass AGB stars on the (${W}_{\mathrm{RP},\mathrm{BP}-\mathrm{RP}}-{W}_{{K}_{s},J-{K}_{s}}$) versus K_s diagram (with the x-axis value of ∼0.2 mag and the y-axis value of ∼10.8 mag; after calculating its K_s magnitude at the LMC distance from its 2MASS K_s magnitude of ∼1.8 mag at ∼760 pc). This group (O-rich and low-mass) contains stars during the early-AGB and thermally pulsing AGB (TP-AGB) phases with initial masses from ∼0.9 to ∼1.4 M_⊙. Its position (${W}_{\mathrm{RP},\mathrm{BP}-\mathrm{RP}}-{W}_{{K}_{s},J-{K}_{s}}$: ∼0.2 mag and ${M}_{{K}_{s}}$: ∼−7.6 mag) coincides with the evolutionary track (as shown in Figure B1 in Lebzelter et al. 2018) of a star with an initial mass of ∼1.3 to ∼1.6 M_⊙.

It is unlikely that this old, evolved star formed in NGC 2264 (even the most conservative estimates place the duration of star formation at less than 10 Myr, e.g., Chevance et al. 2020), and its relatively fast velocity suggests that it has moved through NGC 2264, having formed elsewhere (we hereupon refer to this star as an "interloper").

Under the assumption that some star-forming regions may host interloping AGB stars, we use N-body simulations and published yields (Karakas 2014; Karakas & Lugaro 2016) to calculate the expected quantities of ²⁶Al and ⁶⁰Fe that could be accreted by a protoplanetary disk in a star-forming region with similar properties to NGC 2264. We assume the AGB star has deposited material from its winds at roughly the same time as star formation takes place, although we test this assumption by varying the size of the volume (and hence density) of AGB ejecta.

The simulations in question contain N_⋆ = 1000 stars drawn from a Maschberger (2013) IMF with a probability distribution of the form

$$\begin{eqnarray}&&p(m)\propto {\left(\displaystyle \frac{m}{\mu }\right)}^{-\alpha }{\left(1+{\left(\displaystyle \frac{m}{\mu }\right)}^{1-\alpha }\right)}^{-\beta }.\end{eqnarray} \tag{ 1 }$$

Here, μ = 0.2 M_⊙ is the average stellar mass, α = 2.3 is the Salpeter (1955) power-law exponent for higher-mass stars, and β = 1.4 describes the slope of the IMF for low-mass objects (which also deviates from the log-normal form; Bastian et al. 2010). We randomly sample this distribution in the mass range 0.1–50 M_⊙.

The resultant total mass for the star-forming regions are of order ∼500 M_⊙ (there is some variation due to the stochastic nature of sampling the stellar IMF), which lies toward the lower end of the observed mass function for star-forming regions (Lada & Lada 2003) and is similar to the mass of NGC 2264, the region in which we have observed an interloping AGB star. Such regions are more common than their higher-mass counterparts (e.g., Westerlund 1, R136), but rarer than low-mass star-forming regions (e.g., Taurus, ρ Oph).

The stars are all single (i.e., we do not include primordial binaries, and this simplification is unlikely to affect the numbers of stars that would encounter the interloping AGB star).

The simulated star-forming regions have initial radii of 1 pc, which results in a median stellar density of 1000 M_⊙ pc⁻³. Prestellar cores, and pre-main-sequence stars, are observed in spatially and kinematically substructured distributions (Cartwright & Whitworth 2004; Hacar et al. 2013; André et al. 2014) and for this reason we set up our simulations as self-similar fractal distributions (Goodwin & Whitworth 2004), where the stars' positions and velocities are correlated according to a fractal dimension D = 2.0. Similarly, the star-forming regions are initially set up with subvirial velocities where the virial ratio α_vir = T/∣Ω∣, and T and ∣Ω∣ are the total kinetic energy and total potential energy of the stars, respectively. In this definition, α_vir = 0.5 means the system is in virial equilibrium. We adopt α_vir = 0.3 to mimic the initial velocity distributions in observed (e.g., Foster et al. 2015) and simulated (e.g., Bate 2012) star-forming regions.

We assume AGB stars that are identified as RW and WW stars will have a similar velocity distribution to stars ejected from star-forming regions, and so we select the median velocity from the distribution in Schoettler et al. (2019). We assume the AGB star traverses the star-forming region almost immediately after star formation (such that the protoplanetary disks are yet to form planets). We note that this occurred for the AGB star that passed through NGC 2264; it interloped through the region 3 Myr ago, which is roughly the age of the stars in the region.

The ejecta from the wind of the AGB star is modeled as a cylinder with length 4 times the half-mass radius of the star-forming region, l_ejecta = 4r_1/2 (note that the half-mass radius expands as the star-forming region evolves, meaning that the cylinder of AGB ejecta also increases in length, reducing the density of the ejecta). We experimented with varying the radii of the cylinder of ejecta, r_ejecta; smaller radii (e.g., r_ejecta = 0.1 pc) lead to high-level enrichment of a few stars, whereas larger radii (0.5–1 pc) lead to enrichment of more stars, but at lower values.

A cartoon of the AGB ejecta geometry is shown in Figure 2. The AGB star (green) has traversed the star-forming region, depositing a cylindrical shaped ejecta. Some stars (shown in red) are enriched and then leave the ejecta region, whereas others remain. Blue stars are those that do not capture any of the ejecta (including fore- and background stars that reside outside of the ejecta). The hypothetical young Sun is shown by the yellow star.

Zoom In Zoom Out Reset image size

We determine whether a star crosses through the AGB ejecta, and if it does we calculate the distance it travels, d_trav. We use this, and the radius of the disk around the star, r_disk, to determine the amount of AGB ejecta swept up from the cylinder:

$$\begin{eqnarray}&&{\eta }_{\mathrm{ejecta}}=\displaystyle \frac{{r}_{\mathrm{disc}}^{2}{d}_{\mathrm{trav}}}{{l}_{\mathrm{ejecta}}{r}_{\mathrm{ejecta}}^{2}}.\end{eqnarray} \tag{ 2 }$$

This equation assumes that none of the AGB ejecta is deflected by the disk (i.e., the relative velocities are low, unlike when a disk encounters a supernova blast wave, Ouellette et al. 2007). Furthermore, this represents the maximum possible amount of material swept up, as we do not account for the inclination angle of the disk (although the average angle is likely to be in the region of 60°, Lichtenberg et al. 2016, which would only reduce the enrichment by a factor of 2).

We then determine the mass of ²⁶Al, ${m}_{{}^{26}\mathrm{Al}}$, swept up by the protostellar disk by dividing the total ²⁶Al yield from the AGB star, ${m}_{{}^{26}\mathrm{Al},\ \mathrm{AGB}}$, by the time spent in the AGB phase t_AGB and then multiplying this by the time taken for the star to travel d_trav, which is Δt:

$$\begin{eqnarray}&&{m}_{{}^{26}\mathrm{Al}}={\eta }_{\mathrm{ejecta}}\displaystyle \frac{{m}_{{}^{26}\mathrm{Al},\ \mathrm{AGB}}}{{t}_{\mathrm{AGB}}}{\rm{\Delta }}t.\end{eqnarray} \tag{ 3 }$$

We perform a similar calculation to determine the mass of ⁶⁰Fe swept up by the disk, ${m}_{{}^{60}\mathrm{Fe}}$:

$$\begin{eqnarray}&&{m}_{{}^{60}\mathrm{Fe}}={\eta }_{\mathrm{ejecta}}\displaystyle \frac{{m}_{{}^{60}\mathrm{Fe},\ \mathrm{AGB}}}{{t}_{\mathrm{AGB}}}{\rm{\Delta }}t,\end{eqnarray} \tag{ 4 }$$

where ${m}_{{}^{60}\mathrm{Fe},\ \mathrm{AGB}}$ is the amount of ⁶⁰Fe produced by the AGB star, and the time variables are as above.

In order to calculated the amount of ²⁶Al and ⁶⁰Fe relative to the stable isotopes in the protoplanetary disks (²⁷Al and ⁵⁶Fe), we make some assumptions about the masses and sizes of the protoplanetary disks.

For each planet-hosting star (masses between 0.1 < M_⋆/M_⊙ < 3 M_⊙), we assign it a disk of mass

$$\begin{eqnarray}&&{M}_{\mathrm{disc}}=0.1\,{M}_{\star },\end{eqnarray} \tag{ 5 }$$

and a radius r_disk = 400 au, commensurate with the observed disks in nearby star-forming regions (e.g., Andrews et al. 2013; Ansdell et al. 2016; Barenfeld et al. 2017; Tazzari et al. 2017; Eisner et al. 2018; Cieza et al. 2019). We assume that the disks do not lose mass, nor do their radii evolve inwards due to external photoevaporation from massive stars, or outwards due to viscous spreading. We assume a gas-to-dust ratio of 100:1 such that the mass of solids in the disk is

$$\begin{eqnarray}&&{m}_{\mathrm{dust}}=0.01\,{M}_{\mathrm{disc}}.\end{eqnarray} \tag{ 6 }$$

The amount of dust that is ²⁷Al is given by

$$\begin{eqnarray}&&{m}_{{}^{27}\mathrm{Al}}=8500\times {10}^{-6}{m}_{\mathrm{dust}},\end{eqnarray} \tag{ 7 }$$

and the amount of dust that is ⁵⁶Fe is given by

$$\begin{eqnarray}&&{m}_{{}^{56}\mathrm{Fe}}=1828\times {10}^{-4}{m}_{\mathrm{dust}},\end{eqnarray} \tag{ 8 }$$

(Lodders 2003). We then use the mass of SLR swept up by the disk to calculate the yields of ²⁶Al and ⁶⁰Fe thus:

$$\begin{eqnarray}&&{Z}_{\mathrm{Al}}=\displaystyle \frac{{m}_{{}^{26}\mathrm{Al}}}{{m}_{{}^{27}\mathrm{Al}}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{Z}_{\mathrm{Fe}}=\displaystyle \frac{{m}_{{}^{60}\mathrm{Fe}}}{{m}_{{}^{56}\mathrm{Fe}}}.\end{eqnarray} \tag{ 10 }$$

In Figure 3 we show the ²⁶Al/²⁷Al ratio as a function of the ⁶⁰Fe/⁵⁶Fe ratio, with the initial ²⁶Al/²⁷Al ratio measured in the solar system shown by the horizontal dashed line (Thrane et al. 2006). The measurement of the abundance of ⁶⁰Fe in the early solar system is more controversial, with estimates varying by several orders of magnitude, from ⁶⁰Fe/⁵⁶Fe ∼10⁻⁸ − 10⁻⁶ (Tang & Dauphas 2012; Mishra et al. 2016; Cook et al. 2021; Kodolányi et al. 2022). The maximum range of values for the initial solar system ⁶⁰Fe/⁵⁶Fe ratio is shown by the vertical dotted–dashed line (Mishra et al. 2016), and vertical dotted line (Tang & Dauphas 2012).

Zoom In Zoom Out Reset image size

We plot the ²⁶Al/²⁷Al ratio and the ⁶⁰Fe/⁵⁶Fe ratio for different AGB progenitor masses; the black points are for 7 M_⊙, the yellow points for 6 M_⊙, the blue points for 5 M_⊙ and the green points are for 3 M_⊙ progenitor masses. All other parameters (density of AGB ejecta, radius of accreting protoplanetary disk, etc.) are kept constant. The lighter colored points indicate stars of roughly solar mass (0.5–1.5 M_⊙). Figure 3 indicates that a reasonably high initial AGB progenitor mass (6–7 M_⊙) is required to provide solar system–like abundances, although these values can also change if we increase the stellar density in the star-forming region, or the radius of the accreting protoplanetary disk. Note that the AGB star we have found interloping through NGC 2264 has a lower progenitor mass (∼1.5 M_⊙), and such stars are more common than the 6–7 M_⊙ star needed to enrich the solar system. However, we merely wish to demonstrate here that AGB stars can interlope through star-forming regions, and a higher-mass AGB star can provide the enrichment found in the solar system.

At this stage, we do not account for the radioactive decay of the SLRs. This is because we assume the AGB deposits material at the instant of star formation in our simulations, and due to the relatively high stellar density, most (>90%) of the AGB ejecta is swept up in the first 0.5–1 Myr.

The presence of short-lived radioisotopes provides an additional heat source in the interiors of planets, especially if the SLRs are incorporated early in the formation of the planetary system, i.e., before differentiation within the individual bodies has occurred. The contribution from ²⁶Al and ⁶⁰Fe dominated the radiogenic heat budget of the early Earth (McDonough 2020), with the amount of heating calculated at Q(t) = 3.5 × 10⁻⁷ W kg⁻¹, assuming ²⁶Al/²⁷Al = 5.85 × 10⁻⁵ (Thrane et al. 2006) and ⁶⁰Fe/⁵⁶Fe = 1 × 10⁻⁶ (Mishra et al. 2016).

For each star in our simulations, we use the abundance of ²⁶Al (defined by the ratio ²⁶Al/²⁷Al) and the abundance of ⁶⁰Fe (defined by the ratio ⁶⁰Fe/⁵⁶Fe) to calculate the geophysical heating (Moskovitz & Gaidos 2011) at a given time after the AGB star has deposited the SLRs in the star-forming region (here, we do account for the decay of the SLRs) using

$$\begin{eqnarray}&&Q(t)={f}_{\mathrm{Al},\mathrm{CI}}{Z}_{\mathrm{Al}}\displaystyle \frac{{E}_{\mathrm{Al}}}{{\tau }_{\mathrm{Al}}}{e}^{-t/{\tau }_{\mathrm{Al}}}+{f}_{\mathrm{Fe},\mathrm{CI}}{Z}_{\mathrm{Fe}}\displaystyle \frac{{E}_{\mathrm{Fe}}}{{\tau }_{\mathrm{Fe}}}{e}^{-t/{\tau }_{\mathrm{Fe}}},\end{eqnarray} \tag{ 11 }$$

where f_Al,CI is the fraction of Al in chondrites (Lodders 2003), E_Al = 3.12 MeV is the decay energy of ²⁶Al, and τ_Al = 0.717 Myr is the radioactive half-life of ²⁶Al (Castillo-Rogez et al. 2009). Similarly, f_Fe,CI is the fraction of Fe in chondrites (Lodders 2003), E_Fe = 2.712 MeV is the decay energy of ⁶⁰Fe (Castillo-Rogez et al. 2009) and τ_Fe = 2.6 Myr is the half-life of ⁶⁰Fe (Wallner et al. 2015). The initial solar system heating value is calculated from these values to be Q_SS = 3.4 × 10⁻⁷ W kg⁻¹.

We calculate the geophysical heating at various times after star formation. We assume that this is the same length of time after the production of AGB ejecta, and the resultant distributions are shown in Figure 4. The different histograms correspond to the heating after 1 (blue), 5 (red), and 10 Myr (gray), respectively. The vertical dotted line indicates the heating value, Q(t), calculated for the early solar system (Moskovitz & Gaidos 2011). The early heating of protoplanetary disks enriched by AGB ejecta is consistent with the value calculated for the solar system, and we find that several combinations of our initial conditions (simulation density, AGB progenitor mass, protoplanetary disk radius, density of AGB ejecta) are within this range.

Zoom In Zoom Out Reset image size

The main caveat to our results is that the proportion of ⁶⁰Fe in comparison to ²⁶Al is toward the high end of the range of values for the measured ratio in the solar system. This could be alleviated if the enriching AGB star has a lower progenitor mass (in some cases the ratio of ²⁶Al/⁶⁰Fe is then higher, e.g., the green points in Figure 3, which are the yields from a 3 M_⊙ progenitor star), and there is an additional source of ²⁶Al, e.g., from the wind(s) of massive stars in the star-forming region. However, the latter scenario would then require the Sun to form in a star-forming region also containing massive stars, which somewhat limits the advantages of enrichment from AGB stars. This would also return us to the arguments against the direct enrichment of the solar system by massive stars (only a small fraction of star-forming regions contain massive stars, and those massive stars FUV and EUV radiation fields could evaporate planet-forming disks).

Our results depend slightly on the volume density of the AGB ejecta; the lower the density of the ejecta, the lower the enrichment. In our default simulations, the ejecta is dispersed in a cylinder of 0.1 pc, but solar system levels of enrichment can occur when the cylinder has a larger radius. Furthermore, the most likely AGB progenitor mass is between 5 and 8 M_⊙, but we note that a lower progenitor mass is possible if the star-forming region is initially very dense (as is the case for NGC 2264; Schoettler et al. 2022). We do not know the initial density of the star-forming region in which the solar system formed, but based on the frequency of dynamical interactions that would truncate the Sun's protosolar disk and disrupt the early solar system, the initial stellar density could be slightly higher than those we adopt in our simulations (>10⁴ M_⊙ pc⁻³; Pfalzner & Vincke 2020).

We also ran test simulations where we lowered the stellar density to 100 M_⊙ pc⁻³. In this scenario, a similar amount of material is swept up, but the enrichment occurs later, after around 5 Myr, once the star-forming region has attained its maximum density. We would therefore expect some of the ²⁶Al and ⁶⁰Fe to have already decayed, unless the AGB star deposited the material several Myr after star formation.

If we instead reduce the number of stars (but keep the stellar density constant at 1000 M_⊙ pc⁻³), then for a single star-forming region fewer stars are enriched overall, but the fraction of stars that are enriched is very similar to the simulations with higher N_⋆ (because the stellar densities are similar).

The size of the protoplanetary disk in our simulations is fixed at 400 au, commensurate with the disk sizes measured in nearby star-forming regions. As expected, reducing the initial disk radii reduces the cross section for enrichment, but we note that viscous expansion of smaller disks means that our adopted disk radii could be conservative estimates after several Myr of evolution (Concha-Ramírez et al. 2019).

In summary, we have shown that evolved stars can encounter forming planetary systems, and the yield of ²⁶Al and ⁶⁰Fe from AGB stars can account for the enrichment and subsequent geophysical heating in the solar system.

We thank the anonymous referee for a prompt and helpful report. R.J.P. acknowledges support from the Royal Society in the form of a Dorothy Hodgkin Fellowship. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation program (grant agreement No. 853022, PEVAP).