Evidence of Rocky Planet Engulfment in the Wide Binary System HIP 71726/HIP 71737

Binary stars are supposed to be chemically homogeneous, as they are born from the same molecular cloud. However, high-precision chemical abundances show that some binary systems display chemical differences between the components, which could be due to planet engulfment. In this work, we determine precise fundamental parameters and chemical abundances for the binary system HIP 71726/HIP 71737. Our results show that the pair is truly conatal, coeval, and comoving. We also find that the component HIP 71726 is more metal-rich than HIP 71737 in the refractory elements such as iron, with Δ[Fe/H] = 0.11 ± 0.01 dex. Moreover, HIP 71726 has a lithium abundance 1.03 dex higher than HIP 71737, which is the largest difference in Li detected in twin-star binary systems with ΔT eff ≤ 50 K. The ingestion of 9.8−1.6+2.0 M ⊕ of rocky material fully explains both the enhancement in refractory elements and the high Li content observed in HIP 71726, thereby reinforcing the planet-engulfment scenario in some binary systems.


Introduction
Chemical tagging is one of the main goals of large surveys such as GALAH (De Silva et al. 2015), APOGEE (Majewski et al. 2017), Gaia-ESO (Gilmore et al. 2012), etc. As per this method, we could identify the birthplaces of stars in the Galaxy using simply their chemical composition. The technique is mainly based on two assumptions: (i) the chemical homogeneity of the proto-cloud from which stars form; (ii) the unique chemical signatures of formed stars compared to other star groups. In this context, both open clusters and wide binaries are ideal laboratories for testing these assumptions, since they are considered coeval and conatal stars (e.g., De Silva et al. 2007;Ting et al. 2012;Andrews et al. 2017;Ness et al. 2018;Hawkins et al. 2020).
Despite the large number of wide binary systems identified in the past decade (e.g., Tokovinin 2014; Andrews et al. 2017;El-Badry et al. 2021), only a few works have dedicated to studying their chemical homogeneity using spectra of high resolving power. In this sense, Laws & Gonzalez (2001) and Gratton et al. (2001) were pioneers in determining chemical abundances of wide binaries using the differential analysis and spectra of very high resolving power (R = 58,000-160,000) and signal-to-noise ratio (S/N = 100-400). Other studies have focused on analyzing a few elements such as iron, vanadium, and lithium (Martín et al. 2002;Desidera et al. 2004Desidera et al. , 2006. Most of these wide binaries are chemically homogeneous, i.e., the metallicity difference between the components (Δ[Fe/H]) is ∼0.0 dex; however, other pairs have Δ[Fe/H] ∼0.1 dex.
A more comprehensive chemical analysis of a large sample of wide binaries was carried out by Andrews et al. (2019), who used APOGEE (R = 22,500) spectra to determine the abundance of 14 elements. They found that the chemical abundances of the systems are consistent to the level of 0.1 dex; however, the precision decreases significantly to ∼0.3 dex for some elements like S (see their Figure 10). Hawkins et al. (2020), using spectra of higher resolving power (R = 60,000), obtained a better precision and showed that the metallicity differences between the components of their sample of wide binaries is centered at [Fe/H] ∼0.0 dex with a scatter of 0.05 dex, and extreme values from about −0.1 to +0.1 dex. There is still a debate concerning whether wide binaries with large separations are still bound. According to Jiang & Tremaine (2010), binary systems with separations  10 5 au are unbound; however, Kamdar et al. (2019) demonstrated through simulations that comoving pairs with 3D differential velocities (Δv 3D ) < 2 km s −1 are probably bound. This is supported by Nelson et al. (2021), who found that chemical abundances of comoving pairs with large separations are also consistent to the level of ∼0.1 dex. More recently, Espinoza-Rojas et al. (2021) showed that binary systems, even those with some metallicity difference, are useful to validate chemical clocks. All these results suggest that wide binaries can be excellent laboratories for testing the chemical-tagging hypotheses. However, it will only be possible if spectra of high resolving power and S/N are employed in the differential analysis to achieve high-precision chemical abundances.
When one of the components in a binary system is richer in refractory elements than its companion, the effect can be directly attributed to engulfment of planetary material (Pinsonneault et al. 2001). This scenario was suggested for several binary systems with Δ[Fe/H] 0.02 dex such as HD 219542, 16 Cygni, XO-2, WASP-94, and HAT-P-4 (Gratton et al. 2001;Ramírez et al. 2011;Biazzo et al. 2015;Ramírez et al. 2015;Teske et al. 2016aTeske et al. , 2016bSaffe et al. 2017).
Most twin-star binaries analyzed with high resolving power (R = λ/Δλ > 60,000) have Δ[Fe/H] up to ∼0.05 dex. The system HIP 34407/HIP 34426 (Ramírez et al. 2019; see also Nagar et al. 2020;Espinoza-Rojas et al. 2021) and HD 240430/HD 240429 (Oh et al. 2018) is the exception with Δ[Fe/H] ∼ 0.2 dex. The Fe-rich component in these systems is enhanced in refractory elements relative to its companion. Notably, the Li abundance (A(Li)) in HD 240430 is larger than HD 240429 by ∼0.5 dex, likely due to a planet accretion of 15 M ⊕ of rocky material. In particular, this scenario had already been studied by Gratton et al. (2001), who analyzed the system HD 219542A-B (Δ[Fe/H] = 0.09 dex) and reported a very high enhancement in Li abundance (ΔLi = 1.35 dex) in HD 219542A.
Based on the results discussed above, the components of binary systems with Δ[Fe/H] 0.02 dex (or somewhat higher, depending on the achieved precision) can be affected by mechanisms that change their original composition, thereby questioning the viability of chemical tagging for these pairs. In this work, we report the discovery of abundance anomalies in the wide binary system HIP 71726/HIP 71737 (hereafter HIP 71726-37), whose components of G-type stars show an unusual differential chemical abundance that suggest a rocky planet engulfment by the Fe-rich component (HIP 71726).

Sample Selection and Observations
We selected our pair from the sample of wide binaries of El-Badry et al. (2021), after applying the color constraints given in Yana Galarza et al. (2021), Table 1, which were specially established for searching solar twins. We estimated the projected separation in ∼16764 au using the ASTROPY package, which employs coordinates and distances as input parameters. The coordinates and parallaxes are adopted from Gaia EDR3 (Gaia Collaboration et al. 2021). The distance (∼55 pc) is calculated inverting the parallax. We also determined the Galactic orbits for HIP 71726-37 using the GALA 4 code (Price- 2 (e.g., Ramírez et al. 2007, see their Figure 3). The 3D velocity difference (Δv 3D ) is 0.34 ± 0.07 km s −1 , which takes into account the space velocities of each component. These results indicate that the binary system is bound and truly comoving, since the components also share identical radial velocities (see Table 1).
The spectra of our comoving pair were obtained using the Robert G. Tull Coudé Spectrograph (TS23; Tull et al. 1995) on the 2.7 m Harlam J. Smith Telescope at the McDonald Observatory during the months of 2021 February and March. The solar spectra were obtained through the reflected light from the Moon. We reduced the TS23 spectra employing our own scripts 5 based on PYRAF, which perform the standard procedure, i.e., creation of the master flat, flat-field correction, sky subtraction, order extraction, etc. The resulting spectra have high resolving power R = 60,000 with S/N = 450 for the Sun and S/N ∼350 for each component of the binary system.

Fundamental Parameters
The radial/barycentric velocity correction and normalization of the TS23 spectra were carried out using our own semiautomatic PYTHON scripts described in Yana Galarza et al. (2021). The equivalent widths (EWs) were measured manually through Gaussian fits to the line profile using the KAPTEYN kmpfit Package (Terlouw & Vogelaar 2015), which is also implemented in our scripts. The technique is based on line-by-line measurements choosing pseudocontinuum regions of 6 Å in order to achieve a high precision. We adopted the line list of Meléndez et al. (2014), which includes information about the excitation potentials, oscillator strengths, laboratory gf log values, and hyperfine structure corrections (McWilliam 1998;Cohen et al. 2003;Asplund et al. 2009). The adopted line list is composed of weak-to-moderate-strength lines (EW < 130 mÅ).
We determined the effective temperature (T eff ), surface gravity ( g log ), metallicity ([Fe/H]), and microturbulence velocity (v t ) using the automatic q 2 ) code adopting the spectroscopic equilibrium technique. The q 2 was configured to employ the Kurucz ODFNEW model atmospheres (Castelli & Kurucz 2003) and the 2019 local thermodynamic equilibrium (LTE) code MOOG (Sneden 1973). The results are summarized in Table 1, where the first two rows are the stellar parameters of the pair relative to the Sun and the last two are the stellar parameters of HIP 71737 and the Sun relative to the primary component (HIP 71726). Both differential methods result in identical stellar parameters for HIP 71737, and in the case of the Sun, we fully recover the solar stellar parameters (T eff = 5777 K, We also computed the ages and masses employing the Yonsei-Yale isochrones of stellar evolution (Yi et al. 2001;Demarque et al. 2004), which is already implemented in the q 2   (3)). b Stellar parameters for the binary system HIP 71726-37 determined through the differential technique relative to the Sun. c Stellar parameters for HIP 71737 and the Sun, both relative to the primary component HIP 71726.
code. Unlike the traditional isochronal age calculations that use only absolute V magnitude or spectroscopic g log as input parameters, we use both parameters in the q 2 code. The parallaxes are adopted from Gaia EDR3, and we corrected them by subtracting −15 ±18 μas as suggested by Stassun & Torres (2021).
Additionally, we estimated the activity indices (S HK,TS23 ) measuring the flux in the Ca II H and K lines following the prescription given in Wright et al. (2004). The S HK,TS23 was converted to the Mount Wilson system (S MW ) using Equation (4) of Yana Galarza et al. (2021). Then, we determined the chromospheric indices ( ¢ R log HK ) following the procedure of Noyes et al. (1984). Finally, the chromospheric ages are estimated from the age-mass-metallicity activity correlation derived by Lorenzo-Oliveira et al. (2016, see their Equation (1)), which takes mass and metallicity biases into account. The results are shown in Table 2.
Both the isochronal and chromospheric ages indicate a difference of ∼1 Gyr between the components, which could signify that the system is not coeval. However, if the Fe-rich component (HIP 71726) was enhanced with refractory elements, its iron content does not reflect the metallicity of the entire star, but only that of the convective zone, which artificially changes the age because an overestimated metallicity is being used to obtain isochronal and chromospheric ages. Therefore, to assess the coevality of our system, we performed a simple experiment based on estimating ages and masses using the mean and the lowest [Fe/H] of the pair, assuming that these are the "real" metallicities. As shown in Table 2, the isochronal ages of the components are identical between them in both cases, with a difference of only 0.1 Gyr. Similar results are found for the chromospheric ages. Therefore, we can consider that our binary system is truly coeval.

Chemical Composition
The chemical abundances for 26 elements (C, O, Na, Mg, Al, Si, S, K, Ca, Sc, Ti, V, Cr, Mn, Co, Ni, Cu, Zn, Sr, Y II, Ba II, La II, Ce II, Nd II, Sm II, and Eu II) were estimated using the q 2 code, which was configured to perform corrections of NLTE for the O I triplet (Ramírez et al. 2007) and takes into account the hyperfine structure for V, Mn, Co, Cu, Y II, Ba II, La II, and Eu II . The internal precision for most differential abundances is ∼0.01 dex, which was calculated following the prescription given in Epstein et al. (2010). It is important to highlight that there is only one spectral line available for the elements K I, Sr I, and Eu II; therefore, their abundances were estimated assuming various pseudocontinuum levels.The adopted abundance for each of these three elements is the average of all the measurements while the observational uncertainty is the standard deviation. Special attention was paid to regions where telluric lines can affect the measurements (e.g., K I at 7698.97 Å). The differential chemical abundances for the components of the pair relative to the Sun are summarized in Table 3. We also estimated the elemental abundance differences between HIP 717237 and HIP 717226 in order to improve the precision. This process was performed adopting HIP 717226 as a reference. The result is shown in Figure 1, where it can be clearly seen that the differential abundance pattern of HIP 71737 is deficient in refractory elements (with high condensation temperature) relative to HIP 71726. Carlos et al. (2016Carlos et al. ( , 2019 showed that the surface Li abundance (A(Li)) in solar twins decreases with age, but with some scatter, probably due to somewhat different stellar masses and metallicities. We determined the A(Li) for the binary system strictly following the procedure given in Yana Galarza et al.  Table 3, while in Figure 2 the Li lines are depicted in the observed spectra of HIP 71726-37. Conspicuously, the Fe-rich star has A(Li)1.03 dex higher than HIP 71737. To our knowledge, this is the largest difference in A(Li) observed in twin-star wide binaries with ΔT eff 50 K. The binary HD 219542A-B analyzed by Gratton et al. (2001) has a larger difference in Li (1.35 dex), but it has ΔT eff ∼300 K.
It is very well known that A(Li) increases with T eff (e.g., Ramírez et al. 2012;Bragaglia et al. 2018), and also its depletion is higher for stars with larger metallicities (Castro et al. 2009;Ramírez et al. 2012;Carlos et al. 2019, see their Figure 5 for thin-disk stars with [Fe/H] > 0). Based on this evidence, the Fe-rich star HIP 71726 should have depleted more Li than its companion; however, this is not the case. According to the observational trend of Li with mass by Carlos et al. (2019), the different masses (0.09 M e ) of the pair suggest a difference of 0.3 ± 0.1 dex in A(Li) (see also Ramírez et al. 2012;Delgado Mena et al. 2015;Castro et al. 2016), which only explains in part the observed Li in HIP 71726. Another possibility is different amounts of rotational mixing between  Santos et al. 2016;Yana Galarza et al. 2016), we found similar projected v i sin velocities for the components (see Table 1). The only plausible scenario that remains to explain the Li difference is the ingestion of rocky material by the Fe-rich star.

Planet Engulfment
The abundance patterns (versus condensation temperature) found by Meléndez et al. (2009) (2015) in solar twins and twin-star binary systems have been linked to the formation of planets (for an alternative explanation regarding dust cleansing, see Gustafsson 2018aGustafsson , 2018b. In particular, it is argued that the formation of rocky planets removes refractory elements from the hosted star. Thus, in a binary system, the component with more rocky planets is likely to be deficient in refractory elements. Based on this hypothesis, the deficiency of refractory elements of HIP 71737 would result from host rocky planets (see Figure 1), although no planets have been detected yet in this system.
An alternative to the metal depletion scenario is a rocky planet engulfment, where the Fe-rich star HIP 71726 would have ingested planetary material rich in refractory elements. To validate this scenario, we used our code TERRA 7 , which is updated with the isochrones of stellar evolution of Yi et al. (2001) and the more recent solar abundances of Asplund et al. (2021). In brief, the methodology consisted in first determining the convective mass of a star through interpolations in metallicity and mass. Then, the code estimates the amount of rocky material necessary to reproduce the observed abundance pattern by adding meteoritic (Wasson & Kallemeyn 1988) and terrestrial (Allègre et al. 2001) abundances patterns into its convective zone (for more details see the appendix of Galarza et al. 2016).
The observed differential abundance pattern of HIP 71726 can be explained by a rocky planet engulfment of -+ 9.8 1.6 2.0 M ⊕ , which is a mixture of 7.0 M ⊕ of chondrite-like composition and 2.8 M ⊕ of Earth-like composition. As can be seen in Figure 3, both the observed (open squares) and predicted (red circles) differential abundances of HIP 71726 are in good agreement.
The A(Li) is also fully predicted by the planetary engulfment (see Figure 3); however, special attention must be paid to this point, since determining the total depleted Li abundance is a challenging task because we do not have information about the initial Li content at the time of accretion or the efficiency at If the planet-engulfment scenario discussed above is correct, HIP 71726ʼs rotation would be slightly faster (Oetjens et al. 2020). However, there are no significant differences in the v i sin of the components (see Table 1). This result does not invalidate our scenario because the planetary engulfment could happen when the star was young and its rotation was high, thus the small increase in angular momentum added by the planet has already been lost. Alternatively, the mass of the planet is too small to cause a significant increase in the star's angular momentum. In addition, according to dos Santos et al. (2016), a typical v i sin for a star with ∼4.9 Gyr is ∼1.8 km s −1 (e.g., HIP 1954 therein). Therefore, the very low values of the v i sin of the components (0.17 km s −1 and 0.10 km s −1 ) could be only a projection effect (ĩ sin 0). This means that the components could have different rotational velocities, but due to their inclinations, the difference cannot be detected.

Summary and Conclusions
The precise stellar parameters and chemical abundances for the binary system HIP 71726-37 were estimated through a differential analysis relative to the Sun. To increase our precision, we also determined the stellar parameters of HIP 71737 and the Sun relative to the primary component (HIP 71726).
The binary system has Δ[Fe/H] ∼0.1 dex, which means that the pair is chemically inhomogeneous. We estimated the Galactic orbits for the system and found identical eccentricity, perigalacticon, and apogalacticon and maximum excursion vertical values for both components. The 3D velocity difference is 0.96 km s −1 , suggesting that the pair is truly comoving. The ages and masses were calculated using isochrones and the age-activity relation of Lorenzo-Oliveira et al. (2016). The slightly different ages (∼1 Gyr) of the components seemed to suggest that they are not coeval. However, when we use the mean or the lower [Fe/H] of the pair, an identical age is obtained (see Table 2). This suggests that HIP 71726 may have been contaminated by some material rich in refractory elements. Therefore, based on the above results, we conclude that the binary system HIP 71726-37 is truly conatal, coeval, and comoving.
The chemical abundance pattern (versus condensation temperature) of HIP 71737 is deficient in refractory elements relative to the Fe-rich HIP 71726 (see Figure 1). Following Meléndez et al. (2009), HIP 71737 could have formed rocky exoplanets. However, no planets have been detected in this system yet. An alternative explanation is a planet-engulfment