© The Authors 2023

1. Introduction

The enhanced depletion of lithium in the Sun, discovered almost 70 years ago (Greenstein & Richardson 1951), remains the epitome of the solar Li puzzle. The surface Li¹ abundance in the Sun (A(Li) = 1.05 ± 0.10 dex, Asplund et al. 2009) is about 140 times less than the initial protosolar abundance, which is assumed to be the meteoritic value (A(Li) ∼ 3.30, e.g., Lodders 2019). Nevertheless, we note that very recently, on the basis of a new 3D non-thermodynamic equilibrium analysis, Wang et al. (2021) have found a Solar surface Li abundance of A(Li) = 0.96 ± 0.05 dex. However, the base of the surface convective layer of the Sun is not hot enough for nuclear reactions to destroy Li, which is destroyed by proton and alpha capture reactions at temperatures higher than ∼2.5 × 10⁶ K. This Li abundance anomaly is not unique to the Sun; large dispersion in Li abundances is observed in solar-type main-sequence (MS) stars of similar age, mass, and metallicity (see, e.g., Lambert & Reddy 2004), which is difficult to explain within the standard models of stellar evolution. Canonical stellar evolution predicts MS stars with Li-depleted interiors (see for example D’Antona & Mazzitelli 1994; Piau & Turck-Chièze 2002; D’Antona & Montalbán 2003; Montalbán & D’Antona 2006), but almost no surface depletion during the MS, nor the dispersion observed in MS stars with very similar stellar parameters. As a consequence, some nonstandard mixing mechanisms have been proposed to explain this observational result such as rotation, mass loss, diffusion, and/or gravity waves (e.g., Pinsonneault et al. 1992; Swenson & Faulkner 1992; Michaud 1986; Umezu & Saio 2000; Charbonnel & Talon 2005), but a satisfactory solution has not yet been found. This has opened up the possibility of an external origin as the cause for this depletion. Tidal interaction with sub-stellar objects such as brown dwarfs (BD) and/or planets, or even their engulfment, are some of the most frequent alternative hypotheses. This latter scenario would not be uncommon, as planetary engulfment seems to be the fate awaiting a significant number of planetary systems. In fact, transit surveys indicate that hot Jupiters (HJs)² are found orbiting ∼1% of Sun-like stars (e.g., Winn & Fabrycky 2015).

Alexander (1967) first suggested that the Li depletion observed in exoplanets’ host stars could be directly related to the mechanism of formation and/or migration of planets in the corresponding stellar system. This idea has been explored in many observational studies (see e.g., King et al. 1997; Gonzalez 2000, 2015; Chen & Zhao 2006; Israelian et al. 2009; Gonzalez et al. 2010; Delgado Mena et al. 2015). The potential discrepancy between lithium abundances in MS stars that host or don’t host exoplanets is, however, controversial and still divides the observer community. Most of the arguments against the existence of this discrepancy rely on a wrong interpretation of observational data or on an observational bias (e.g., Ryan 2000; Luck & Heiter 2006; Baumann et al. 2010). The reason for this is, in part, the difficulty in finding a large enough homogeneous sample of MS stars having equal properties, such as equal age, mass, metallicity, and rotational velocity. Furthermore, the derived Li abundances in many of these observational studies are not accurate enough to extract any solid conclusion.

A way to partially overcome this issue is carrying out detailed differential chemical analyses in solar-like twin stars in wide binary systems under the hypothesis that they were formed at the same time, in the same birth cloud, and with a nearly identical composition. Nowadays, this type of analysis is able to determine chemical abundances with unprecedented (relative) precision (≲0.01 dex for some elements, see e.g., Oh et al. 2018; Liu et al. 2018; Maia et al. 2019; Nagar et al. 2020). Unfortunately, the planet population around most of these solar-like twin stars is unknown, and in the few cases where planets have been detected, these may orbit or not (apparently) a star showing a chemical anomaly. As far as we know, only in three of the solar-like binary systems studied so far has Li been measured accurately enough: the pairs 16 Cyg A-B, Kronos-Krios, and HIP 71726–HIP 71737 (Galarza et al. 2021; Maia et al. 2019; Oh et al. 2018). Surprisingly, Li is found to be more abundant in the most metal-rich star of the pair (for instance, in 16 Cyg A with respect to B), the star which presumably accreted rocky material (i.e., a sub-stellar object). This seems to contradict the conclusion from the Li surveys in exoplanets’ host stars mentioned above, which suggests an extra Li depletion in the host star. Therefore, from the observational point of view, the impact of the engulfment of a sub-stellar object by an MS star on the surface Li abundance remains unclear.

On the theoretical side, planetary accretion may occur from the pre-main sequence through the white dwarf stage (see e.g., Kaiser et al. 2021; Kruckow et al. 2021), and different physical processes drive the accretion at different stages (Jackson & Carlberg 2018). Nevertheless, only a few studies exist in the literature on the interaction between an MS star and a sub-stellar object (see e.g., Metzger et al. 2012; Yamazaki et al. 2017; Jia & Spruit 2018), since the overwhelming majority of them are devoted to the interaction occurring during the cool giant phases of the stellar evolution, in which the probability of planetary engulfment increases (Siess & Livio 1999a,b; Villaver et al. 2014; Privitera et al. 2016a,b; Aguilera-Gómez et al. 2016, 2020; Staff et al. 2016; Jackson & Carlberg 2018; MacLeod et al. 2018; Soares-Furtado et al. 2021). Most of these theoretical studies are 1D analytical simulations and usually assume that: (a) the accreted sub-stellar object or rocky-material dissolves completely in the stellar envelope at the point where the Virial temperature of the accreted material equals that in the envelope (Siess & Livio 1999a); and (b) the engulfment does not significantly alter the internal structure and evolution of the host star. Obviously, the former hypothesis is critical for evaluating the impact on the surface chemical composition of the host star. On this basis, Montalbán & Rebolo (2002) first studied planet accretion onto solar-type stars and showed that it may significantly change the surface abundance of Li if it takes place during the MS evolutionary stage. Nevertheless, they concluded that there is a narrow range of stellar masses and metallicities where Li can be efficiently preserved in the envelope. This was contradicted in part by Baraffe & Chabrier (2010) who found that episodic accretion onto young solar-type stars can produce objects with significantly higher central temperatures than non-accreting stars of the same mass and age. As a consequence, Li depletion is larger in these objects. This could explain the extra Li depletion observed in some young stars belonging to open clusters. Later on, Théado & Vauclair (2012) and Deal et al. (2015) showed that the accretion of rocky material onto the host star creates an inverse gradient of molecular weight, which leads to thermohaline convection. As a consequence, heavy elements are mixed downward until the mean molecular weight gradient becomes nearly flat. In most cases, a tiny or no signature of the accreted heavy elements remains at the surface after a few Ma. They concluded that the accretion of heavy matter cannot lead to any increase in the Li abundance at the surface, but conversely, it may lead to a reduction.

To our knowledge, the sole 3D studies on the accretion of a sub-stellar object by an MS star are those by Sandquist et al. (1998, 2002). These authors performed hydrodynamic simulations to study the giant-planet (Jupiter and Saturn-like planets) consumption by a solar-like host star and, as a consequence, its possible Li pollution. They showed that depending on the details of the planetary interior, the partial or total dissolution of giant planets can result in significant enhancements in the metallicity of the host star with mass in the range of 1.0 ≲ M/M_⊙ ≲ 1.3. Nevertheless, their calculations had a very limited spatial resolution and did not simulate the full structure of the MS star. Additionally, their simulations lacked a mechanism to transport angular momentum from the planet to the star’s envelope.

From the discussion above it follows that both observational and theoretical studies are far from being conclusive and many questions still remain open in the engulfment scenario. These questions include the fate of the sub-stellar object, and whether it survives or if it is destroyed by ablation or tidal disruption. And if some sub-stellar objects are destroyed, the fraction of those that are is not known. It is further unknown if the interaction leads to an alteration of the MS star properties, such as spin-up, mass loss, luminosity, and size, or a change in the surface abundances, in particular, that of Li.

For the present study, we performed 3D hydrodynamical simulations of the engulfment of a BD into an MS star of ∼1 M_⊙ and solar composition by using the smoothed particle hydrodynamics (SPH) method. Our main aim is to answer some of the questions above mainly concerning the final fate of the sub-stellar object and the consequences that the engulfment may have on the host star in terms of its mass loss and rotation, and in particular on the possible change in the surface Li abundance. We focus our study on the relevance of two parameters: the mass of the BD and the type of interaction. In a forthcoming study, we will address the same issues but we will change the evolutionary status of the host star to a red giant branch and/or a horizontal branch star. The structure of the paper is the following. In Sect. 2 we describe the numerical setup and the simulated engulfment scenarios together with the main approximations that were made. Section 3 discusses the results of these simulations, particularly the impact that the engulfment provokes on the global properties of the hosting star and on its surface Li abundance. Finally, in Sect. 4 we discuss the main results, present our conclusions, and suggest possible ways for this study to be continued.

2. Numerical setting

Traditionally, SPH codes have been used to simulate planetary impacts and stellar collisions (Slattery et al. 1992; Sills et al. 1997) due to their intrinsic conservation properties and ability to follow complex geometries. To simulate the scenarios presented in this study, we used the hydrodynamic code SPHYNX (Cabezón et al. 2017). SPHYNX³ is a state-of-the-art SPH code that includes most of the latest upgrades in the SPH technique, including the Integral Approach to Derivatives (García-Senz et al. 2012; Cabezón et al. 2012; Rosswog 2015; Valdarnini 2016), the sin c family of interpolating kernels (Cabezón et al. 2008), artificial viscosity switches (Cullen & Dehnen 2010; Read & Hayfield 2012), and generalized volume elements (Saitoh & Makino 2013; Hopkins 2013; García-Senz et al. 2022). It also includes a two-level dynamic load-balancing implementation that considerably reduces the load imbalance that might appear in a distributed calculation (Mohammed et al. 2020). Overall, this code is well suited to coping with the highly dynamic, geometrically distorted evolution of stellar collisions, and the development of hydrodynamical instabilities, even in the subsonic regime. The ability of the SPH technique to handle free surfaces is also particularly suitable for this scenario, where we can simulate both full stars and follow the fluid that may escape the system. Additionally, the spatial resolution is adaptive (which in SPH translates to variable smoothing length), increasing where density increases and reaching resolutions as high as 10⁸ cm (i.e., ∼1.4 × 10⁻³ R_⊙) in our configuration. Finally, SPH codes have excellent energy and momentum conservation properties, which is of critical importance when evaluating angular momentum transfer in dynamic scenarios such as those studied here.

2.1. Simulated scenarios

In Table 1 we list the scenarios that we simulated. These are divided into two groups depending on the mass of the BD: the A models correspond to M_BD = 0.019 M_⊙, while the B models use M_BD = 0.01 M_⊙. With the former mass choice, we wanted to maximize the possible effects of the interaction of a sub-stellar object on the host star and, with this relatively low BD mass, we ensured that the original Li content in the BD was fully preserved. We note that previous 1D studies – although using simplifying assumptions – revealed the existence of a critical mass below which the colliding sub-stellar object is evaporated within the envelope. This mass limit was found to be ∼0.02 M_⊙ (e.g., Livio & Soker 1984); hence, we are in the position to test this result in a 3D setup. The second choice is close to the theoretical mass limit threshold separating a planet from a BD, ∼13 M_J (see e.g., Spiegel et al. 2011). On the other hand, the idea behind using the lower-mass BD model is to get closer to the planet engulfment scenario but avoiding the need to simulate rocky structures, which may imply changing the equation of state. For both models (A and B), we explored three types of interactions: a head-on collision, a grazing collision, and a merger. The first two are more violent, have shorter timescales, and are less realistic than the latter interaction, due to their extremely low cross section. Therefore, we expect to see significant differences in the outcome of these interactions. The 00 and 05 models in Table 1 are collisions: head-on and grazing, respectively. The M models, which probably represent a more realistic scenario, are mergers induced via setting the BD in an orbit around the MS star with an orbital velocity artificially reduced by 20% to induce the merger⁴. The last column of Table 1 shows the approximate amount of CPUh employed for each calculation. These values are purely indicative. Despite all models being calculated with 224 cores, the simulations were performed in different facilities with different generations of CPUs and network interconnections. The wall-clock time for each simulation was 30 − 40 days for the collisions and 100 days for the mergers.

Figure 1 shows a snapshot of the initial positions of the stars for each different interaction. We note in Table 1 that the difference in the number of particles (i.e., SPH fluid elements) between the BD and the MS stars in the A models is a factor of ∼52, and it is about a factor of 100 for the B models. This is due to the fact that we used the same mass for all the particles to avoid unwanted numerical effects arising from the interaction between particles with different masses.

Fig. 1. Onset of each simulated scenario: head-on collision (model A00, left panel), grazing collision (A05, middle), and merger (AM, right). We show here the particles in a thin slice along the orbital plane (X − Y plane) for the A models. The B models have the same setup, correspondingly. The axes are in units of 1010 cm and density is color-coded.

2.2. Initial conditions

The initial configuration of the MS star, which mimics the current Sun, was achieved by randomly distributing SPH particles in 3D, following 1D density, temperature, and composition profiles of the Sun at age ∼4.5 Ga from Bahcall (2005). According to this solar model, the bottom of the convective envelope is located at r ∼ 0.71 R_⊙ and contains a mass of ∼0.024 M_⊙. After relaxing the star, we obtained a configuration that was very close to the initial 1D profiles. We used a stellar Helmholtz equation of state (EOS), which includes the contributions of radiation in the form of Planckian photons, an ideal gas of ions with Coulomb corrections, and an electron-positron gas with an arbitrary degree of relativistic regimes and degeneracy (Timmes & Swesty 2000). This EOS has a range of applicability wide enough to allow us to simulate both the MS and the BD stars, and their interaction. It is important to stress that while the MS is close to an ideal gas, the BD is mostly supported by nonrelativistic degenerate electrons. Having the same EOS for both interacting bodies renders unnecessary any complex treatment of interacting or mixing matter coming from different stars. This would not be the case if the secondary object were a planet.

Finally, we adopted a lithium surface abundance (i.e., in the convective envelope) similar to that observed on the surface of the current Sun (A(Li) = 0.96), while in the interior of the MS star, the Li abundance profile was that for the current age of the Sun obtained in a full evolutionary calculation of a 1 M_⊙ and solar composition stellar model. The computed profile shows that the Li abundance becomes negligible immediately below the convective envelope. We note that the depth and mass of the convective envelope do not change significantly from the zero-age main sequence to the current Sun. The relative evolution of the depth of the convective envelope and the total solar radius is less than 0.2% per 10⁹ a, while the convective envelope mass would decrease from the initial 0.0329 M_⊙ to the current value of 0.02415 M_⊙ (Bahcall et al. 2001). Therefore, if the engulfment occurred during the very early MS evolution, the outcome would not be very different.

To obtain the initial model of the BD, we integrated the mechanical equilibrium configuration using the Helmholtz EOS. At each integration step, we calculated the temperature profile assuming that the degeneracy parameter (η) was constant, fixed with the central density (ρ_c), temperature (T_c), and composition (via the mean molecular weight per electron, μ_e),

$$\begin{array}{c}\hfill {\eta }_{\mathrm{c}}=\frac{{(3{\pi }^2{ħ}^3)}^{2/3}}{2m_{\mathrm{e}}{k}_{\mathrm{B}}{T}_{\mathrm{c}}}{\left(\frac{N_{\mathrm{a}}{\rho }_{\mathrm{c}}}{{\mu }_{\mathrm{e}}}\right)}^{2/3}\sim \frac{3.018\times {10}^5}{T_{\mathrm{c}}}{\left(\frac{{\rho }_{\mathrm{c}}}{{\mu }_{\mathrm{e}}}\right)}^{2/3},\end{array}$$(1)

where the c subscripts mean central values and all other constants have their usual meaning. Hence, the temperature radial profile is obtained as

$$\begin{array}{c}\hfill T(r)=\frac{\eta (r)}{{\eta }_{\mathrm{c}}}{T}_{\mathrm{c}}.\end{array}$$(2)

This approximation is valid for most of the BD structure for a given mass, composition, and elapsed time, as the degeneracy parameter is approximately constant along adiabats (Burrows & Liebert 1993). This method produces reasonable density profiles for BDs given the use of a detailed EOS, such as the Helmholtz EOS used in this work. Although the resulting temperature profile is only qualitatively correct, especially close to the surface of the BD, the details of the initial temperature profile are not very relevant, owing to the strong heating caused by the collision.

In a second step, we mapped these profiles in a 3D random distribution of SPH particles and relaxed the system, letting it achieve a stable configuration. With the BD being a long-lived, cold, and fully convective star, we considered a constant Li abundance throughout, and equal to the accepted current value in the interstellar medium (ISM), A(Li) ∼ 3.0. We used 100 000 (A models) or 52 400 particles (B models) for the BD so that all SPH particles in the simulation had the same mass. Figure 2 shows the initial density profile for the corresponding objects. All SPH particles are shown as dots jointly with the corresponding 1D theoretical density profiles as solid lines. The low dispersion of the SPH profiles and their similarity to the 1D profiles show that the 3D initial conditions adequately represent the theoretical profiles.

Fig. 2. Density profiles for the MS (left panel) and BD (right panel) stars after relaxation. All SPH particles are represented with purple dots. Solid green lines show the initial 1D theoretical density profiles.

Once we have the initial configurations of the individual objects, we can arrange the onset of the interaction. Many studies have addressed the physics of star-planet interactions, some of them focusing on how a planetary orbit changes under the action of tides between the star and the planet (see e.g., Livio & Soker 1984; Siess & Livio 1999a,b; Villaver & Livio 2007; Nordhaus et al. 2010; Bear & Soker 2011, Villaver et al. 2014; Privitera et al. 2016a,b; Stephan et al. 2020). These studies have provided many interesting clues about the initial conditions for an engulfment to occur in terms of the mass of the star, the mass of the planet, the initial distance between the planet and the star, the rotation, and the importance of other physical ingredients such as the mass-loss rates or convective overshooting. However, the details of the previous orbiting evolution of the system and the ultimate reason that may have caused the collision is beyond the scope of this study. Therefore, depending on the simulated scenario, we set both objects apart, with an initial distance d₀. Then, we assigned an initial velocity v_0, rel to both, which corresponded to the free-fall velocity (v_ff) at d₀ for the collisions, and 80% of the orbital velocity (v_orb) for the merger. Figure 1 shows the configurations with respect to the center of the MS for simplicity. Nevertheless, the simulations were performed with respect to the center of mass frame of reference, in which both bodies had nonzero initial velocities.

3. Results

In this section, we present the results obtained in our simulations. Using the A cases (i.e., M_BD = 0.019 M_⊙) as nominal scenarios, we first focus on the dynamics of the system, exploring in detail the (more extreme) cases of head-on and grazing collisions, and finally the smooth merger. Secondly, we discuss the mixing and distribution of BD material after the engulfment, particularly how its Li content is distributed in the host star⁵. Thirdly, we comment on the angular momentum transfer in the case of a grazing collision and the merger scenario. Finally, we discuss the effects of reducing the mass of the BD (the B cases, M_BD = 0.01 M_⊙).

3.1. Analysis of the dynamics

3.1.1. Head-on collision

Model A00 represents our nominal case, corresponding to a head-on collision with zero impact parameter. In the following, we describe the evolution of the collision as it resulted from our 3D simulation.

Figure 3 shows several snapshots taken at different times during the collision (model A00). Figure 4 (left) shows the velocity evolution of the denser region of the BD in comparison with the local speed of sound of the surrounding MS star material. Shortly, after first contact, a bow shock is formed as the BD penetrates the MS star supersonically (with a peak Mach number of ℳ ≃ 11) and it speeds up toward the core of the MS star. After a short ballistic phase, the BD losses part of its mass through a combination of direct mass stripping and ablation. Nevertheless, it manages to cross the core of the MS star, while it heats up and compresses. Due to the extreme friction experienced by the BD, its kinetic velocity is leveraged to heat the material. This slows down the BD, which is not able to break out through the MS star antipodes, making the BD move subsonically until its direction is inverted (t ≃ 2700 s). At this point, the fragmentation of the BD increases the ratio surface-to-volume, which makes it prone to dissolve. Mixing of the BD material with that of the MS star is facilitated by the Kelvin–Helmholtz (KH) instability, which has a characteristic growth time of

$$\begin{array}{c}\hfill {\tau }_{\mathrm{KH}}=\frac{\lambda ({\rho }_1+{\rho }_2)}{\sqrt{{\rho }_1{\rho }_2}\left|{\mathit{v}}_{\mathrm{rel}}\right|},\end{array}$$(3)

Fig. 3. Snapshots of the particle distribution for the head-on collision scenario with η = 0 (model A00) at different times, increasing from left to right and from top to bottom. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and the temperature is color-coded.

with ρ₁ and ρ₂ being the densities of the interacting fluids, λ the size of the growing mode of the KH instability, and v_rel the relative velocity. Assuming that ρ₁ ∼ ρ₂, v_rel ∼ 5 × 10⁷ cm s⁻¹, and λ ∼ R_BD being the mode responsible of the BD destruction, we can estimate τ_KH ∼ 350 s. According to numerical simulations of the destruction of a homogeneous sphere by a supersonic wind, approximately t ∼ 3τ_KH is needed to fully destroy the sphere (see e.g., Agertz et al. 2007; Hopkins 2013; Frontiere et al. 2017; García-Senz et al. 2022). This pushes our estimation for the destruction of the bubble to the order of t ∼ 1000 s after the bow shock is formed (t ∼ 1800 s in Fig. 3), namely, t ∼ 2800 s. In our simulation, neither the density of the BD nor the density of the MS star is constant, and nor is the relative velocity of the fluid while the BD traverses the MS star, which results in longer survival of the BD t ∼ 3500 s, but still similar to our estimation. At this point, we cannot talk anymore of a BD per se, but of blobs of fluid that previously belonged to the BD. Thanks to the Lagrangian nature of SPH, we can follow the fluid elements of the MS star and of the BD independently, even in the case of extreme mixing. In that respect, the remaining fragments of the BD are not in dynamical equilibrium, and they experience again a gravitational pull toward the center of the MS star. This time they completely dissolve because of the combined effect of hydrodynamic instabilities and heat losses through their surface, effectively mixing with the MS star. All these mixing episodes take place in relatively high-density regions of the MS star, keeping the fresh material of the BD in deeper layers and bound to sink to the center of the star. We can see two repetitions of this process (see the first two peaks in the BD velocity curve of Fig. 4) until the BD material is completely mixed with the MS star. A shockwave is launched during the first crossing episode of the BD, which eventually breaks through the outer layers of the MS star, inducing some amount of mass ejection (see Table 2). To calculate a lower limit of the mass lost by the system, we added the mass of all SPH particles whose radial velocity was positive and greater or equal to the escape velocity evaluated at their radius, considering the total enclosed mass and correcting by the position and velocity of the center of mass of the MS star. About 1.9 × 10⁻³ M_⊙ are ejected from the system. Of those, only 7.8 × 10⁻⁶ M_⊙ were originally from the BD. This amount of ejected mass is quite similar to the total mass in the Solar System, excluding the Sun.

Fig. 4. BD velocity evolution (green line) compared with the local speed of sound of the surrounding MS star’s material (purple line). The BD Mach number evolution is shown in blue. The top panels show the A models, while the bottom panels are for the B models.

3.1.2. Grazing collision

We also performed a variant of the previous scenario, this time with a nonzero impact parameter (model A05 with η = 0.5 R_⊙). This scenario is more likely than a head-on collision of the BD, and it may modify the angular momentum of the MS star, which could be detected observationally (see Sect. 3.3 below).

Figure 5 shows snapshots of the distribution of the particles at different times, with the temperature being color-coded. The initial part of the interaction is similar to the A00 case. The BD impacts supersonically (peak ℳ ≃ 12) against the outer layers of the MS star, forming a bow shock that crosses the entire host star. The BD is compressed as it gets closer to the denser core of the host star. Nevertheless, due to the nonzero impact parameter, the BD misses the center of the MS star and spirals around it. In this case, it takes a longer time than case A00 for the BD to be able to slow down. As a consequence, a sizable fraction of the initial spiral trajectory of the BD happens supersonically (see the top central panel in Fig. 4). This launches a spiral shockwave during the first part of the interaction. When the BD reaches more diluted regions of the MS star, it expands and triggers the growth of hydrodynamical instabilities during the next part of the orbit, which mix part of its material (up to 70%) with that of the MS star. The remaining material of the BD reaches the inner core of the MS star half of an orbit later and slows down. At the end of the simulation, about 27% of the BD material is in the inner core of the MS, most of the BD material (45%) is set in a dense toroidal structure around the MS core, and the remaining 28% is in a tight spiral arm. From here on, both structures are slowly pulled in by the gravitational force of the MS core. As in the A00 case, most of the interaction happens within relatively deep layers of the MS star, and hence most of the mass of the BD will end up in the densest regions of the host star.

Fig. 5. Snapshots of the particle distribution for the grazing collision scenario with impact parameter η = 0.5 (model A05) at different times, increasing from left to right and from top to bottom. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and the color represents the temperature.

The outer layers of the host star expand as they receive the impact of the spiral shockwave and a slightly higher amount of material is ejected from the system than in the case of A00 (∼2.2 × 10⁻³ M_⊙). In this case, the amount of BD mass that is lost is almost one order of magnitude smaller than in the A00 model (∼9.5 × 10⁻⁷ M_⊙) for reasons that are explained in Sect. 3.2. It is worth noting the low-temperature material that can be seen at the center of the MS core in the last panel of Fig. 5. This is just a transitory situation. Once the BD loses kinetic energy, its material will thermalize with that of the MS core, but on a longer timescale than that simulated in this work. Therefore, despite its low temperature, this material will be, on all counts, burned anyway.

3.1.3. Merger

As mentioned in Sect. 2.1, in order to induce the merger of the BD, we set it in an orbit (d_orb = 1.5 R_⊙) with an orbital velocity artificially reduced (20% less than the velocity needed for a stable orbit at that distance). This led the BD into a trajectory that brought it close to the external layers of the MS star in a quarter of an orbit. Figure 6 shows the evolution of the merger at different stages with color-coded temperatures. The first interaction with the material of the MS star happens supersonically (ℳ ≃ 8) and due to the relatively low density in the outer layers of the MS star, the BD is effectively sped up by the gravitational potential of the host star during the next half orbit. Part of the BD material (∼10%) is stripped out via ablation and it mixes with the outer layers of the MS star, but the stronger effect is due to tidal forces. By the time the BD reaches the antipodes of the entry point, it has been noticeably distorted and stretched (third panel of the top row in Fig. 6). Due to the combined work of the tidal forces and the brief slingshot acceleration, the BD is stripped into a long spiral arm that is partially accreted on top of the MS core, while the remaining material orbits around the MS star (central row of Fig. 6).

Fig. 6. Snapshots of the particle distribution for the merger scenario (model AM) at different times. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and color represents temperature.

As the interaction proceeds, the spiral tightens and keeps accreting on the outer layers of the MS star (bottom row of Fig. 6) until about 60% of the BD mass is in the innermost regions of the MS star and the remaining 40% is in the spiral arm. The final state of this simulation is still not stationary (last snapshot in Fig. 6), but the timescale needed in order to allow the spiral arm to settle down on the MS star surface is much longer than what can be affordable in a simulation limited by the dynamical timescale.

Nevertheless, the amount of the total mass lost by the system is smaller than in the more violent head-on and grazing collisions: ∼5 × 10⁻⁴ M_⊙ (see Table 2), and no material belonging to the BD escapes the system, meaning that the full spiral arm will eventually fall back. Figure 7 (right panel) shows that this material would form a circumstellar disk in the plane of the collision, with possible observational consequences in terms of infrared excess, provided that this structure lasts for a long enough time.

Fig. 7. Particle distribution of the BD projected in the 2D cylindrical plane for the A scenarios, at the end of the simulations: t = 10 689 s and t = 8059 s in the head-on and grazing collision cases (left and central, respectively) and t = 38 402 s for the merger scenario (right). The semicircles represent the original solar radius (solid black line) and the bottom of the convective zone at ∼70% of the solar radius (solid purple line). The density is color-coded.

3.2. Li abundance

In our scenarios, the BD fluid elements act as tracers of Li-rich material. Because of the Lagrangian nature of SPH, we can directly track where the BD material is located throughout the whole simulation, simply by following the SPH particles of the BD. Determining where this material ends within the MS star’s structure allows us to predict the expected observable Li abundance after the collision or merger. Admittedly, we can only simulate several thousands of seconds with a full 3D hydrodynamics code, as the collision is very dynamic and we are limited by the Courant condition. To give a final answer on where the Li-rich material ends, it would be necessary to continue this calculation up to millions of years with a 3D stellar evolution code, which nowadays is impossible. Due to the highly axial asymmetry of the final outcome of the simulations, the common procedure of mapping the final 3D structure into a shell-averaged 1D profile and following the evolution with a spherically symmetric hydrostatic code is not a trivial task. Nevertheless, we can provide some constraints on the fate of the BD material and guidelines on what to expect from scenarios such as those simulated here.

We start by discussing how the MS and BD materials are spatially distributed at the end of the simulations. Figure 7 shows a cylindrical projection of all SPH particles of the BD at the end of each simulation for the A cases. The first noticeable detail is that the BD material is not homogeneously distributed after the collision. There is a higher concentration along the equatorial plane (namely, the collision plane) than along the polar plane. This density inhomogeneity may have an impact on the subsequent evolution of the MS star. However, this issue must be addressed with a 3D stellar evolutionary code, which is beyond the scope of the present study. As mentioned above, the final outcome of our simulations has a fully inhomogeneous non-axisymmetric structure in density, temperature, and composition. Therefore, the extraction of any averaged 1D structure carries the risk of losing the information regarding the collision with the MS star. For the cases of head-on and grazing collisions, most of the BD material is confined in the deep interior of the MS star, while only a small fraction remains above the original convective zone and is still gravitationally bounded (1.32% and 0.35% for the A00 and A05 models, respectively), which we assume remains unperturbed at r ∼ 0.7 R_⊙, as well as the original stellar mass above this point (∼0.02 M_⊙, see also Table 2) after the collision has settled down to its final configuration. As we cannot follow the evolution of the final object until it reaches equilibrium⁶, the location of the new convective envelope might be different. This rough factor of two difference in the surface pollution between both scenarios can be understood when looking in detail at the dynamics of each collision. The A00 case has a clear preferred direction of movement for the ablated material of the BD: the impact direction. As the BD supersonically carves its path through the MS star’s material, it leaves behind a temporary low-pressure funnel that allows the BD material to escape toward the MS star’s surface. In the case of model A05, because the BD spirals in, the low-pressure region behind the BD follows the same in-spiral path, away from the impact line, diluting any preferred direction for the BD matter to escape. This effect can be seen in Fig. 8, where we present a 2D cut along the collision plane for both collision scenarios of case A. The color of the particles represents their radial velocity. It is clear that a preferred direction with higher outward radial velocities exists in the case of model A00, while the material moving outward in the case of model A05 has to push against layers of the inert MS star’s material.

Fig. 8. 2D cut of the 3D particle distribution for the A00 (left) and A05 (right) cases at t = 2600 s. The radial velocity is color-coded.

The case of the merger is completely different. As the BD spends a relevant amount of time in the diluted regions of the MS star, the ablated BD material remains in these outer layers and it never reaches denser and hotter regions. It is only when the BD moves closer to the MS core that tidal forces are strong enough to rip the BD apart. Once this happens, the BD material is pulled in through a thin stream of matter, which, due to angular momentum conservation, spirals around the MS core. The flow of BD material is not able to reach the central core of the MS, as it does in the collision cases. Nevertheless, it is able to deposit Li-rich material on the surface of the inner core of the MS, where it will be destroyed, nonetheless. A sizable fraction (∼40%) of the BD material is left on the remainder of the spiral arm, whose evolution happens on a longer timescale and it could eventually be allocated on the outer layers of the MS, within the convective region and its surface.

Figure 9 shows the time evolution of the percentage of the BD mass that has radius r ≥ 0.7 R_⊙ and is still gravitationally bound, for all simulated models. Solid lines correspond to the A models. Here we can see that we reached a steady state on all simulations regarding the amount of BD mass that could eventually remain above the original MS convective region. The trend follows the results discussed above: head-on and grazing collisions have a lower capability to pollute the external layers of the host star than a merger.

Fig. 9. Time evolution of the percentage of gravitationally bound BD mass with r ≥ 0.7 R⊙. The lower x-axis is for collisions, while the upper x-axis is for mergers. The solid lines are for the A models and the dashed lines are for the B models.

It is important to note that in the merger case, about 50% of the BD mass that has r > 0.7 R_⊙ and is still gravitationally bound, shows an inversion of the entropy gradient. This inversion is the result of a mixture of ablated material that has been strongly heated by the shockwave during the first phases of the interaction and material expelled by simple angular momentum transfer during the destruction of the BD by tidal forces. Figure 10 characterizes the significant heating that the BD material suffers, owing to the relatively high density of the environment, and the influence on its entropy profile. The panel on the left shows the initial profiles of temperature and entropy (in k_B per baryon) for both stars. The panel on the right shows the final profiles for the same magnitudes. We note the change in temperature of the BD material (about two orders of magnitude closer to the MS star center), and the overall entropy increase in the outer layers of the MS star (∼ × 2) and the whole BD material (∼ × 2 − 4). In the latter, the entropy inversion is clearly seen for r ≳ 0.7 R_⊙, denoted by the solid vertical line. As a result, this entropy profile points to a convective unstable layer. The long-term effect of this instability is unclear, as the end state of all simulations is still too dynamic to predict the actual final location and state of those BD particles. Nevertheless, some mixing is expected to happen on a timescale longer than can be simulated with a 3D code nowadays.

Fig. 10. Initial (left) and final (right) temperature and entropy radial profiles for both stars in the AM case. We show all SPH particles in all curves. The profiles in the left panel are centered on the center of mass of each star, while in the right panel they are centered on the center of mass of the remnant. The vertical black line in the right panel shows the bottom of the convective envelope.

Figure 11 shows the final coarse particle distribution of our merger simulations at a large scale. Matter from the MS star is represented by purple points, while BD matter is in green. Blue points are the particles that will escape the system. As no matter from the BD reaches escape velocity (see ejected BD mass in Table 2), it is clear from this plot that the final distribution of Li-rich material would be mostly located around the collision plane in an inhomogeneous way, forming a relatively thin spiral arm. Additionally, the outer layers of the MS have expanded considerably, although most of this matter is still gravitationally bound and will fall back. This figure gives a qualitative idea of how the Li-rich material is distributed at the end of the simulation, but not of its quantity, which was discussed in Sect. 3.1.3.

Fig. 11. Spatial distribution for the MS star (purple) and the BD (green) material at the end of the AM (left column) and BM (right column) merger scenarios. Particles in blue have radial velocity larger than the local escape velocity and all of them belong to the MS star. Every row shows a thin slice in a different plane of the same 3D distribution. All panels have a size of 2 × 1012 cm per side and the black circle in the center of each panel shows the original radius of the MS star.

Besides the spatial distribution of the BD material, it is important to track their thermal evolution during the simulated time to assess if this matter ever reaches Li-burning temperatures (≳2.5 × 10⁶ K) and if it does, for how long. To this extent, we tracked the maximum temperature of each BD SPH particle and the amount of time that they spent at temperatures larger than or equal to the previous temperature threshold. We show the results in Fig. 12. The points show the maximum temperature (T_7, max, in 10⁷ K units) reached by each particle of the BD throughout the whole simulation in function of its radial coordinate at the end of the calculation. The color represents the amount of time that a specific particle spent with temperatures higher than or equal to the Li burning temperature, shown by the solid line. In Fig. 13 we show the characteristic burning time of Li for the two main channels of proton captures in function of temperature. The color of the lines represents the reaction channel and the type of line is for different density values that cover most of the range in our simulations.

Fig. 12. Maximum overall achieved temperature of the BD material in function of the radial distribution at the end of all simulated models. The color represents the amount of time for which a specific particle had temperatures higher than 2.5 × 106 K (i.e., the assumed Li burning temperature limit, shown with a solid green line).

Fig. 13. Characteristic timescales for the two main channels of Lithium burning (7Li(p,α)4He in purple and 7Li(p,γ)7Be in green) for different densities and in function of the temperature in 107 K units (T7).

In the collision cases, almost the full mass of the BD reaches temperatures much higher than T₇ = 0.25, but a fraction of those (1.4% of the BD mass for the A00 and B00 cases, and 0.5% for the A05 and B05 cases) spends less than 1 h with such temperatures and cools down to below the Li-burning temperature limit within 1000 − 3000 s. In particular, these include all those particles that end up at r ≥ 0.7 R_⊙. The 05 cases (middle panels of Fig. 12) take a longer time to interact with the denser layers of the MS star and, consequently, the BD material reaches lower temperatures than the 00 cases, and they spend less time above the Li-burning limit. This points to the fact that, despite reaching very high temperatures, a fraction of the BD material is able to cool down relatively fast, and remains unburnt in the outer layers of the MS star. As mentioned in Sect. 3.1.2, the decrease in temperature for the BD material at a low radius is transitory, and this material is considered to be burnt on a longer timescale. Cases A05 and B05 show very similar behavior.

On the other hand, in the merger scenarios, almost 40% of the BD mass never reaches the Li-burning temperature limit and an additional 4% spends less than 1 h at such temperatures. As it can be seen in Fig. 12 (right panel), most of this material remains at radius ≥0.7 R_⊙. The reason for this is due to the longer time that the BD spends interacting with the outer layers of the MS star. This interaction is enough to peel off material of the BD and develop a long thin spiral arm, which neither reaches the deeper layers of the MS star nor is substantially heated. The BD material in the merger simulations reaches minimum radii that are approximately one order of magnitude larger than in the collision scenarios. This suggests that merger scenarios are more prone to result in larger lithium contamination of the outer layers of the MS star.

On the other hand, all BD material with r ≤ 0.7 R_⊙ experiences temperatures much higher than T₇ = 0.25 for timescales comparable to the total simulated time in all simulations. This points to the fact that their lithium will most likely be completely burned on a relatively short timescale, as expected for the deeper regions of the host star. This limit between burnt and unburnt material, which can be clearly seen in Fig. 12 around r ∼ 0.7 R_⊙, independently of the simulated scenario, appears naturally. This led us to think that, despite the grand display of material being ejected from the system and the BD being destroyed, the mass change is very low in comparison with the total mass of the MS star. Therefore, it is very likely that the location of the bottom of the convective envelope does not significantly change, and hence our assumption of considering it still at r ∼ 0.7 R_⊙ throughout all the simulations.

Considering the above-quoted percentages of the BD material that remains within the former convective region, it is straightforward to estimate the new Li surface abundance in the MS star after the collision and merger scenarios (see e.g., Jackson & Carlberg 2018, their Eq. (10)). As mentioned above, a sizable fraction of this material is marginally exposed to temperatures higher than T₇ = 0.25 during the simulations (Fig. 12). For the typical density and temperature values existing at the bottom of the convective envelope in the current Sun (∼0.75 g cm⁻³ and T₇ ∼ 0.2, respectively), the characteristic time for Li burning via the ⁷Li(p,α)⁴He reaction is ∼2 × 10¹¹ a. Then, assuming the lithium abundance currently observed in the Sun as the surface abundance before the collision, the new surface abundance values would be A(Li) = 1.33 and 1.09 for the A00 and A05 models, respectively. These increments in the absolute abundance of Li are small and it would be difficult to differentiate from the typical Li abundance observed in solar-like stars. However, in the case of the merger, the Li enhancement is considerably larger due to the significant amount of the BD mass deposited in the envelope of the MS star: in the AM model, the new Li abundance would increase by a factor larger than twenty (A(Li) ∼ 2.46). This huge increase in the Li abundance would be easily detectable observationally, provided that this surface Li enhancement remained for a long enough time. It should be noted that the trigger of the thermohaline convection (see Sect. 1), which would eventually erase this Li enhancement, would not develop in the scenarios studied here, as the surface metallicity is not altered due to the engulfment; the MS star and the BD share identical metallicity, and only differ in their Li content.

3.3. Angular momentum transfer

Except for the case of a head-on collision, all interactions should transfer angular momentum to the MS star. This is an important observable, which might hint at a history of engulfment. For instance, 1D simulations performed by Carlberg et al. (2009), Privitera et al. (2016a,b) showed that the interaction between a red giant and a planet followed by planet engulfment can indeed accelerate the surface rotation for a sufficiently long time to produce fast-rotating red giants with rotational surface velocities above ∼8 km s⁻¹. Rapid rotator red giants exist and usually show some Li enhancement compared to the slow rotator red giants, although the origin of this Li enhancement is still unclear (see e.g., Carlberg et al. 2013; Mallick et al. 2022). Also, recent 1D simulations for MS stars by Oetjens et al. (2020) suggest that the dynamical interaction between an MS star (in the range 0.8 − 1.0 M_⊙) and its companion (1 − 20 M_J) can induce an increase in the rotational velocity of the host star. This occurs while it is still on the MS, from 1 to over 40 km s⁻¹ if an engulfment occurs, depending on the initial semimajor axis, the planet-to-star mass ratio, and the metallicity of the host star.

We did not include any initial intrinsic rotation of the colliding bodies in our simulations. A rotation such as that of the present Sun (∼2 km s⁻¹) is extremely small in comparison with the highly dynamic interaction of the stars. For the longest simulations (mergers), which are on the order of several 10⁵ s, the outer layers of the MS star would have rotated less than 10⁶ km. This is three orders of magnitude smaller than the spatial resolution in the lower-density regions of the MS star. The rotational profile of a BD is less known and high rotational velocities may be expected (v sin i > 10 km s⁻¹, see e.g., Konopacky et al. 2012). However, these values are still too small to produce significant changes within dynamical timescales. Nevertheless, we note that three extremely fast rotating BDs have been recently observed (Tannock et al. 2021); one of them with the highest v sin i ∼ 100 km s⁻¹ value ever reported for an ultra-cool BD.

We measured the evolution of the angular momentum of the fluid within 1 R_⊙ from the corresponding center of mass for all simulated models. From that, and calculating the moment of inertia, we could derive an average angular velocity of the fluid at all times. We restricted the calculation to a sphere with radius r = 1 R_⊙, which will provide a lower boundary to the final angular velocity, as we expect this to increase when the expelled (but still gravitationally bounded) matter will fall back onto the host star. We also conducted the same analysis with a larger radius (r = 1.2 R_⊙), which yielded similar results. Figure 14 shows the evolution of the angular velocity for all simulations and Table 2 shows the corresponding final period of the MS star. As expected, the head-on cases produce almost no change in the rotation, while grazing collisions and mergers have a remarkable effect. The first peak and subsequent dip in the evolution of the mergers for t = 5000 − 12 500 s is a direct consequence of how we calculated the angular momentum. Because we restricted it to r = 1 R_⊙, the value changes drastically when the BD crosses that limit, in and out, during its first orbit. After that, the trajectory of the accreted BD material always remains within that radius and the evolution of the angular velocity is smoother (increasing slowly as material keeps being accreted). Model A05 increases the rotational velocity of the MS star to 46 km s⁻¹, while model AM increases it to ∼52 km s⁻¹ at the end of the simulation. For a given choice of the initial relative velocity and impact parameter (Table 1), this increment factor would be similar for a similar mass ratio between the BD and the MS star, in agreement with the 1D simulations by Oetjens et al. (2020; see Sect. 3.4). Regarding the radial profile of the angular velocity for the final objects of the A05 and AM models, we found that they are in a transitory situation. The core of the MS star (∼0.2 M_⊙) exhibits a strong differential rotation, while the remainder of the star shows a rotation closer to rigid body, increasing only in the outer layers, where mass is still accreting. The fact that the average angular velocity shown in Fig. 14 is stable for most of the last part of all the simulations, despite its radial variability, suggests a constant redistribution of angular momentum so that the provided angular velocity is a possible observable. Additionally, the fast-rotating central core might have observational implications for asteroseismology studies through the detection of g-modes in stars that are suspected to have recently suffered a sub-stellar engulfment (see e.g., Fossat et al. 2017; Eggenberger et al. 2019).

Fig. 14. Average angular velocity evolution for the A (solid lines) and B (dashed lines) models. The calculation is restricted to the fluid within R = 1 R⊙.

Current observational measurements of the rotational velocity in MS field-star samples are still too small to fully test our model. Nevertheless, the values obtained can be compared with the equatorial angular velocity of the Sun (∼2 km s⁻¹) and with the v sin i values measured in a few MS stars of similar spectral type that are suspected to have suffered a planet engulfment and/or the accretion of rocky material. For instance, it has been found that a significant fraction among a sample of F stars hosting a HJ show average v sin i values larger than those stars in their sample (in the same effective temperature region) without a HJ (Delgado Mena et al. 2015). These moderate rotators show, in addition, lower Li abundances (by 0.14 dex), suggesting that rotational mixing (possibly induced by planet interaction and/or engulfment) might be the cause of a greater depletion of Li. On the contrary, in the 16 Cyg pair, star A (with a Li excess supposedly due to the accretion of some rocky material) shows a slightly larger rotational velocity (2.23 km s⁻¹) than star B (1.27 km s⁻¹; Davies et al. 2015). Similar values are observed in the co-moving pair Kronos and Krios, and the star HD 240430, with Li abundance ∼ + 0.5 dex larger, shows also a larger rotational velocity that HD 240429 (2.5 vs. 1.1 km s⁻¹, Oh et al. 2018).

All these observed rotational velocity values are compatible with those expected for MS stars of their masses and age (Skumanich 1972; Tassoul 2000), and contrast with the increase of more than one order of magnitude in the angular velocity expected from interactions such as those presented here, in the A models. The observations that exist so far do not show MS field stars with masses and ages similar to the Sun, with v sin i ≳ 10 km s⁻¹ (Donahue et al. 1996; Stauffer 2004). Such large rotational velocities are only observed in some solar-like MS stars belonging to very young open clusters (e.g., in the Pleiades, Queloz et al. 1998; Terndrup et al. 2000; Rebull et al. 2016; Praesepe, Douglas et al. 2019; or NGC 6811, Curtis et al. 2019). We note, however, that for Sun-like cool stars, the torque triggered by the magnetic braking mechanism increases with the increasing angular velocity, implying, as a consequence, a spin-down rate $\dot{\mathrm{\Omega }}\sim ({\mathrm{\Omega }}_o/2t_o){({\mathrm{\Omega }}_{\star }/{\mathrm{\Omega }}_o)}^3$; with t_o = 2 × 10⁸ a and Ω_o = 10 km s⁻¹/R_⋆ (see e.g., Skumanich 1972; Soderblom 1983; Pantolmos & Matt 2017). Therefore, for R_⋆ ∼ 1 R_⊙ and the surface Ω_⋆ values quoted above, the rotational velocity would be reduced in a time considerably lower than ∼10⁹ a. This would make the observational detection of any rotational velocity excess induced by an engulfment, such as those studied here, difficult, considering the relatively low frequency of these interactions. Nevertheless, there is observational evidence for a possible dependence of the braking law on the stellar mass, resulting in K dwarfs spinning down much more slowly than F and G dwarfs, and that this braking may stall for an extended period of time, the duration of which increases toward lower masses and cooler temperatures (e.g., Curtis et al. 2019). In any case, one should wait for more extensive and dedicated surveys on the rotational velocity in solar-like field stars before a definite conclusion is reached.

3.4. Changing the mass of the BD

To explore the impact of the mass of the engulfed star, we reduced the mass of the BD adopted in the A models (0.019 M_⊙) to 0.01 M_⊙ in the B models (i.e., closer to a giant planet mass), and we recalculated the same three interactions with the same MS stellar mass and initial configurations. The relaxation process for this lighter BD was the same as for the A models, and it used the same EOS, assuming the same Li abundance. In the right panel of Fig. 2, we show the initial density profile of the 3D particle distribution after relaxation, compared with the original 1D theoretical profile. Being a factor of 2.8 less dense, this lighter BD is expected to be destroyed more easily than the heavier one.

In the case of the head-on and grazing models (B00 and B05), the dynamics are similar to those already discussed regarding models A00 and A05. The main difference is that the BD is destroyed faster by its interaction with the MS star, as expected. The B collision models also eject three to four times less mass than the A collision models (see Table 2), owing to the lower kinetic energy stored in the less massive BD at the time of the impact. Interestingly, model B00 is able to eject three times more mass from the BD than model A00, most likely due to its lower compactness, which makes the ablation process more efficient. The case of the merger (model BM) also presents a similar evolution to that of model AM (see Fig. 15). The BD is again torn apart by tidal forces in half of an orbit, but the resulting spiral arm is more diluted in this case and it mixes rapidly with the surrounding material of the MS star. Overall, the merger produces a similar amount of ejected mass as in case A, belonging exclusively to the MS in both cases.

Fig. 15. Snapshots of the particle distribution for the merger scenario with the lower-mass BD (model BM) at different times. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and the temperature is color-coded.

The BD mass distributions at the end of the simulations (see Fig. 16) follow the same pattern as in the A models. In fact, the fraction of the BD mass that remains in the outer layers of the MS star is very similar to that in the A models (see Table 2). This can also be seen in the final periods of the MS star, which follow the same hierarchy as in the A models, albeit with slightly lower angular velocities due to the lower total angular momentum of the system with the smaller BD.

Fig. 16. Particle distribution of the BD projected in the 2D cylindrical plane for the B scenarios, at the end of the simulations: t = 11 478 s and t = 11 923 s in the collision cases (left and central, respectively), and t = 29 821 s for the merger scenario (right). The semicircles represent the original solar radius (solid black line) and the bottom of the convective zone at ∼70% of the solar radius (solid purple line). The density is color-coded.

Figure 9 shows the time evolution of the percentage of the BD mass with r ≥ 0.7 R_⊙ for the B cases with dashed lines. In these cases, the results show that the material of the BD can be ablated more easily when the BD mass is lower, as expected because of its lower gravitational energy content. Figure 16 shows the cylindrical projection of all BD particles at the end of each simulation. Comparing it with Fig. 7, the same conclusions as for the A models apply here. The main difference appears in model BM with respect to model AM, as the BD mass in the outer regions is distributed in a thinner disk than in case AM, more concentrated around the collision plane but still with a similar density at the equator. Therefore, similarly to model AM, in model BM, a circumstellar disk of matter around the MS star might also form, which eventually may produce an infrared excess. We suggest performing dedicated studies for the detection of infrared excess in MS stars with anomalous Li abundances. We note that large infrared excesses have been found in K-giants with enhanced Li abundances, which, in addition, are usually fast rotators (Rebull et al. 2015), although the actual origin of the Li enhancement is still unclear (Zhou et al. 2019).

As in the A models, it is straightforward to estimate the new Li abundance in the envelope of the MS star at the end of the B model simulations. We obtained A(Li) ∼ 1.31, 1.13, and 2.24 for models B00, B05, and BM, respectively. These values are very similar to those obtained in the A models. Therefore, the same conclusions on the observational detectability of this Li increase mentioned in Sect. 3.2 hold in the case of a lower BD mass. Although the fraction of the BD mass deposited in the envelope is larger in the B models, the net mass deposited is lower.

Regarding the angular momentum transfer to the MS star, Fig. 14 shows the time evolution of the angular velocity of the fluid with r ≤ R_⊙. It is not surprising that there is a dependence on the mass of the BD. The factor between the final angular velocity (see Table 2) of models with different mass is ∼2, which is very similar to the mass ratio of the BDs of the A and B models. Finally, the same trend observed in the A models between types of collisions and the final angular velocity is observed in the B models. Given that very little BD mass is lost by the system (see Table 2, third column), a linear relationship between the final angular velocity and mass of the impacting body is expected. Therefore, an impact similar to the impacts simulated here, with a Jupiter-mass planet instead of a BD, should lead to angular velocities about eight times lower than the B cases just after the event takes place, namely v sin i ∼ 3.1 km s⁻¹ for a grazing collision with an impact parameter of 0.5 R_⊙, and v sin i ∼ 3.6 km s⁻¹ for a merger-like scenario. These rotational velocities are of the same order as those typically observed in solar-like stars. Therefore, the rotational velocity would not serve as an observational hint to detect singular collisions and/or engulfments of a planetary mass similar to that of Jupiter or less by a solar-like MS star.

Finally, the fate of the Li-rich material has a very similar trend as in the A models. As Fig. 12 (second row) shows, all BD material with r ≲ 0.7 R_⊙ will be inexorably burnt. Only the material in the outer layers never reaches the Li-burning temperature, or if it does, it is for a short time span.

Looking at Fig. 12, it is interesting to note that a clear boundary appears naturally in all simulations between the burnt and unburnt material at r ∼ 0.7 R_⊙, which is traditionally assumed as the location of the bottom of the convective envelope. No previous assumption was made in this respect in our code, yet the boundary appeared naturally, independently of the mass of the BD. Despite the highly dynamic interaction and the fact that the final body was still not in equilibrium at the end of the simulations, this result reinforces our premise that the location of the bottom of the convective envelope will not change considerably.

4. Summary and conclusions

In this work, we have explored the effects of the engulfment of a BD by a solar-like MS star. Our main goal was to study if this interaction might significantly alter the lithium abundance on the surface of the host star and induce changes in its main properties, with possible observational consequences. To do so, we conducted 3D hydrodynamical simulations, exploring different types of interaction (head-on and grazing collisions, and mergers) and the impact of changing the BD mass. Our main findings are summarized as follows:

• In all the explored scenarios, the majority of the BD mass remains in the inner parts of the MS star, even reaching the stellar core. This is at odds with previous suggestions, based on 1D simulations, that BD with masses ≤0.02 M_⊙ would mostly be diluted in the envelope of the MS star. In our M_BD = 0.019 M_⊙ case, collisions produced very little pollution in the outer shells of the MS star (∼0.3 − 1.4% of the BD mass). Even with a zero impact parameter, where a low-pressure funnel eases the flow of BD material to external layers, almost all BD material will be confined in the high-density region of the host star. When decreasing the mass of the BD to 0.01 M_⊙, the percentage of its mass that remains in the outer layers increases only slightly (∼0.9 − 2.4%). Thus, on the basis of the present simulations, we cannot elucidate if a complete dilution of the impacting object in the MS star’s envelope would occur at much lower masses, such as that of a giant planet.

• In the case of mergers, a sizable fraction (≃40% for both models) of the BD mass settles in the outer layers of the MS star and never reaches Li-burning temperatures.

• We find a Li enhancement in the envelope of the MS star at the end of the simulated interactions. Such enhancement is very small in the head-on and grazing collision scenarios for the two BD masses studied. Nevertheless, the enhancement would be close to a factor of 20 − 30 in the case of a merger. In this latter case, assuming that the enrichment remains long enough in the envelope, it could be observationally detected as a Li enhancement in surveys of solar-like MS stars. It is worth noting that these values are indicative, as the final state in our simulations is still evolving (albeit on a timescale longer than the dynamical scale). Therefore, additional work should be done to explore the long-term evolution of the system.

• There is a clear increase in the rotational velocity of the MS star after the interaction, being larger in the case of a merger (∼52 km s⁻¹ for AM, and ∼29 km s⁻¹ for BM), although quite similar to that produced by a grazing collision (∼46 km s⁻¹ for A05, and ∼25 km s⁻¹ for B05). This signal is nevertheless transitory. Even if the rotational velocity can increase even more due to the fallback of the MS and BD material that is still gravitationally bound, mechanisms such as magnetic breaking would spin-down the star in short timescales (≲10⁹ a).

• As most of the BD mass remains gravitationally bound to the MS star after the interaction, the angular momentum transfer is directly proportional to the mass of the BD and the dynamical initial conditions. Assuming a similar kinetic content and type of impact, this allows us to extrapolate the rotational velocity of the host star in the case of an interaction with a Jupiter-like mass planet, being in the range of 2 − 5 km s⁻¹, depending on the type of interaction. This would discard rotational velocity as a proxy to select solar-like MS stars that may have suffered an engulfment of a planet-mass object. High rotational velocity values might be used instead as a probe to identify planet and BD interactions with red giant stars, as shown by Privitera et al. (2016a,b). This will be the subject of a similar forthcoming study.

• A significant amount of mass is lost from the system (in the range 10⁻⁴ − 10⁻³ M_⊙) in all our simulations. The overwhelming majority of this mass belongs to the host star. In the case of a merger, this mass may form a disk of circumstellar matter, which may have observational consequences such as an infrared excess. This agrees with previous results from 1D studies on the interaction of a sub-stellar object with an MS (e.g., Stephan et al. 2020) and/or a giant star (e.g., Soker 1998, 2001).

It has to be emphasized that, although the evolving timescale at the end of the simulation is much larger than the dynamical timescale, the final fluid distribution is still not in complete equilibrium. To follow the simulation until all gravitationally bound matter falls back and the final object reaches a stable configuration is out of reach for a full 3D hydrodynamical simulation limited by the Courant condition. To know the final distribution of matter, a stellar evolution code should be used. Nevertheless, we want to stress the asymmetric nature of the simulated interactions and their lack of hydrostatic equilibrium. In order to follow its evolution, we would need a 3D hydrodynamical stellar evolution code, which is not currently available.

In summary, our 3D simulations show that the outcome of the engulfment of a sub-stellar object by an MS star depends on the geometry of the engulfment and the mass of the colliding object. The main differences are the induced rotational velocity and the amount of mass of the sub-stellar object, which is either diluted in the outer layers or in the interior of the MS star. As a result, some of these interactions yield a considerable mass loss and a measurable change in the surface Li abundance, both of which could be observationally detected.

Acknowledgments

This work has been supported by the European Research Council (FP7) under ERC Advanced Grant Agreement No. 321263 – FISH, by the Swiss Platform for Advanced Scientific Computing (PASC) projects DIAPHANE and SPH-EXA, and by the SKACH consortium through funding by SERI. It has also been supported by the Spanish projects PGC2018-095317-B-C21, PID2021-123110NB-100 financed by the MCIN/AEI FEDER “Una manera de hacer Europa”, and by the MINECO Spanish project PID2020-117252GB-100. The authors acknowledge and thank the invaluable support of the scientific computing core facility sciCORE at the University of Basel (https://scicore.unibas.ch/), the Swiss National Supercomputing Center CSCS (https://cscs.ch), the Centro de Servicios Informáticos y Redes de Comunicación CSiRC at the University of Granada (https://csirc.ugr.es/), and the Instituto de Astrofísica de Canarias (https://www.iac.es/). All these institutions provided infrastructures and resources to perform these calculations.

References

All Tables

All Figures

	Fig. 1. Onset of each simulated scenario: head-on collision (model A00, left panel), grazing collision (A05, middle), and merger (AM, right). We show here the particles in a thin slice along the orbital plane (X − Y plane) for the A models. The B models have the same setup, correspondingly. The axes are in units of 1010 cm and density is color-coded.
In the text

	Fig. 2. Density profiles for the MS (left panel) and BD (right panel) stars after relaxation. All SPH particles are represented with purple dots. Solid green lines show the initial 1D theoretical density profiles.
In the text

	Fig. 3. Snapshots of the particle distribution for the head-on collision scenario with η = 0 (model A00) at different times, increasing from left to right and from top to bottom. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and the temperature is color-coded.
In the text

	Fig. 4. BD velocity evolution (green line) compared with the local speed of sound of the surrounding MS star’s material (purple line). The BD Mach number evolution is shown in blue. The top panels show the A models, while the bottom panels are for the B models.
In the text

	Fig. 5. Snapshots of the particle distribution for the grazing collision scenario with impact parameter η = 0.5 (model A05) at different times, increasing from left to right and from top to bottom. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and the color represents the temperature.
In the text

	Fig. 6. Snapshots of the particle distribution for the merger scenario (model AM) at different times. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and color represents temperature.
In the text

Fig. 7. Particle distribution of the BD projected in the 2D cylindrical plane for the A scenarios, at the end of the simulations: t = 10 689 s and t = 8059 s in the head-on and grazing collision cases (left and central, respectively) and t = 38 402 s for the merger scenario (right). The semicircles represent the original solar radius (solid black line) and the bottom of the convective zone at ∼70% of the solar radius (solid purple line). The density is color-coded.

In the text

	Fig. 8. 2D cut of the 3D particle distribution for the A00 (left) and A05 (right) cases at t = 2600 s. The radial velocity is color-coded.
In the text

	Fig. 9. Time evolution of the percentage of gravitationally bound BD mass with r ≥ 0.7 R⊙. The lower x-axis is for collisions, while the upper x-axis is for mergers. The solid lines are for the A models and the dashed lines are for the B models.
In the text

Fig. 10. Initial (left) and final (right) temperature and entropy radial profiles for both stars in the AM case. We show all SPH particles in all curves. The profiles in the left panel are centered on the center of mass of each star, while in the right panel they are centered on the center of mass of the remnant. The vertical black line in the right panel shows the bottom of the convective envelope.

In the text

Fig. 11. Spatial distribution for the MS star (purple) and the BD (green) material at the end of the AM (left column) and BM (right column) merger scenarios. Particles in blue have radial velocity larger than the local escape velocity and all of them belong to the MS star. Every row shows a thin slice in a different plane of the same 3D distribution. All panels have a size of 2 × 1012 cm per side and the black circle in the center of each panel shows the original radius of the MS star.

In the text

	Fig. 12. Maximum overall achieved temperature of the BD material in function of the radial distribution at the end of all simulated models. The color represents the amount of time for which a specific particle had temperatures higher than 2.5 × 106 K (i.e., the assumed Li burning temperature limit, shown with a solid green line).
In the text

	Fig. 13. Characteristic timescales for the two main channels of Lithium burning (7Li(p,α)4He in purple and 7Li(p,γ)7Be in green) for different densities and in function of the temperature in 107 K units (T7).
In the text

	Fig. 14. Average angular velocity evolution for the A (solid lines) and B (dashed lines) models. The calculation is restricted to the fluid within R = 1 R⊙.
In the text

	Fig. 15. Snapshots of the particle distribution for the merger scenario with the lower-mass BD (model BM) at different times. We only show particles within a thin slice along the orbital plane. Each snapshot shows a box with side 2 × 1011 cm and the temperature is color-coded.
In the text

Fig. 16. Particle distribution of the BD projected in the 2D cylindrical plane for the B scenarios, at the end of the simulations: t = 11 478 s and t = 11 923 s in the collision cases (left and central, respectively), and t = 29 821 s for the merger scenario (right). The semicircles represent the original solar radius (solid black line) and the bottom of the convective zone at ∼70% of the solar radius (solid purple line). The density is color-coded.

In the text