C3PO – IV. Co-natal stars depleted in refractories are magnetically more active – possible imprints of planets

ABSTRACT

Chemical abundance anomalies in twin stars have recently been considered tell-tale signs of interactions between stars and planets. While such signals are prevalent, their nature remains a subject of debate. On the one hand, exoplanet formation may induce chemical depletion in host stars by locking up refractory elements. On the other hand, exoplanet engulfment can result in chemical enrichment, and both processes potentially produce similar differential signals. In this study, we aim to observationally disentangle these processes by using the Ca ii infrared triplet to measure the magnetic activity of 125 co-moving star pairs with high signal-to-noise ratio, and high-resolution spectra from the Magellan, Keck, and VLT (Very Large Telescope) telescopes. We find that co-natal star pairs in which the two stars exhibit significant chemical abundance differences also show differences in their magnetic activity, with stars depleted in refractories being magnetically more active. Furthermore, the strength of this correlation between differential chemical abundances and differential magnetic activity increases with condensation temperature. One possible explanation is that the chemical anomaly signature may be linked to planet formation, wherein refractory elements are locked into planets, and the host stars become more active due to more efficient contraction during the pre-main-sequence phase or star–planet tidal and magnetic interactions.

1 INTRODUCTION

It has long been observed that compared to |${\sim}85$| per cent of solar twins, the Sun exhibits a depletion in the elemental abundance of refractory elements (such as Fe, Ti, and Al),¹ and this depletion is correlated with the condensation temperature (e.g. Meléndez et al. 2009; Bedell et al. 2018; Rampalli et al. 2024). Similar phenomena have also been observed and validated in other stars (e.g. Ramírez et al. 2010; Ramírez, Meléndez & Asplund 2014; Liu et al. 2020, 2024; Flores et al. 2024).

The connection between this phenomenon and planets has been investigated by the following studies. For example, Meléndez et al. (2009) and Ramírez, Meléndez & Asplund (2009) proposed that the sequestration of these refractory elements by terrestrial planets may explain this depletion. Additionally, Booth & Owen (2020) suggested that giant planets, such as Jupiter, could contribute by creating dust traps, thereby preventing the fall of refractory-rich dust on the host star. Recently, Hühn & Bitsch (2023) conducted numerical simulations on the abundance differences resulting from such dust traps and found that the measurements of the elemental abundance differences for a wide binary system, which hosts a giant planet, align reasonably well with their simulations.

Unlike chemical depletion that may be caused by exoplanet formation processes, chemical enrichment can result from the ingestion of planetary material (e.g. Chambers 2010; Hühn & Bitsch 2023). Based on this hypothesis, Liu et al. (2024) reported evidence of planetary ingestion in at least 8 per cent of stars, a rate comparable to observations in solar twins (Ramírez et al. 2009) and N-body simulations (Bitsch & Izidoro 2023). This conclusion is supported by the consistency between observed elemental abundance differences and theoretical models of planetary engulfment.

In this work, we revisit the data set as in Liu et al. (2024) but focus on magnetic activity, which is directly derived from the spectra. Magnetic activity is closely linked to stellar rotation, with stars exhibiting increasingly stronger activity as rotation rates increase. This activity reaches saturation at rapid rotation, characterized by rotation periods of less than |${\sim}3$| d (see Işık et al. 2023, for a recent review). In close binaries, due to tidal interactions, stars can spin more rapidly when their rotation becomes synchronized with orbital periods. Consequently, these binary stars, even at the old ages, may exhibit enhanced activity, as seen in red giants (e.g. Gehan, Gaulme & Yu 2022). Therefore, it is likely that hot Jupiters and their hosts may undergo analogous interactions.

Tidal interactions between planets and their host stars have been extensively studied in the literature (see recent reviews by Ogilvie 2020; Lanza 2022), and these interactions manifest primarily in two mechanisms: equilibrium tides and dynamical tides (e.g. Ogilvie & Lin 2004; Mathis 2015). Equilibrium tides result from the lag between the tidal bulge raised on the star by the planet and the planet’s position in its orbit. This lag transfers angular momentum between the orbit and the star’s spin, altering the star’s rotation over time (Ogilvie 2014). Complementing equilibrium tides are dynamical tides, which excite stellar oscillations through the gravitational potential of the planet and the Coriolis force of the star (Ogilvie & Lin 2007). The efficiency of dynamical tides depends on stellar properties, particularly the thickness of the outer convective zone (Mathis 2015).

In accordance with theoretical predictions on tides, observations by Pont (2009) and Brown et al. (2011) revealed that certain host stars of hot Jupiters exhibit unexpectedly rapid rotation rates. Additionally, Poppenhaeger & Wolk (2014) observed that host stars expected to experience strong tidal interactions with their planets display enhanced X-ray activity compared to their counterparts in wide binary systems devoid of known planets. Furthermore, direct instances of magnetic star–planet interactions have been found, where stars hosting short-period planets exhibit significantly higher flare rates (e.g. Kavanagh et al. 2021; Feinstein et al. 2024), or excess flaring aligned with the planet’s phase or close to periastron (e.g. Shkolnik et al. 2005; Pillitteri et al. 2011; Maggio et al. 2015; Ilin et al. 2024).

In addition to tidal and magnetic interactions, giant planets can further shape the rotational and magnetic evolution of their host stars. As Booth & Owen (2020) suggested, the formation of dust gaps can magnetically decouple the host star from the planet, enabling more efficient stellar contraction during the pre-main-sequence (PMS) phase. This process shortens the lifetime of the disc and causes the host star to spin-up more rapidly. As a result, these stars exhibit faster rotation rates and stronger magnetic activity during the main sequence (e.g. Monsch et al. 2023). Thus, the presence of giant planets may influence not only stellar chemistry but also the rotation and magnetic activity of their host stars. However, few studies have investigated the connection between magnetic activity and the observed chemical anomalies in these stars. The motivation for this work is to observationally explore whether magnetic activity can help disentangle the closely related differential signals seen in stellar chemical anomalies.

2 DATA AND METHODS

2.1 Data and sample selection

Understanding the relationship between stellar chemical anomalies and planetary processes requires careful analysis of stars that share similar initial conditions. Co-moving and spatially close stellar pairs, often classified as co-natal, provide a valuable opportunity for such studies (Kamdar et al. 2019; Nelson et al. 2021). Co-natal stars are thought to have originated from the same molecular cloud, implying a shared chemical history and evolutionary pathway. By comparing differential elemental abundances in these stars, we can isolate variations that may arise from external factors such as interstellar medium (ISM) inhomogeneity and atomic diffusion. This investigation forms the basis of the Complete Census of Co-moving Pairs Of stars (C3PO) program, which is dedicated to high-precision, homogeneous differential elemental abundance analysis of co-moving star pairs (Yong et al. 2023; Liu et al. 2024; Sun et al. 2025a).

The C3PO program collected high-resolution, high signal-to-noise ratio (⁠|${\rm S/N} \approx 250$| pixel⁻¹ at 600 nm) spectra using three major instruments: the Magellan Inamori Kyocera Echelle (MIKE) spectrograph on the Magellan Telescope over 7 nights (⁠|$R \approx 50\, 000$|⁠), High Resolution Echelle Spectrometer (HIRES) on the Keck Telescope over 1 night (⁠|$R \approx 72\, 000$|⁠), and Ultraviolet and Visual Echelle Spectrograph (UVES) on ESO’s Very Large Telescope (VLT) over 26.4 h (⁠|$R \approx 110\, 000$|⁠). Each pair of stars in our sample has similar colour (⁠|$|\Delta ({\rm BP} - {\rm RP})| \le 0.15$|⁠) and absolute magnitude (⁠|$|\Delta M_G| < 1$|⁠) in the Gaia G band. This similarity between pairs facilitated the line-by-line differential analysis, a crucial approach for achieving high precision in stellar parameters and relative elemental abundances (e.g. Meléndez et al. 2009; Sun et al. 2020). Additionally, this method allows us to cancel out major sources of systematic uncertainties, such as those arising from model atmospheres and atomic line data (Bedell et al. 2014; Nissen & Gustafsson 2018).

For this study, we revisit 125 co-moving star pairs selected from the C3PO program (Yong et al. 2023; Liu et al. 2024; Sun et al. 2025a). Each pair consists of a reference star and an object star, both of which exhibit similar fundamental properties. The majority of the pairs show differences of |$\Delta \rm{T_{\rm{eff}}} < 300 \, {\rm K}$|⁠, |$\Delta \rm{\log g} < 0.3$|⁠, and |$\Delta [{\rm Fe/H}] < 0.3$|⁠. The distinction between reference and object stars is arbitrary, meaning the reference star can have either higher or lower elemental abundances than the object star. Fig. 1 illustrates the distribution of these stars on the Kiel diagram, with each pair connected by a purple line, demonstrating their proximity in stellar parameters. Most of the stars lie along the main sequence, indicating that this sample represents a relatively young stellar population. The median ages, determined through isochrone fitting, are 3.25 Gyr for the object stars and 3.50 Gyr for the reference stars. The median absolute age difference between the object and reference stars is 1.30 Gyr, which is comparable to the median age uncertainty (1.68 Gyr), suggesting that their ages are equal to within uncertainties and that the paired stars are consistent with being co-natal. These ages were derived using the isoclassify pipeline (Huber et al. 2017; Berger et al. 2020), based on |$T_{\text{eff}}$|⁠, |$\log g$|⁠, and [Fe/H] measurements, along with their respective uncertainties.

Following the methodology of Liu et al. (2024), we classify 91 of the 125 co-moving systems as co-natal, based on 3D separations |$\Delta s< 10^6$| au and 3D velocity differences |$\Delta v < 2.0$| km s⁻¹. Among these co-natal systems, 65 exhibit homogeneous chemistry, defined as |$|\rm{\Delta [\rm {Fe/H]}}| < 3\sigma _{\rm{\Delta [\rm {Fe/H]}}}$|⁠. The remaining 26 systems show chemical anomalies, characterized by |$|\rm{\Delta [\textrm {Fe/H]}}| \ge 3\sigma _{\rm{\Delta [\rm {Fe/H]}}}$|⁠. Here, |$\Delta$|[Fe/H] represents the metallicity difference between the two stars in each pair, and |$\sigma _{\rm{\Delta [\textrm {Fe/H]}}}$| corresponds to the uncertainty in |$\Delta$|[Fe/H].

We adopt stellar parameters and elemental abundances from Liu et al. (2024), which were derived using the high-precision spectra described above. The typical uncertainties in the differential stellar parameters are around 15 K for |$T_{\text{eff}}$|⁠, 0.035 for |$\log g$|⁠, and 0.012 dex for [Fe/H]. Based on these differential stellar parameters, elemental abundances for 21 elements (C, O, Na, Mg, Al, Si, S, Ca, Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Zn, Sr, Y, Ba, and Ce) were calculated on a line-by-line basis. The average errors in differential abundances for nearly all co-moving pairs are approximately 0.015 dex (3.5 per cent).

In Fig. 2, we present four co-moving systems that exhibit strong correlations between elemental abundance difference (⁠|$\Delta$|[X/H]) and condensation temperature. The slopes of these correlations vary between positive and negative due to the arbitrary designation of reference and object stars in each pair. These correlations align with previous observations for several binaries of stellar twins hosting planets or debris discs (see Melendez & Ramirez 2016, for a review). In the next section, we introduce a differential activity index (i.e. the activity difference between the two stars in each pair) as an additional diagnostic to further investigate the chemical anomalies of these co-natal stars.

2.2 Activity index

We employ the residual equivalent width of the Ca ii infrared triplet lines as a chromospheric activity indicator in this study, denoted as |$\Delta W_{\textrm {IRT}}$|⁠. This metric is defined as the sum of the flux differences between the flux-normalized spectra of the paired stars, calculated as the flux of the object star minus that of the reference star, consistent with the ordering used in computing |$\Delta$|[X/H]. The summation is evaluated over a narrow window, 0.75 Å wide, centred around each of the Ca ii infrared triplet lines (8498, 8542, and 8662 Å). This |$\Delta W_{\textrm {IRT}}$| metric has been validated in the literature, where it has demonstrated a strong correlation with the chromospheric activity diagnostic |$\log R^{\prime }_{\text{HK}}$|⁠, derived from the Ca ii H&K lines (Busà et al. 2007; Martin et al. 2017). It has also been widely used to study stellar activity in large-scale surveys such as RAVE (Radial Velocity Experiment; Žerjal et al. 2013), Gaia Radial Velocity Spectrometer (RVS; Lanzafame et al. 2023), and LAMOST (Large Sky Area Multi-Object Fiber Spectroscopic Telescope; Huang et al. 2024). Figs 3 and 4 show examples of the spectra and the integration windows for two-star pairs – one with minimal difference and one with a pronounced difference in their line cores.

There are 82 systems with spectra covering all three triplet lines, 19 systems with spectra covering only the first two lines, and 24 systems that lack spectral coverage of the triplet lines entirely. To account for this variation, we sum the flux differences within the windows and then calculate the average by dividing the total flux difference by the number of lines measured for each star. Importantly, we verified that the number of lines used for calculating |$\Delta W_{\textrm {IRT}}$| does not affect our final conclusions, as confirmed by tests using a single, randomly chosen Ca ii infrared triplet line. The window width used in our analysis, 0.75 Å, is half of that employed by Lanzafame et al. (2023), chosen to minimize the influence of line wings, which are unrelated to stellar activity. However, testing with a wider window of 1.5 Å showed no significant impact on our results. All spectra are shifted to the rest frame using radial velocities derived from this study by fitting a composite model, which consists of a line and a Gaussian for each absorption feature. We find this method to be particularly effective for precisely identifying local features, such as the Ca ii triplet lines, and it differs from the commonly used cross-correlation function methods, which are typically applied over broader wavelength ranges. Initial radial velocity estimates were taken from Liu et al. (2024) for the fitting process.

To estimate the uncertainty in |$\Delta W_{\textrm {IRT}}$|⁠, we propagate the flux errors in the two spectra for the reference and object stars. Given two spectra with flux values |$f_1(\lambda)$| and |$f_2(\lambda)$| at each wavelength |$\lambda$|⁠, the flux difference is defined as

The corresponding uncertainty in |$\Delta f(\lambda)$| follows standard error propagation:

where |$\sigma _{f_1}(\lambda)$| and |$\sigma _{f_2}(\lambda)$| are the respective flux uncertainties of the two spectra. Given that flux uncertainties are not available for the C3PO spectra (Yong et al. 2023), we estimate the S/N by selecting a continuum region near each line of the Ca ii infrared triplet (8474–8490, 8518–8532, and 8626–8646 Å). Within each of these three regions, we compute the median flux |$F_{\mathrm{cont}}$| and the |$1\sigma$| scatter of the flux from the 16th and 84th percentiles of the flux distribution. The S/N is then estimated by |$F_{\mathrm{cont}} / \sigma$|⁠. With S/N, the flux uncertainties can be expressed as

where |$\mathrm{(S/N)}_1$| and |$\mathrm{(S/N)}_2$| are the signal-to-noise ratios of the two spectra. Since |$\Delta W_{\textrm {IRT}}$| is defined as the sum of flux differences over a given wavelength range:

the total uncertainty in |$\Delta W$| is given by

Substituting |$\sigma _{\Delta f}(\lambda)$|⁠, we obtain

Similar to the calculation of |$\Delta W_{\textrm {IRT}}$|⁠, when multiple lines in the triplet are observationally available, we compute the final uncertainty |$\sigma _{\rm{\Delta W_{\rm {IRT}}}}$| by averaging over all N lines, namely,

To validate our |$\Delta W_{\textrm {IRT}}$| measurements, we calculate another commonly used chromospheric activity index, |$S_{\textrm {HK}}$|⁠, based on the Ca ii H&K lines (3968 and 3934 Å, respectively), following standard methodology (see e.g. Zhang et al. 2020; Gehan et al. 2022; Cordoni et al. 2024; Yu et al. 2024). It is important to note that |$\Delta S_{\textrm {HK}}$| is used solely for validation purposes and |$\Delta W_{\textrm {IRT}}$| serves as the primary metric for our subsequent analysis. This choice is due to wavelength coverage limitations in our spectra as determined by the observation settings – Magellan/MIKE (3350–9300 Å), Keck/HIRES (4200–8500 Å), and VLT/UVES (4800–6800 Å, whose spectra were thus not used for measuring activity indexes; see Yong et al. 2023). With this spectral coverage, we measure |$\rm{\Delta S_{\rm {HK}}}$| for 40 per cent of the systems and |$\Delta W_{\textrm {IRT}}$| for 81 per cent of the star pairs.

To calculate |$\Delta S_{\textrm {HK}}$|⁠, we first identify the Ca ii H&K lines by redetermining the radial velocities using the |${\rm H} \, {\alpha }$| line (⁠|$\lambda _0 = 6563$| Å in air), applying the same method used for identifying the Ca ii triplet lines. We calculate the integrated emission line fluxes within the H and K bandpasses, using 1.09 Å full width at half-maximum (FWHM) triangular windows. The continuum fluxes within the R and V bandpasses (centred at 4001 and 3901 Å, respectively) are summed using 20 Å rectangular windows. We then calculate |$S_{\textrm {HK}}$| as follows (Karoff et al. 2016):

where the S-index calibration factor is |$\alpha = 1.8$| (Hall, Lockwood & Skiff 2007). The |$S_{\textrm {HK}}$| measurements are then used to calculate the differential activity index, |$\Delta S_{\textrm {HK}}$|⁠, which is the difference in |$S_{\textrm {HK}}$| between the object and reference stars.

Fig. 5 demonstrates a clear correlation between our two differential activity indices (⁠|$\Delta W_{\textrm {IRT}}$| and |$\Delta S_{\textrm {HK}}$|⁠) across the full temperature range of our sample. We note that the typical uncertainties in |$\Delta W_{\textrm {IRT}}$| (median = 0.09) are smaller than the symbol size – primarily due to the high S/N of the C3PO spectra – indicating that additional factors contribute to the scatter in the correlations. One possible source of this scatter is the variation in formation regions within the chromosphere. Indeed, similar levels of scatter in activity indices also based on the Ca ii H&K and infrared triplet lines have been reported in previous studies (e.g. Lanzafame et al. 2023). Nevertheless, this correlation validates the use of |$\Delta W_{\textrm {IRT}}$| as a reliable measure of stellar activity. Moreover, |$\Delta W_{\textrm {IRT}}$| measurements primarily reflect genuine activity differences between paired stars and are minimally affected by light contamination from close companions. This is supported by our ongoing LCOGT (Las Cumbres Observatory Global Telescope) monitoring campaign of 14 randomly selected C3PO pairs, which has revealed only one potential unseen companion among all monitored stars.

We acknowledge that one possible systematic error in measuring differential activity could arise from variations in calcium abundance ([Ca/H]) between paired stars (Carrera et al. 2007; Da Costa 2016). Noteworthy is that [Ca/H] reported in Liu et al. (2024) have been obtained using lines other than Ca ii infrared triplet or the Ca ii H&K discussed here. For instance, in Fig. 4, the object star (red curves) has a measured [Ca/H] value 0.08 dex (18 per cent) lower than its reference star (black curves), which could potentially result in shallower triplet lines. However, we find this effect to be negligible based on two key observations. First, the emission features in the Ca ii H&K line cores show substantial differences between the paired stars (Fig. 4c), with core flux variations reaching a factor of 2 – far too large to be explained by the small 0.08 dex difference in [Ca/H]. Second, the strong correlation between |$\Delta W_{\textrm {IRT}}$| and |$\Delta S_{\textrm {HK}}$| across our sample (Fig. 5) confirms that the observed variations primarily reflect differences in magnetic activity rather than calcium abundance.

3 DIFFERENTIAL ACTIVITY CORRELATES WITH ELEMENTAL ABUNDANCE DIFFERENCE

In Fig. 6, we present the relationship between the differential activity index, |$\Delta W_{\textrm {IRT}}$|⁠, and the slope of the condensation temperature versus elemental abundance difference, |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial T_{\rm {cond}}}$| (the same slope seen in Fig. 2), for 21 co-natal star pairs with chemical anomalies and available |$\Delta W_{\textrm {IRT}}$| measurements. The other co-natal pairs, which do not display chemical anomalies, have |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial T_{\rm {cond}}}$| values consistent with zero by definition and are therefore not shown. This figure reveals a correlation between |$\Delta W_{\textrm {IRT}}$| and |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial T_{\rm {cond}}}$| in these pairs, as indicated by the best-fitting linear model. Specifically, pairs with high ||$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial T_{\rm {cond}}}$|| values tend to exhibit strong differential activity. Conversely, four pairs with lower |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial T_{\rm {cond}}}$| values ( |$-0.3 < ({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial T_{\rm {cond}}} < 0.3$|⁠) also show low differential activity (⁠|$-0.7 < \Delta W_{\textrm {IRT}} < 0.7$|⁠), suggesting minimal or no exoplanetary impact. Among the 21 systems, four pairs identified by Liu et al. (2024) as potential cases of exoplanet engulfment (marked with green arrows), where observed elemental abundance differences match theoretical models of planetary engulfment, align with this trend.

Fig. 7 shows a more detailed analysis of the relation between |$\Delta$|[X/H] and |$\Delta W_{\textrm {IRT}}$| for four elements, namely carbon (C, volatile), sodium (Na, intermediate refractory), iron ii (Fe ii, refractory), and scandium ii (Sc ii, refractory). The condensation temperatures of these elements range from 40 to 1659 K, with their values annotated in Fig. 7. While no clear correlation is observed for C among the pairs with chemical anomalies, an inverse correlation is evident for Na, with a slope of |$-0.010 \pm 0.001$|⁠. This inverse correlation strengthens for Fe ii, with a slope of |$-0.018 \pm 0.001$|⁠, and is most pronounced for Sc ii, with a slope of |$-0.026 \pm 0.002$|⁠. These progressively stronger inverse correlations suggest a dependence of the |$\Delta W_{\textrm {IRT}}$|–|$\Delta$|[X/H] slope on condensation temperature for the chemical anomaly sample.

How does the differential activity of star pairs with homogeneous chemistry (i.e. |$|\rm{\Delta [\textrm {Fe/H]}}| < 3\sigma _{\rm{\Delta [\textrm {Fe/H]}}}$|⁠) compare to those with chemical anomalies? In Fig. 7, we show the differential activity for systems with homogeneous chemistry, represented by grey squares. As expected, these systems show no correlation between |$\Delta$|[X/H] and |$\Delta W_{\textrm {IRT}}$| due to their chemical homogeneity. Notably, differential activity among chemically homogeneous samples is concentrated around zero, while that of chemically anomalous samples (red circles) is more broadly and evenly distributed. Quantitatively, approximately 78 per cent of the chemically homogeneous pairs have |$\Delta W_{\textrm {IRT}}$| values between −2 and 2, with a standard deviation of 0.8. For the remaining 22 per cent, this increased differential activity may be due to phase differences in short-term (stellar rotation) and long-term (activity cycles) activity variations, even though stars within each pair generally have similar stellar properties, suggesting comparable activity levels.

We also perform a two-sample Kolmogorov–Smirnov (KS) test to evaluate whether the two |$\Delta W_{\textrm {IRT}}$| distributions shown in Fig. 7 are statistically similar. As the same sample of 21 stars is analysed for each element, the |$\Delta W_{\textrm {IRT}}$| distributions are identical across the four elements shown, except for C (top panel), where the sample consists of 20 stars. The test yields a p-value of 0.008, indicating a rejection of the null hypothesis at the 1 per cent significance level and confirming that the two |$\Delta W_{\textrm {IRT}}$| distributions are significantly different. Overall, Fig. 7 suggests that pairs without chemical anomalies tend to exhibit weaker differential activity, consistent with our finding that chemically anomalous pairs with potential planetary influences show stronger differential activity (see Fig. 6).

After examining elemental abundance differences, we now explore pairs with chemical anomalies using another metric, the differential abundance ratio, |$\Delta [X$|/O] = [X/O]|$_{\rm obj}$| − [X/O]|$_{\rm ref}$|⁠. Here, the differential represents the difference between the object and reference stars, and the abundance ratio is defined between element X and oxygen (O). We focus on |$\Delta [X$|/O] ratios, particularly |$\Delta$|[refractory/volatile], as they provide a direct comparison of elements that may be influenced differently by planetary processes, with refractory elements being more affected and volatile elements largely unaffected (e.g. Chambers 2010; Booth & Owen 2020). Similar to Fig. 7, we again analyse X = C, Na, Fe ii, and Sc ii in Fig. 8. We find that |$\Delta$|[C/O], which represents the ratio between two volatile elements (C and O), shows an insignificant correlation with differential activity. In contrast, the correlation becomes progressively stronger from |$\Delta$|[Na/O] to |$\Delta$|[Fe ii/O] and |$\Delta$|[Sc ii/O].

Expanding the correlation analysis to all 21 elements, Fig. 9 illustrates the significance of the inverse relationship between |$\Delta$|[X/H] and |$\Delta W_{\textrm {IRT}}$| (i.e. |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial \Delta W_{\textrm {IRT}}}$|⁠, represented as the slopes of the dashed lines shown in Fig. 7) as a function of condensation temperature. Systems with chemical anomalies (red circles) exhibit stronger anticorrelations for elements with higher condensation temperatures. That is to say, these anticorrelations are more pronounced for refractory elements (e.g. Ti and Al) compared to volatile elements (e.g. C and O). In contrast, the star pairs without chemical anomalies (black squares) show no strong correlations. Quantitatively, the slope of the red dashed line is |$\left(-1.89\pm 0.13\right) \times 10^{-5}$|⁠, while the slope of the black dashed line is |$\left(2.54\pm 1.19\right) \times 10^{-6}$|⁠, indicating that the correlation for systems with anomalies is much stronger.

Intriguingly, Earth-like elements such as C, O, Mg, Si, Ca, Ti, Cr, and Fe (e.g. Xu et al. 2014; Harrison et al. 2021), whose abundances are available in this work, follow the overall trend (see Fig. 7). This suggests that rocky planets or giant planets with rocky cores may be linked to the observed chemical anomalies.

Is the correlation with condensation temperature shown in Fig. 9 merely a variant of the correlation with atomic number, given that refractory elements are generally heavier than volatile elements? To check this, in Fig. 10, we examine the significance of the correlation between |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial \Delta W_{\textrm {IRT}}}$| and atomic number and calculate the slope of the red dashed line to be |$\left(-1.10\pm 0.22\right) \times 10^{-4}$|⁠. To properly compare this value with that of the red dashed line in Fig. 9, we rescaled the condensation temperature by dividing it by a factor of 28.6 (the maximum condensation temperature divided by the maximum atomic number) to match the scale of atomic number. After this adjustment, we find that the slope of the correlation between |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial \Delta W_{\textrm {IRT}}}$| and the rescaled atomic number is one-fifth of the slope of the correlation shown in Fig. 9. Additionally, the scatter (standard deviation) of the red circles around the best fit, represented by the red solid line in Fig. 10, is 0.004, which is twice as large as the scatter shown in Fig. 9. Therefore, the weaker and less tightly constrained correlation with atomic number compared to condensation temperature (Fig. 9 versus Fig. 10) suggests that condensation temperature is a stronger correlate of the observed abundance trends.

It is noteworthy that there is no clear correlation between |$\Delta W_{\textrm {IRT}}$| and |$\Delta \rm{T_{\text{eff}}}$| or between |$\Delta W_{\textrm {IRT}}$| and |$\Delta \rm{\log g}$| (see Fig. 11), indicating that these observed trends are not driven by subtle differences in stellar structure among the star pairs. To quantify this, we perform Kendall’s tau test, a non-parametric statistical method for assessing the correlation between two variables, with the null hypothesis being the absence of correlation. The test yields p-values of 0.571 for |$\Delta W_{\textrm {IRT}}$| versus |$\Delta \rm{T_{\rm{eff}}}$| and 0.432 for |$\Delta W_{\textrm {IRT}}$| versus |$\Delta \rm{\log g}$|⁠. Since these p-values exceed the commonly used significance thresholds of 0.01 and 0.05, the results indicate that we fail to reject the null hypothesis, concluding that there is no significant correlation in either case. Additionally, the inverse relationships between |$\Delta$|[X/H] and |$\Delta W_{\textrm {IRT}}$| remain significant when we restrict the analysis to star pairs with solar-like |$T_{\text{eff}}$| in the range |$5572 < \rm{T_{\rm{eff}}}< 5972 \, {\rm K}$|⁠.

4 DISCUSSIONS

4.1 Possible planet imprints suggested by differential activity

Liu et al. (2024) analysed elemental abundances in the C3PO sample of co-moving stellar pairs and found that at least 8 per cent of the pairs exhibit chemical anomalies. Their Bayesian analysis demonstrated consistency between the observed abundance differences and theoretical models of planetary engulfment, suggesting that planetary ingestion offers a better explanation for these anomalies compared to models of random abundance scatter or atomic diffusion effects.

In this work, we focus on differential magnetic activity directly calculated from stellar spectra. The inverse correlation between magnetic activity and chemical anomalies revealed in this study provides new insights into this phenomenon. While only 8 of the 21 chemically anomalous pairs show strong Bayesian evidence favouring planet engulfment models in Liu et al. (2024), we chose to analyse the complete set of chemically anomalous pairs for several reasons.

First, the conservative Bayesian evidence threshold used in Liu et al. (2024) was specifically calibrated through mock simulations to identify chemical signatures consistent with planet formation or ingestion of Earth-like composition planets. The choice of Earth-like composition, while somewhat arbitrary, reflects the limitation that Earth remains the only planet whose detailed chemical composition is well understood. Therefore, other pairs could still harbour planet-related signatures, albeit with different compositional patterns.

Second, since the C3PO program targets stellar twins specifically to mitigate stellar evolution effects, and given that ISM inhomogeneity is unlikely at scales below |$10^6$| au (the separation cut-off in this study), planet-related processes remain among the few viable explanations for the observed chemical anomalies. Importantly, our key findings, such as the trend shown in Fig. 6, remain robust even when considering only the pairs with high Bayesian evidence.

According to Guo (2023), stars that experience planetary engulfment exhibit shorter rotation periods compared to those without planets due to angular momentum transfer, with this difference increasing over time. While these stars maintain shorter rotation periods post-engulfment, for solar-mass host stars – such as those in our sample (median mass of 1.05 |${\rm M}_{\odot }$|⁠) – this rotational difference rapidly diminishes, independent of the planet’s mass or the star’s initial rotation period. Given the well-established correlation between rotation rate and stellar activity (e.g. Işık et al. 2023), planetary engulfment scenarios would predict either a weak or positive correlation between |$\Delta$|[X/H] and |$\Delta W_{\textrm {IRT}}$|⁠. However, our observation of a strong negative correlation appears inconsistent with this expectation.

On the other hand, exoplanet formation processes may deplete elemental abundances in host stars, as refractory elements can become sequestered during the formation of terrestrial planets, reducing their presence in the stellar atmosphere (e.g. Chambers 2010). Furthermore, close-in planets could enhance host star activity through tidal and magnetic interactions (Poppenhaeger & Wolk 2014; Kavanagh et al. 2021, also see the introduction for more references). This scenario predicts a negative correlation between |$\Delta W_{\textrm {IRT}}$| and |$\Delta$|[X/H], both of which are calculated using the same designation of object and reference stars within each pair, aligning with our observations.

Another mechanism in the exoplanet formation scenario that could enhance stellar activity involves shortened circumstellar disc lifetimes. Planets form within these discs during the first few million years of stellar evolution. During this early stage, the central star typically remains magnetically coupled to the disc through ‘disc locking’, preventing stellar spin-up through accretion and contraction toward the main sequence (Koenigl 1991; Collier Cameron & Campbell 1993; Edwards et al. 1993; Shu et al. 1994). However, Booth & Owen (2020) suggested that giant planets like Jupiter can create a gap in the gas disc, forming a pressure trap that prevents refractory-rich dust from accreting onto the host star. As a result, this gap can lead to the depletion of refractory elements in the star. Meanwhile, we expect that this gap decouples the outer disc from the host star, potentially enhancing stellar contraction during the PMS phase and shortening the disc’s lifetime. This aligns with the findings of Rosotti, Ercolano & Owen (2015), who showed that giant planets create gas gaps that reduce the accretion rate, leading to shorter disc lifetimes. Furthermore, theoretical studies by Gallet & Bouvier (2015) and Roquette et al. (2021) demonstrate that shorter disc lifetimes are associated with faster stellar rotation. Consequently, the host star may spin-up due to this enhanced contraction, resulting in increased rotation and activity during the main sequence (e.g. Monsch et al. 2023).

Given these considerations, we can hypothesize that signatures of exoplanet formation gradually weaken during a star’s main sequence lifetime. As shown in Fig. 12, systems with strong (⁠|$5\sigma$|⁠, where |$\sigma$| is the uncertainty in |$({\partial \Delta \left[X \operatorname{/}\rm {H} \right]})/{\partial T_{\rm {cond}}}$|⁠) correlations between |$T_{\rm {cond}}$| and |$\Delta$|[X/H] (cyan circles, 45 per cent of systems with chemical anomalies) exhibit a trend between stellar age and absolute differential activity, where younger stars show stronger activity differences. In contrast, systems without strong |$T_{\rm {cond}}$|–|$\Delta$|[X/H] correlations (dark purple circles in Fig. 12a) display a much weaker trend between |$\Delta W_{\textrm {IRT}}$| and age (light purple circles in Fig. 12b). These results hint at a connection between stellar age and the longevity of planetary chemical signatures. Given the relatively large uncertainties (1.6 Gyr) in age measurements, future studies focusing on acquiring additional data and improving age precision will be valuable for confirming and further refining these trends.

If the |$T_{\rm {cond}}$| signal is indeed due to planetary signatures – an assumption that warrants further verification by ruling out other possibilities (see next section for a discussion) – then planet formation, rather than planet ingestion, appears more likely to create both the observed chemical patterns (especially those with high Bayesian evidence) and the negative correlation with magnetic activity.

From a theoretical perspective, Booth & Owen (2020) argued that the 0.04 dex depletion in refractories observed in the Sun relative to solar twins due to planet formation can be readily explained by factors such as large disc masses (e.g. |$M_{\rm D} \approx 0.2 \, {\rm M}_\odot$|⁠), low photoevaporation rates, the shrinking convective zone after 2 Myr, and earlier planet formation, which collectively enhance dust depletion in the system. However, recent simulations in Hühn & Bitsch (2023) suggest that planet formation should not produce signals of the magnitude observed in C3PO (out of 21 pairs, two with |$\Delta \left[\text{Na/H}\right]$| and seven with |$\Delta \left[\rm{Sc\,\small{II}/H}\right]$| values ranging from 0.1 to 0.2 dex). Moreover, if the chemical patterns and rotation-induced magnetic activity arise from planet formation, maintaining such signals up to 4 Gyr requires further theoretical studies.

Resolving this theoretical discrepancy is beyond the scope of this observational paper, which aims to report the correlation between magnetic activity and chemical abundance patterns rather than arrive at theoretical explanations. Nevertheless, if planet-induced signals can persist for several Gyr, this may help explain the long-observed chemical anomalies of the Sun compared to other solar analogues (Bedell et al. 2018).

4.2 Alternative explanations

Here, we examine potential alternative mechanisms that could produce our observed signals. First, we consider stellar activity cycles. According to Spina et al. (2020), activity cycles affect elemental abundance measurements through magnetic field variations in the stellar atmosphere, which alter absorption-line equivalent widths. While such cycle-induced variations are real and could potentially complicate abundance analyses, Spina et al. (2020) demonstrated that these abundance variations show no correlation with condensation temperature. This is a crucial distinction from our observations, as the |$T_{\rm {cond}}$| correlation (Fig. 9) represents a key signature in our data set. This difference suggests that the correlation we detect between abundance variations and activity is not driven by the physics of activity cycles, consistent with the idea that exoplanets may be the underlying cause.

Internal stellar mixing processes present another possible mechanism that warrants careful consideration. In stars, atomic diffusion or gravitational settling causes heavier elements to sink while lighter ones rise (see Alecian & Deal 2023, for a recent review). Rotational mixing can counteract this process by transporting elements from the stellar interior to the surface, with faster rotation leading to more efficient mixing (e.g. Deal et al. 2020). This mechanism would predict that slower rotating stars show greater depletion of heavy elements, as reduced mixing efficiency would allow more settling to occur. However, our observations reveal the opposite trend: magnetically more active stars (which rotate faster) exhibit greater chemical depletion. Additionally, atomic diffusion effects are expected to produce systematic trends with stellar evolutionary state, particularly surface gravity. The C3PO survey’s focus on stellar twins specifically minimizes such evolutionary effects, as our pairs have nearly identical surface gravities. Indeed, we find (1) a weaker and less tightly constrained correlation of differential activity with atomic number compared to condensation temperature (Fig. 9 versus Fig. 10) and (2) no correlation between magnetic activity and stellar surface gravity in our sample, arguing against mixing processes as the primary driver of our observed chemical patterns. These multiple lines of evidence suggest that stellar mixing processes alone cannot explain the observed relationship between magnetic activity and chemical abundances.

Another widely discussed factor influencing abundance trends with condensation temperature is galactic chemical evolution (GCE), as abundance patterns can be influenced by variations in stellar ages (e.g. Bedell et al. 2018; Sun et al. 2025b). However, Spina, Meléndez & Ramírez (2016a) suggested that GCE effects could be disentangled from planetary signatures by removing the [X/Fe] versus age trends. Furthermore, other studies indicate that while GCE can account for some of the observed variations, it cannot fully explain the entire trends observed in many stars (Nissen 2016; Spina et al. 2016b; Bedell et al. 2018). For example, the Sun still exhibits refractory depletion even after GCE corrections are applied (Nissen 2016; Spina et al. 2016b). Therefore, GCE is likely a contributing factor but not the primary driver of these chemical anomalies. The co-natal stars in our study further minimize the impact of GCE because of their shared formation history.

5 CONCLUSIONS

We investigated the differential activity of 21 co-natal star pairs with chemical anomalies, selected from the C3PO program, which provides the largest and most homogeneous sample of co-moving systems to date. Each pair consists of a reference star and a companion star with similar fundamental properties (see Section 2.1). The typical (median) age of these stars is approximately 3.2 Gyr, with a standard deviation of 3.5 Gyr. Using high-resolution, high-S/N spectra from the Magellan and Keck, we derived a differential activity index, |$\Delta W_{\textrm {IRT}}$|⁠, based on the flux difference in the Ca ii infrared triplet lines between the two stars in each pair, which is validated with |$\Delta S_{\textrm {HK}}$| (see Fig. 5 and Section 2.2). Using elemental abundance measurements from Liu et al. (2024), we then examined the relationships between differential activity (⁠|$\Delta W_{\textrm {IRT}}$|⁠), elemental abundance differences (⁠|$\Delta$|[X/H]), and dust condensation temperature. Our main findings are as follows:

• Pairs exhibiting strong correlations between condensation temperature and elemental abundance differences also show pronounced differential activity (Fig. 6).

• We discover significant inverse correlations between |$\Delta W_{\textrm {IRT}}$| and |$\Delta$|[X/H] for refractory elements (Na, Fe ii, and Sc ii), indicating that chemically depleted stars are magnetically more active. No such correlation exists for volatile elements (Fig. 7).

• The strength of these inverse correlations increases systematically with condensation temperature (Fig. 9).

• In systems with strong correlations between condensation temperature and elemental abundance difference, younger stars exhibit larger differential activity, suggesting an age-dependent effect (Fig. 12).

These patterns add a new dimension to understanding the origin of stellar chemical anomalies. If the chemical signatures are indeed planet related, then the inverse correlation between chemical depletion and magnetic activity, coupled with the condensation temperature dependence, points toward planet formation rather than engulfment as a possible mechanism (see Section 4). Booth & Owen (2020) found that the 0.04 dex depletion observed in the Sun relative to solar twins due to planet formation can be readily explained. However, Hühn & Bitsch (2023) predict chemical signatures substantially smaller than those we observe. Understanding this discrepancy and the longevity of these signals requires further investigation.

Several follow-up studies have the potential to provide more insights. First, radial velocity monitoring, including multi-epoch Gaia RVS spectra from upcoming data releases (e.g. Gaia Data Release 4), will provide valuable data to investigate whether these bright stars host giant planets. This is particularly relevant for pairs with significant differential abundances and activity (Fig. 7) and offers a potential test of the Booth & Owen (2020) hypothesis. Second, rotation period measurements from TESS (Transiting Exoplanet Survey Satellite) light curves (e.g. Colman et al. 2024) could illuminate the connection between differential activity and rotation, though careful consideration of binary contamination effects will be necessary. Finally, theoretical work is needed to reconcile the magnitude of observed chemical signatures with planet formation models. These combined efforts will advance our understanding of how planets shape the chemical and magnetic evolution of their host stars.

ACKNOWLEDGEMENTS

JY acknowledges the Joint Research Fund in Astronomy (U2031203) under a cooperative agreement between the National Natural Science Foundation of China (NSFC) and the Chinese Academy of Sciences (CAS). Y-ST was supported by the National Science Foundation under grant no. AST-2406729. T.N. acknowledges support from the Knut and Alice Wallenberg Foundation. This paper includes data gathered through the C3PO programme with the 6.5-m Magellan Telescope located at Las Campanas Observatory, Chile. Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The observatory was made possible by the generous financial support of the W. M. Keck Foundation. Based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 108.22EC.001. 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.

DATA AVAILABILITY

The spectra used to calculate activity indices in this work were collected with Magellan and Keck and processed by Yong et al. (2023) and Liu et al. (2024). Stellar parameters and elemental abundances were adopted from Liu et al. (2024). To facilitate the reproduction of the main results presented in this work, the data on activity indices will be made available upon reasonable request to the corresponding author.

Footnotes

REFERENCES