At least one in a dozen stars exhibits evidence of planetary ingestion

The study seeks to detect and quantify chemical signatures of planets in stellar atmospheres, specifically abundance patterns correlated with condensation temperature (Tcond) that may arise from ingestion of planetary material or from planet formation processes that deplete refractories from the proto-stellar disc. Previous efforts have been hampered by small signal amplitudes, uncertain occurrence rates, and heterogeneous samples with varying ages. Because co-natal stars form from the same cloud with near-identical initial compositions, abundance differences between co-natal pairs can more cleanly reveal planet-related signatures. Using Gaia astrometry to identify co-moving, co-natal pairs and obtaining high-precision differential abundances, the authors test for ingestion signatures and separate them from other processes such as atomic diffusion by employing a Bayesian evidence framework applied to multiple competing models.

Prior work established that stellar surface abundances can be modified by accretion of planetary material and by planet formation depleting refractories (e.g., Pinsonneault et al. 2001; Hühn & Bitsch 2023; Meléndez et al. 2009; Booth & Owen 2020). Condensation temperature trends have been reported in solar twins and binaries (e.g., Adibekyan et al. 2014; Nissen 2015; Ramírez et al. 2015; Saffe et al. 2017; Oh et al. 2018; Nagar et al. 2020; Galarza et al. 2021). However, distinguishing planet signatures from Galactic chemical evolution, inhomogeneous ISM, and atomic diffusion is difficult, especially in heterogeneous samples. Gaia enables identification of co-moving, potentially co-natal pairs (Kamdar et al. 2019; Nelson et al. 2021). Atomic diffusion can drive subtle abundance changes (Dotter et al. 2017), complicating interpretation. Previous compilations suggested 20–35% of Sun-like stars may show ingestion signatures based on Fe anomalies (Spina et al. 2021), whereas Behmard et al. (2023) argued detections are rarer (<4.9%) when thermohaline mixing is efficient. There remains debate over the frequency and detectability of ingestion and the reliability of Tcond trends as indicators without rigorous model comparison.

Sample selection: Using Gaia EDR3 photometry and astrometry, 125 co-moving pairs (250 stars) were selected with cuts on 3D separation (Δs < 30 pc), 3D velocity separation (Δv < 2 km s−1), colour similarity (0.65 ≤ (B−PRP) ≤ 1.15 mag; |Δ(B−PRP)| ≤ 0.15 mag), absolute magnitude difference (ΔMG < 1 mag), brightness (G < 10), and by excluding larger groups (>4 members) via friends-of-friends to remove clusters and moving groups. Pairs were categorized as close/co-natal (Δs < 10^6 AU; 91 pairs) and far/control (Δs ≥ 10^6 AU; 34 pairs). Observations and data reduction: High-resolution, high S/N spectra were obtained with Magellan/MIKE (R≈50,000), Keck/HIRES (R≈72,000), and VLT/UVES (R≈110,000), achieving S/N≈250 per pixel at 600 nm. The sample primarily comprises late-F and G dwarfs. A strictly line-by-line differential approach was used between stars in each pair to minimize systematics. Stellar parameters: Effective temperature, surface gravity, microturbulence, and [Fe/H] were derived via a two-step differential method relative to the Sun and then pairwise, enforcing differential excitation/ionization balance and Fe I–Fe II equality, with a grid-search algorithm. Typical differential uncertainties: 15 K in Teff, 0.035 dex in log g, 0.012 dex in [Fe/H]. Abundance analysis: Equivalent widths were measured for 21 elements (C, O, Na, Mg, Al, Si, S, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Zn, Sr, Y, Ba, Ce), focusing on clean, intermediate-strength lines. Abundances were computed with MOOG (1D LTE) and ODFNEW atmospheres; hyperfine structure corrections were included for Sc, V, Mn, Co, Ba; 3D non-LTE corrections were applied to the O I 777 nm triplet; limited non-LTE checks for Na, Mg, Al, Mn indicated minor effects (~0.01–0.015 dex) in the differential analysis. Typical differential abundance uncertainties are ~0.015 dex across most species. Bayesian modeling: To interpret abundance differences Δ[X/H] within each pair, three models were fit using dynesty nested sampling to estimate posterior parameters and Bayesian evidence: (1) Planetary ingestion model, parameterized by the mass of ingested bulk-Earth composition material (ME) added to the surface convection zone of one star and an intrinsic scatter term σscatter. The convection zone mass fraction fcz was estimated from MESA-based grids as a function of Teff and stellar mass. Priors: ME ∈ [0.5,80] M⊕; σscatter ∈ [0,1]. (2) Flat model (null hypothesis) with an overall abundance offset Δ[M/H]offset and σscatter, testing for uniform shifts without Tcond structure. (3) Atomic diffusion model using MIST predictions for 6–8 elements (O, Na, Mg, Si, S, Fe, Ti; O and S not always available), fitting Δ[X/H]model,i = Δ[X/H]atom,i + Δ[M/H]offset with σscatter. For each pair, the Bayesian evidence ln(Z) was computed for each model; the key discriminants are Δln(Z) = ln(Z)planet − ln(Z)flat and Δln(Z)atom = ln(Z)planet − ln(Z)atom. Mock data tests: Synthetic mock “signal” (noise-free model-based) and mock “noise” (zero signal plus observational errors) samples of 125 pairs were generated to calibrate Δln(Z) thresholds and assess false positives. Tcond trend analysis: Linear weighted fits of Δ[X/H] versus Tcond were performed for each pair (various fitting ranges tested), and slopes with uncertainties were derived. However, due to degeneracies with inhomogeneous ISM and diffusion, Tcond slopes alone were not used as the primary indicator. Detection criteria: The study adopted three criteria for strong ingestion evidence: (i) Tcond slope > 4×10^5 dex K−1 with significance > 30; (ii) Δln(Z) > 3.5 (planet vs. flat); (iii) Δln(Z)atom > 2.5 (planet vs. diffusion).

• Sample and precision: 125 co-moving pairs analyzed (91 close/co-natal; 34 far/control). Achieved ~0.015 dex (3.5%) typical differential abundance uncertainties over 21 elements.
• Bayesian evidence distributions: Close/co-natal pairs show systematically higher ln(Z)planet and Δln(Z) than far/control pairs. Mock tests confirm Δln(Z) discriminates signal from noise.
• Detections: 11 co-natal pairs have Δln(Z) > 3.5; 7 co-natal pairs meet all three stringent criteria, constituting robust ingestion detections. Reported occurrence rate of ingestion among co-natal Sun-like stars: 8% ± 3% (Poisson uncertainty).
• Example case (Pair 124; HD 185726/HD 185689): Δln(Z) (planet−flat) = 12.3; Δln(Z)atom (planet−diffusion) = 5.7; Required accreted mass ME = 3.07 (+0.26, −0.22) M⊕ for HD 185689 to reproduce observed Δ[X/H]; strong positive Δ[X/H]–Tcond trend (slope 13.85 ± 0.75 ×10−5 dex K−1).
• Diffusion discrimination: In far/control pairs, over half have Δln(Z)atom < 0 and all have Δln(Z)atom < 2.5, consistent with diffusion contributing to observed differences; 10/11 strongest ingestion candidates have Δln(Z)atom > 2.5, indicating diffusion cannot explain their patterns.
• Tcond trends: Distributions of Tcond slopes are similar between close and far samples; Tcond trends alone cannot disentangle ingestion from other effects and may yield false positives.
• Accreted masses: For the seven candidates, average ingested mass is 4.3 ± 0.8 M⊕ (standard deviation 2.1 M⊕). Larger Δln(Z) generally correlates with larger inferred ME, modulated by fcz.
• Overall Tcond behavior: Across pairs, mean Tcond slope ≈ (3.5 ± 0.3) ×10−5 dex K−1; ingestion candidates have larger average slopes ≈ (6.9 ± 1.0) ×10−5 dex K−1.

The findings provide robust evidence that some Sun-like, co-natal stars have experienced ingestion of rocky, Earth-like material, leaving detectable chemical signatures in their photospheres. By using co-natal pairs and a Bayesian evidence framework to compare ingestion, flat, and atomic diffusion models, the study mitigates confounding effects from initial abundance inhomogeneities and diffusion. The 8% occurrence rate aligns with estimates from solar twins (15 ± 9%) and is below earlier compilations suggesting 20–35%, likely due to the present study’s homogeneous sample, higher precision, and multi-element modeling beyond Fe and C. The rate is marginally higher than recent predictions (<4.9%) that assume efficient thermohaline mixing, but thermohaline efficiency in 1D stellar models is poorly constrained and hydrodynamical studies suggest low efficiency, leaving room for detectable signatures. The magnitude of observed abundance anomalies (>0.05 dex up to 0.15 dex) favors late-time planet engulfment over subtler signatures from disc processes or terrestrial planet formation alone, which typically produce smaller differences (0.02–0.04 dex). N-body simulations predict most early instabilities (<100 Myr), but such early events may be diluted over Gyr; thus, the detections likely reflect later accretion triggered by outer perturbers (e.g., cold giants, stellar flybys) or long-term atmospheric erosion leading to instabilities. Given that 30–50% of Sun-like stars host inner super-Earths, the inferred ingestion rate suggests roughly one in 4–10 such systems might experience late ingestion, informing the long-term dynamical evolution of planetary systems. Alternative planet-related scenarios (disc accretion, refractory sequestration by terrestrial planets, or architectures with distant giants inhibiting accretion) may contribute in some systems, but the amplitudes observed here more strongly support engulfment.

This work presents high-precision, homogeneous differential abundances for 125 co-moving pairs, deploying a Bayesian model-comparison framework to distinguish chemical signatures of planetary ingestion from flat abundance offsets and atomic diffusion. At least seven robust ingestion cases are identified among 91 co-natal pairs, implying an occurrence rate of 8% ± 3% and approximately doubling the number of known ingestion instances. The results demonstrate that Tcond trends alone are insufficient and that Bayesian evidence-based indicators are effective for robust detection. The amplitudes of the detected signatures favor late-time engulfment events, offering new observational constraints on star–planet interactions and the long-term dynamical evolution of planetary systems. Future work should combine these abundance constraints with detailed dynamical modeling and improved treatments of stellar mixing (e.g., thermohaline and rotationally induced processes), extend non-LTE/3D analyses across more elements, and expand homogeneous co-natal samples to refine occurrence rates and timescales.

• The control (far) pairs may not be co-natal and thus can include intrinsic initial abundance differences, limiting their utility as a pure control for systematics.
• Atomic diffusion modeling used a limited set of elements (6–8) and relies on specific model grids; while relative evidence is robust, absolute model dependencies remain.
• Abundance analysis is based on 1D LTE (with select 3D/non-LTE corrections), which may leave residual systematics, though the differential approach mitigates many effects.
• Thermohaline and other mixing processes are uncertain; efficient mixing could reduce detectability of older ingestion events, potentially biasing occurrence rate toward more recent events.
• Sample completeness is ~45% at G < 10, and selection is brightness- and colour-limited; ingestion rates may vary with stellar type/age not fully captured here.
• Tcond trends are not uniquely diagnostic and can be mimicked by ISM inhomogeneities and diffusion; robust inference depends on the Bayesian model comparison framework.
• Early engulfment events (<100 Myr) may have signatures diluted over Gyr, so the reported rate likely reflects detectable, later events rather than total historical ingestion.