First star survivors as metal-rich halo stars that experienced supernova explosions in binary systems

The search for the ﬁrst stars formed from metal-free gas in the universe is one of the key issues in astronomy because it relates to many ﬁelds, such as the formation of stars and galaxies, the evolution of the universe, and the origin of elements. It is not still clear if metal-free ﬁrst stars can be found in the present universe. These ﬁrst stars are thought to exist among extremely metal-poor stars in the halo of our Galaxy. Here we propose a new scenario for the formation of low-mass ﬁrst stars that have survived until today and observational counterparts in our Galaxy. The ﬁrst stars in binary systems, consisting of massive- and low-mass stars, are examined using stellar evolution models, simulations of supernova ejecta colliding with low-mass companions, and comparisons with observed data. These ﬁrst star survivors will be observed as metal-rich halo stars in our Galaxy. We may have identiﬁed a candidate star in the observational database where elemental abundances and kinematic data are available. Our models also account for the existence in the literature of several solar-metallicity stars that have space velocities equivalent to the halo population. The proposed scenario demands a new channel of


Introduction
The first stars in the universe draw much attention for the understanding of the star formation history in the early Galaxy. Such attention involves the critical questions: do they still exist? how can they be found? Answers to these questions have a strong impact on the study of how our Sun was formed and where the terrestrial elements originated. This is because the first stars give a hint towards the formation and evolution of stars, and the synthesis of elements in stars, starting from elemental abundances equivalent to those in the Big Bang nucleosynthesis.
The existence of the first stars in the current universe is still controversial. The basic understanding is that they are typically too massive to exist for a period of more than 10 Myr from their birth. The preference for the formation of massive stars is due to the lack of elements that work as a coolant to compress star-forming gas (Bromm & Larson 2004). However, there are some arguments about the initial mass of the first stars. There is a possibility of the formation of a star with initial mass below the solar mass that has a lifetime comparable to or larger than the age of the universe (Clark et al. 2011;Susa et al. 2014;Stacy et al. 2016). In particular, recent simulation studies of the first stars focus on the formation of such low-mass first stars around massive stars (Susa et al. 2014) or the formation of the first binary systems (Stacy et al. 2016;Sugimura et al. 2020).
The effort to find the first stars in our Galaxy has revealed hundreds of stars with [Fe/H] −3 (Bond 1980;Beers et al. 1985;Ryan et al. 1991;McWilliam et al. 1995;Ryan et al. 1996;Fulbright 2000;Norris et al. 2001;Carretta et al. 2002;Johnson 2002;Cayrel et al. 2004;Cohen et al. 2004;Honda et al. 2004; Barklem et al. 2005;Aoki et al. 2007;Lai et al. 2008;Caffau et al. 2013;Norris et al. 2013;Roederer et al. 2014b;François et al. 2018). These stars are called extremely metal-poor (EMP) stars; among them the currently most iron-poor star has [Fe/H] < −7 (Keller et al. 2014). The origin of EMP stars has been controversial mainly in explaining the large fraction of carbonenhanced metal-poor (CEMP) stars. Previous scenarios can be classified into three theoretical models, all of which try to explain the abundance patterns of CEMP stars. In the three major scenarios, CEMP stars are formed (1) by binary mass transfer from the former AGB stars (Suda et al. 2004), (2) from the gas in the interstellar medium (ISM), contaminated by the ejecta from the first-generation supernova (Umeda & Nomoto 2003), and (3) from the gas ejected from fast-rotating massive stars (Meynet et al. 2006).
In this paper, we explore another possibility for the formation of the first stars that can be verified by observations, i.e., the evolution of the first stars in a binary system consisting of a massive star and a low-mass star. In such a system, the ejecta of a supernova explosion in the massive star collide with the low-mass companion, which will either strip away or be gravitationally confined to the surface of the companion, or both may happen due to a wide range of the speed of the ejecta. This scenario provides a new pathway to look for evidence of polluted first stars among known halo stars in the Milky Way Galaxy, where there are more than 500 stars with detailed chemical abundances derived from high-resolution spectra.
The paper is organized as follows. In the next section, we provide the overview of our scenario. The details of the models are described in section 3. The results of our simulations are discussed in section 4. Section 5 provides the implications from our scenario. Conclusions follow in section 6.

New scenario to identify the survivors of the first stars
We propose that some low-mass first stars should have survived supernovae in binary systems. They can be observed in our Galaxy as metal-rich halo stars if they were contaminated by the supernova ejecta in close binary systems. A sufficiently small binary separation is required for surviving stars to change their surface chemical composition by the accretion of the ejecta. This is not necessarily guaranteed because typical massive stars evolve to red supergiants, having radii up to ∼1000 R before the core-collapse. If the separation between two stars is shorter than the radius of the evolved massive star, the binary system will undergo the common envelope phase where the companion star will be embedded in the envelope of the massive star. Such a system will have a short life span due to the mass transfer and cannot be a candidate for the first-star binaries. The first stars can avoid the evolution to red supergiants due to the initial lack of metals in the centre, more specifically CNO elements. The final radius of a metal-free massive star is reported to be of the order of 10 R (Heger & Woosley 2010;Tanikawa et al. 2020). This result is supported even if stellar rotation is taken into account (K. Takahashi & T. Yoshida 2020 private communication). . The data points are compiled from the abundance data from selected literature using the SAGA database (see appendix 1). However, we made corrections to the abundances by considering the effect of dilution by the convective envelopes in subgiants and red giants to compare models and observations with lithium abundances at the main sequence phase. The plotted data are classified and counted (shown by the numbers next to the labels) according to the evolutionary status as discussed in the text. The data points with circles denote the model results in table 1. (Color online) The upper boundary of the metal content of stars not to evolve into red supergiants depends on the initial mass, where the threshold value is [Fe/H] ∼ −7 for 12 M and ∼ −2 for 20 M (see subsection 3.1). It is to be noted that less massive stars are more abundant than more massive stars, and hence metal-free stars are dominant among massive stars with small final radii.
The evidence of the first supernova binaries can be checked observationally by the lithium and iron abundances of companion stars. Figure 1 shows the lithium abundances of known metal-deficient stars as a function of metallicity. Lithium is a good tracer of the change of the surface chemical composition because it is easily destroyed at layers with temperatures above 2.5 × 10 6 K, which means that the total lithium content is small. Also, the depletion or enhancement of lithium is easily identified by observations because the stars in the main sequence phase have a typical lithium abundance value of A(Li) = 2.1, which is called the Spite plateau (Spite & Spite 1982). Once the surface chemical composition is altered by external pollution such as stripping and/or accretion due to the collision of supernova ejecta, the surface lithium abundance is depleted by factors of several or more from the original value. Observationally, a non-negligible fraction of metal-poor stars shows lithium depletion, which is not explained by the standard stellar evolution models. In this study, we explore the possibility of lithium depletion.
Simulations of the collisions between supernova ejecta and a binary companion are performed using a 3D smoothed particle hydrodynamics (SPH) code, ASURA (Saitoh et al. 2008). The initial conditions for the simulations are calculated by a 1D Lagrangian hydrodynamics code calibrated by SN 1987A (Shigeyama & Nomoto 1990) so that the model reproduces the light curve of the supernova well. The initial masses and explosion energy are set at 15, 20, and 25 M and 10 51 erg, respectively. The remnant mass of the progenitor is assumed to be 1.3 M , following the results of Heger and Woosley (2010). Stellar mass loss before the explosion is ignored, following the previous studies (Heger & Woosley 2010;Tanikawa et al. 2020). The initial condition for the explosion is mapped on 3D simulations by setting the pressure, density, and kinetic energy in the primary star. The structure of the secondary star is constructed by mapping a 1D model star computed with a stellar evolution code (Suda & Fujimoto 2010). We assumed that the convective zone of the surface of the model mixes the materials homogeneously and does not change its mass and structure by the collision of supernova ejecta. This assumption may affect the estimate of the final chemical composition by the simulations. Stellar rotation is Downloaded from https://academic.oup.com/pasj/article/73/3/609/6249457 by guest on 28 July 2021 ignored for the sake of simplicity. In our simulation setups, binary components should be tidally synchronized so that the effect of rotation is not significant.
Only the separation is a relevant binary parameter in our simulations. This is because the orbital velocity is much smaller than the impact velocity of the ejecta and because the simulation time is much shorter than the orbital period. Here we presume that the eccentricity is close to zero to avoid the passage of a companion star into the envelope of a primary. We fixed the supernova model in this study for the following reasons: (1) the simulations by changing the velocity distribution and/or the total kinetic energy of the ejecta will be essentially the same as the models with different binary separations, (2) it will be a big computation cost to run more simulations with different inputs for supernova models, and (3) this project is still in embryo and it will be too much work to test unusual types of supernova such as asymmetric explosions and mixing and fallback models. The impact of the variations in supernova models is addressed later in this paper and such simulations can be a future work.
In the next section, we describe the details of the models and assumptions to perform numerical simulations and compare the results with observations.

Models and assumptions
3.1 Evolution of massive metal-free and metal-poor stars The characteristics of the evolution of metal-free and metal-poor stars are the smaller radii at the ends of their lives compared with stars that evolve through the red supergiant phase. Metal-deficient massive stars are supported by hydrogen-burning in the centre at the early phase without CNO elements (the p-p chain reactions), which requires higher temperatures than hydrogen burning with CNO catalysts (the CNO cycles). The higher the temperature of the nuclear burning regions, the faster the progress of the nuclear burning stages. Therefore, the final nuclear burning finishes before the stars pass through the Hertzsprung gap where the radii increase rapidly.
We have computed the evolution of massive stars with various mass and metallicity, which is displayed in figure 2. We use the same stellar evolution code (Suda & Fujimoto 2010) with the addition of carbon burning reactions. Cross-sections of 12 C + 12 C and 12 C + 16 O are taken from the fitting formulae provided in Caughlan and Fowler (1988). The computations are terminated when the mass fraction of carbon in the centre becomes smaller than The two symbols on the evolutionary tracks correspond to the onset of helium (smaller symbols) and carbon (larger symbols) burning, respectively. Computations were terminated at the end of the carbon-burning phase where the mass fraction of carbon abundance is below 0.02. (Color online) 0.02. The locations of the onset of the helium-and carbonburning phases on the Hertzsprung-Russell (H-R) diagram are shown by smaller and larger circles, respectively. For the sake of simplicity, carbon and oxygen in these reactions are converted to 25 Mg. This approximation is sufficient to see the final location of model stars on the H-R diagram.
In some models with low-metallicity stars, we found numerical instabilities caused by the hydrogen ingestion into helium-burning layers where the temperature is the order of 10 8 K. For instance, 15 M models with the metallicity of [Fe/H] ≤ −8 experience the inward extension of hydrogen-burning convection into the helium core before the onset of carbon burning. This results in a hydrogen flash with a hydrogen-burning luminosity exceeding log (L/L ) > 10, which does not provide any convergence in the stellar evolution. It is not clear why this mixing occurs, and it should be studied in a separate paper. The same phenomenon is reported in previous studies for rotating models with M = 160 M (Takahashi et al. 2018) and M > 200 M (Yoon et al. 2012). We tested models with different time steps and rezoning during these evolutionary phases, and we always encountered the same instabilities. We also tested other models using different stellar evolution codes and found the same phenomena.
To compromise the computations of these models, we artificially prohibited the mixing of the convection into the helium core. This can be justified for the evolution of the star itself because the overall evolution is controlled by the nuclear burning at the centre, not by the shell burning. The evolutionary tracks in figure 2 are apparently consistent with the previous studies (Heger & Woosley 2010;Tanikawa et al. 2020).
The computations were also made for the models of lowmass stars with zero and low metallicities. The model of M = 0.8 M with mass fraction, Z, of elements heavier than helium set at zero was constructed at the age of 10 Myr from the zero-age main sequence, which was used as a target star in the 3D simulations. Other models of M = 0.82 M with various metallicities were computed to check the effect of surface pollution on the low-mass stars. In the following, we present the models of 0.8 M and 0.82 M stars in the same plot, but there are only minor quantitative differences in these models. Figure 3 shows the evolution of the mass of the surface convective zones from the main sequence phase to the beginning of the red giant phase.
To confirm the effect of surface pollution, we computed the models by covering metal-rich (Z = 10 −4 and Z = Z ) materials in the surface convective zones. The depth of the convective zone is almost the same as the models without pollution because the structure of the surface convective zone is determined by the nuclear burning in the center.

Sample selection from the database
The data have been taken from the 2019 December 11 version of the Stellar Abundances for Galactic Archaeology (SAGA) database (Suda et al. 2008(Suda et al. , 2011(Suda et al. , 2017Yamada et al. 2021). The differences in the adopted solar We have removed sample stars that are likely affected by extra mixing during the RGB phase. There is increasing evidence of lithium depletion by extra mixing in red giants above the RGB bump (Charbonnel et al. 2020). These stars are removed from the sample by defining the boundary of the RGB bump, which is estimated by the comparison of stellar models with the observed effective temperature, luminosity, and metallicity. We identified the loca- the majority of extra mixing stars show lithium depletion after the correction with stellar evolution models (see figure 1).

SPH simulations
We conducted a series of hydrodynamical simulations of binary systems of the first stars, to estimate the evolution of the surface lithium abundances and metallicities. In this section, we describe the models, numerical methods, and some comparison results, which are necessary to select the fiducial models. Table 1 provides the model parameters in this study. The model H15F is the fiducial model in this study. The reason for the choice of our numerical setups can be found in the following subsections and the Appendix. The model names that use "S" (for "separation") are used to investigate the dependence on binary separation. The models named with "R" (for "resolution") refer to resolution studies, and the model with "SSPH" corresponds to the investigation with the standard SPH method. The model with "W," which means "whole," adopts the highest resolution in the whole simulation volume for the ejecta without using anisotropic particle-mass distribution, as described in sub-subsections 3.3.3 and 3.3.4. See the appendix for our feasibility studies on resolution and computational methods, corresponding to models named with "R," "W," and "SSPH."

Numerical simulation code
Numerical simulations shown here were carried out by a parallel code ASURA (Saitoh et al. 2008(Saitoh et al. , 2009, which was originally developed for the simulations of galaxy formation. The hydrodynamic equations are solved by a SPH method (Lucy 1977;Gingold & Monaghan 1977) with an improved version of the SPH [density-independent formulation of SPH; DISPH ] implemented. The conventional formulation of SPH assumes the smoothness of the density field and cannot handle fluid instabilities due to the unphysical surface tension that appears at contact discontinuities. We employed the simple equation of state (EoS) that assumes ideal and mono-atomic gas (i.e., the specific heat rate is γ = 5/3). Chemical reactions in the ejecta and the star are not taken into account. Self-gravity is solved by the Tree method (Barnes & Hut 1986). We only considered the monopole moment with the tolerance parameter θ = 0.5. The interactions between particle-particle and particle-monopole moments are computed with Phantom-GRAPE (Tanikawa et al. 2012(Tanikawa et al. , 2013. 2.97 * The columns represent the model name, the number of SPH particles for a companion star and for a supernova ejecta, the binary separation, the unbound mass by the stripping, the bound mass by the accretion, and the final iron and lithium abundances without and with lithium-rich ejecta in the surface of the companion star. See section 5 for the details of estimating the abundances. † The accreted mass is not estimated because we stopped the simulation too early to estimate the effect of accretion. Figure 5 shows the schematic picture of our model of a binary system. It consists of a low-mass star and supernova ejecta. 1D models are mapped on to the particle distributions in a 3D volume. We fixed the mass of the low-mass star and three different supernova ejecta named H15, H20, and H25 which represents 15, 20, and 25 M , respectively, from Heger's pre-supernova models (Heger & Woosley 2010). We tested four different separations to investigate the dependence of the stripping and accretion of the surface of the low-mass star by the collision of ejecta. The impact of the initial separation on the evolution is studied for the H15 model by adopting four different separations. The separations of the H20 and H25 models are the minimum allowable in which the low-mass star and supernova progenitor contact with each other. The minimum separation in the H15 model also has this configuration (see sub-subsection 3.3.5 for further details). We assumed that the convective zone of the surface of the model mixes the materials homogeneously and does not change its mass and structure by the collision of supernova ejecta. This assumption may affect the estimate of the final chemical composition by the simulations. We ignored the rotation of a low-mass star with respect to the supernova progenitors because its timescale is much longer than those of the expanding ejecta. In our simulation setups, binary components should be tidally synchronized so that the effect of rotation is not significant.

Mapping the 1D stellar model on to the 3D space
We used the result of the numerical simulation of 0.8 M with Z = 0 as described above. Figure 6 shows the radial profiles of the 3D model star taken from the 1D model at t = 10 Myr from the zero-age main sequence.
To construct a 3D model with particles, we put the particles in individual shells so that the mass contained in (r in , r out ) satisfies the condition rout rin where ρ(r) is the density profile and m(r in ) is the particle mass at r in . The exact position of a particle in the shell is determined by R(r out − r in ) + r in with a random number R ∈ [0, 1]. The values of the density at individual positions are given by linear interpolation between ρ(r in ) and ρ(r out ). The angular positions, the polar and azimuthal angles, are also randomly chosen. The central particle is distributed by finding r out from the above equation for m(r in ) = 0 and m(0) = 32m p, * , where m p, * is the minimum mass of an SPH particle and is set at 3 × 10 −7 M . For the rest of the mass distribution m(r), we adopted the following formula to reduce the computational cost: where R edge ≈ 0.64 R . The outer layers have finer mass resolutions. The total number of gas particles is 1023991. The star has the chemical composition of X = 0.767, Y = 0.233, Z = 0, where X, Y, and Z represents the Downloaded from https://academic.oup.com/pasj/article/73/3/609/6249457 by guest on 28 July 2021 The circle at the top is the low-mass star companion in the binary system. The arrows represent the ejection of matter from the massive star. The initial separation between the two stars is measured by the distance between the centre of the supernova ejecta and the low-mass companion. (b) The initial structure of the star. The inner region consists of coarse particles to reduce the computational cost. (c) The initial structure of the ejecta. The ejecta are divided by the direction so that the region has the finest resolution for colliding angles.
mass fraction of hydrogen, helium, and other elements, respectively, at the age of 10 Myr from the onset of hydrogen burning in the center. The internal energy of each particle is computed by the linear interpolation of the temperature profile with the assumption of an ideal gas with the adiabatic index of γ = 5/3 and the mean molecular weight of 0.6.
Before starting the simulations of ejecta collisions, we relaxed the particle distribution of a 3D model star following the simulations for Type Ia SNe (Rimoldi et al. 2016). We added the dumping term in the hydrodynamical force and integrated the star for 2 hr, which is sufficient for the relaxation because the dynamical time of the system is estimated to be 30 min. This technique enables us to avoid the oscillations of the stellar surface, which is not suitable for our simulations due to the change of the cross-section of collisional particles.

Mapping the 1D ejecta model on to the 3D space
The 1D profiles of pre-supernova models with Z = 0 are taken from the literature (Heger & Woosley 2010). We have employed the 15, 20, and 25 M models. Figure 7 shows the radial profiles of the three supernova ejecta models. The outer edges of ejecta are 10, 14, and 24 R , and the corresponding times from ignition are 1747, 2195, and 5705 s, respectively. The model of 15 M is prescribed in the model of SN 1987A (Shigeyama & Nomoto 1990).
We have developed a new technique to ensure a high resolution in colliding particles by changing the particle mass depending on the direction of the ejecta particles.
Here we define the direction of the ejecta particles by the angle θ from the centre of the ejecta along the z-axis connecting with the centre of the companion star. The particle mass depends on the initial location of the ejecta on the (θ , z) coordinate as follows: where m p,ej is the minimum mass of an ejecta particle. The positions of the particles are determined in the same way as described in the previous subsection. The radial velocities, internal energies, and densities for individual particles are assigned to reproduce the profile in figure 7. We tested four different choices of m p,ej , i.e., 1 × 10 −6 , 3 × 10 −7 , 1 × 10 −7 , and 3 × 10 −8 M . The third choice (H15) is adopted as the fiducial resolution because that result converges with the result using the highest mass resolution (H15Rc). The corresponding numbers of particles for the four models are 479442, 1599777, 4800569, and 16001929, respectively. Those for H20 and H25, with the fiducial mass resolution, are 6487398 and 8189647, respectively. We have confirmed the validity of using particles with different masses by running a test simulation using equal-mass particles and compared it with the directiondependent mass model for the same condition for the collision. It is confirmed that the stripped mass from the companion star and the accreted mass are consistent with each other. We found that we need ≈ 140000000 particles for the equal-mass model to achieve a similar number of effective particles (tan θ < 0.25, z ≥ 0), i.e., our new technique is successful in reducing the number of particles by 1/30. Thanks to our improved method for the assignment of mass in the SPH particles, more than 130000 particles effectively interact with the companion star despite the small visual angle of 3. • 5 in the H15 model.
The consistency between 1D and 3D models is confirmed by the total kinetic energy just before the impact of the ejecta. We also checked the stripped mass of the companion with the setup of Type Ia supernovae and compared it with the previous study under the same initial mass and separation (Pakmor et al. 2008).

Initial conditions for a binary system
In the fiducial model, we adopt the initial separation of 0.1 au (21.5 R ), which is larger than the radius at supernova explosion according to metal-free star models (Heger & Woosley 2010) and the critical radius of the Roche-lobe overflow based on the empirical formula (Eggleton 1983). We also computed the cases with 0.2, 0.4, and 0.8 au to investigate the dependence on the initial separation. To compare various progenitors, we computed the models of 15, 20, and 25 M stars with the minimum separations that are characterized by the size of the progenitor stars. These correspond to 0.048, 0.064, and 0.112 au.
The orbital motion of the binary system was ignored in this study for the sake of simplicity. This is reasonably validated by the short timescale (∼10 hr) in our simulations compared with the orbital period (a few days to a few months). The consideration of the orbital motion may have some effect on the amount of accretion at the separation of 0.1 au, while the stripping of the surface layers will not be affected because the stripping is dominated by massive fast ejecta in a shorter timescale (a few hours). It is worth investigating the case of accretion with the motion of a binary orbit taken into account, but this is beyond the scope of this paper.
By ignoring the centrifugal force by the orbital motion, a gravitational pull exerts on the ejecta and the star, which Downloaded from https://academic.oup.com/pasj/article/73/3/609/6249457 by guest on 28 July 2021 results in an unwanted motion in the system. Therefore, we filtered a long-range force of gravity for the companion star by ignoring the gravitational interaction with the ejecta at the distance of more than 10 R from the centre of the star. This prescription mimics the orbital motion and is sufficient to avoid the strong gravity caused by the ejecta at the initial position and to consider the effective accretion process after the collision.
Mass loss from the progenitor stars are not considered, which is justified by the metallicity dependence of stellar winds.
The rotation of the star and the ejecta is not considered. The consideration of rotation would be more realistic, especially for the ejecta, but it is too complicated to implement in the simulations. On the other hand, fast rotation in companion stars will not be mandatory in the binary systems studied.
We also employed a particle-split method (Kitsionas & Whitworth 2002;Martel et al. 2006) to retain a high resolution even in large separations. In the run of the H15Sa model, the ejecta particles within 6 R from the center of the low-mass star split into eight smaller particles. In the runs of the H15Sb and H15Sc models, those within 12 R split into eight smaller particles, and those smaller particles divided into eight further smaller particles when they were within 6 R . Figure 8 shows the result of simulations for a binary system consisting of 15 and 0.8 M stars separated by 0.048 au (the H15 model in table 1). The impact of the ejecta on to the binary companion is simulated for 24 hours after the explosion for this model. Once the ejecta material collides with the surface of the companion star, a bow shock is formed near the surface. The basic picture of the outcome of the collision is consistent with previous studies (Pakmor et al. 2008;Liu et al. 2013;Hirai et al. 2018).

Results
There are two possible impacts on the surface of a binary companion, induced by the collision of supernova ejecta. Fast-moving outer ejecta, with a velocity as much as 10000 km s −1 , strip the surface layers of the companion star, while a part of slowly moving inner ejecta, down to 3000 km s −1 , is accreted by the companion. It is found that the ejecta do not strip the whole of the surface convective zone. Therefore, the surface chemical composition does not change due to this effect (top panel in figure 9). The mass of the accreted matter is estimated by calculating the total energy of individual particles representing the ejecta (bottom panel in figure 9).
We found that the particles near the surface of the companion stars experienced a shock heating that gave rise to high temperatures with 10 8 K. At this temperature,  7 Li can be destroyed quickly even for this short timescale event. The density of the shocked region is the order of 10 −4 g cm −3 , which is lower than the density, ∼1 g cm −3 , of the shell where 7 Li burns in the envelope with T ≈ 3 × 10 6 K. On the other hand, the nuclear reaction rates are much larger at 10 8 K than at 3 × 10 6 K, by 16 orders of magnitude (Caughlan & Fowler 1988). The nuclear timescale for lithium destruction by proton-capture reactions is estimated to be ∼300 s. However, the particles affected by shock heating do not accrete on to the surface of the companion star and escape from the system. Therefore it is justified to ignore the change of the surface chemical composition by nuclear reactions in our hydrodynamical simulations. Figure 10 shows the time evolution of the mass of bound ejecta particles grouped by chemical compositions. The inner ejecta, such as the nickel-and silicon-rich layers, result in more efficient accretion than the outer ejecta like the hydrogen-and helium-rich layers. This can be interpreted as strong velocity dependence of the accretion rate, where the amount of accretion is larger for slower ejecta. If the ejecta shells were unmixed and retain their chemical composition, the ejecta shells within the oxygen, neon, and magnesium layers would be accreted on to the companion. This means that the final surface abundances of companion stars would be very metal-rich with unusual abundance patterns. However, supernova ejecta are thought to be well mixed during the shock propagation in the envelope by the Rayleigh-Taylor instability, inferred from the timing of the detection of line gamma-rays and X-rays in SN 1987A (Arnett & Fu 1989;Kumagai et al. 1989). It is to be noted that there are multi-dimensional simulations of the SN 1987A explosions (Utrobin et al. 2019), but the yields of supernova ejecta are not available from their results. Although our models are not able to predict the abundance patterns, it will be interesting to specify the progenitors by characteristic abundance patterns in the first supernovae. If the observed stars reflect abundance patterns by core-collapse supernovae, we may observe deficiencies of odd-Z elements. Also, the abundance ratios of [Na/Mg] or [Ca/Mg] can be diagnoses for pair-instability supernovae (Takahashi et al. 2018), which can be used to discriminate our scenario from the second-generation stars formed out of the ejecta of pair-instability supernovae.
The accretion process is subject to uncertainties in estimating the bound mass. We estimated the expected amount Downloaded from https://academic.oup.com/pasj/article/73/3/609/6249457 by guest on 28 July 2021  of the bound mass by considering the balance between the gravitational binding energy and the kinetic energy of the individual particles of the ejecta. This is because the accretion process is complicated and the phenomenon happens over a much longer timescale than the time of the simulations. To estimate the accurate effect of accretion, we iterated to compute gravitational force exerted on the individual particles after removing the particles that are going to escape from the system, until no particles escape any more.
It is highly uncertain how the accreted gas particles mix with the convective envelope of the companion star. In our simulations, we do not follow the accretion process, which proceeds in a thermal timescale, due to a much longer timescale than the dynamical process in this study, and the technical difficulty to simulate gas accretion on to a stellar surface. The detailed modelling of the accretion process is poorly understood due to the complexity of the formation of accretion discs and the considerations of radiation pressure, magnetic force, and other physical processes on a stellar surface (Marietta et al. 2000). It is still to be established how the convective envelope is reconstructed after the gas accretion. There are no adequate prescriptions to implement the mixing of matter into the convective envelopes either in 3D hydrodynamical simulations or in a 1D hydrostatic stellar evolution code.
To consider the effect of mixing suggested by observations as above, we assumed that the bound ejecta particles are well mixed in the convective envelope of companion stars. We obtained [Fe/H] = −0.2 for well-mixed ejecta in the case of model H15 where the maximum amount of accretion is achieved. The lithium abundance will be reduced by 0.86 dex, from the estimate of the accreted mass and the mass of the surface convective zone for a metal-free 0.82 M model at the time of impact, when 10 7 yr have passed since the onset of hydrogen burning. For longer separations, lithium depletion is much less. In the 15 M case, lithium depletion is negligible at 0.2 au. In larger separations, the amounts of stripping and accretion decrease very rapidly. The results for all the models are summarized in table 1. The lithium abundances will not change after the accretion event since the surface convective zone of low-mass stars recedes in mass during the main sequence evolution (figure 3) and will not mix with the lithium-containing layer.

Discussions
The values of [Fe/H] mixed in table 1 are calculated by the following equation: where M acc denotes the accreted mass taken from the simulation results in table 1. Other parameters are taken from stellar models: M ini , M rem , and M conv is the initial mass, the remnant mass which is set at 1.3 M , and the mass of the surface convective zone of the low-mass companion which is set at 3.73 × 10 −3 M , respectively. The iron yields from supernovae, Y Fe , are taken from the mixed models with the explosion energy of 1.2 × 10 51 erg in table 6 of Heger and Woosley (2010). The abundance parameters, X Fe, and X H, , are taken from the literature (Asplund et al. 2009). The hydrogen abundance of X H in the surface of low-mass first stars is assumed to be the solar value. The final lithium abundances are estimated by the following equation: where M(Li, A Li , Y Li , M ej , and A Li, ini denote the lithium mass in the surface convective zone, the mass number of 7 Li, lithium yields from a supernova, the total mass of supernova ejecta, and the initial lithium abundance, respectively. The values of Y Li and M ej are taken from the same models as above. The initial value of lithium abundance is assumed to be 2.1. The free parameter f depicts the contribution of the lithium yield in the supernova ejecta. In table 1, we provided the final lithium abundances without lithium yields, A(Li) woLi , using the above equation by setting f = 0 and those with lithium yields, A (Li) wLi , by setting f = 1. It is to be noted that the above two equations do not include the stripped mass. This is because the mass of lithium-containing layers is always smaller than the mass of the surface convective zone of the low-mass companion in our simulations. For instance, the interaction between the ejecta and the companion can strip the envelope with a mass of 3.3 × 10 −3 M at most, as shown in figure 9 and table 1. The stripped mass is smaller than the mass of the surface convective zone of a 0.82 M star at the age of ∼10 Myr, which is the lifetime of a 15 M star. In this case, the stripping event will not reduce the surface lithium abundances, although it depends on how we treat the mixing in the surface convective zones, as there are no reliable models or theories on it. Here we expect that the surface convective zones will be reconstructed after the accretion event rather than the stripping event. The duration between the stripping and accretion is much shorter than the time for the surface convective zone to reconstruct. Therefore, the accretion events of ejecta overwrite the effect of the stripping.
Our models predict solar-metallicity stars in the halo of our Galaxy for shorter binary separations. Figure 1 shows the comparisons of lithium abundances. Provided that the first massive stars explode with small radii (figure 2), it is more likely that solar-metallicity stars in the halo are survivors of the first stars in the universe, originally having a chemical composition as a result of big-bang nucleosynthesis. In the case of longer binary separations, our models may account for metal-poor ([Fe/H] ∼ −3) lithium-depleted stars and lithium-rich [A(Li) > 3] stars. Metal-poor ([Fe/H] < −5) stars without lithium depletion are also possible, although we cannot distinguish between our binary scenario and the case with single stars. The abundance patterns of all the measured elements should be checked individually, but they will vary according to the abundances of yields and the degree of mixing in the ejecta. The observed lithium depletion around [Fe/H] = −3 is not explained by our models due to a small amount of lithium depletion for longer separations and too much accretion for shorter separations. Instead, our models predict a small or negligible depletion at [Fe/H] < −5 and the existence of solar-metallicity stars in the halo stars in our Galaxy.
To confirm the existence of the first metal-rich stars in the Galactic halo, we need to exclude the possibility of normal metal-rich halo stars. In our proposed scenario, we have excluded the possibility of stars that survived the common-envelope phase because a companion star in the common envelope phase will accrete sufficient mass to become much more massive than 0.8 M . Such stars can be easily excluded from the position on the H-R diagram. Even if there is a Wolf-Rayet star with a low-mass companion, its supernova ejecta must be too fast to be accreted on to its companion due to the absence of the hydrogen-rich envelope in the progenitor. There is also an argument that the evolutionary path of the SN 1987A progenitor is not established. The most promising scenario for the progenitor is that it is a binary merger (e.g., Morris & Podsiadlowski 2007). If this is the case, we need a third star to make our scenario work. This is obviously an unusual case and such stars cannot be a major source of noise in the search for metal-rich halo stars.
The formation of solar-metallicity stars in metal-free or metal-poor host halos could be another exceptional case to consider. According to the hydrodynamical simulations of star formation triggered by supernovae with initial conditions taken from cosmological simulations, secondgeneration stars after first supernovae have metallicities up to [Fe/H] ∼ −1 (Wise et al. 2012;Smith et al. 2015;de Bennassuti et al. 2017;Chiaki & Tominaga 2020). It is not realistic that solar-metallicity stars would form with a single supernova in a metal-free or metal-poor cloud in its host cloud. However, it may be possible if the kinetic energy of supernovae is very low and the swept-up mass is correspondingly small. It would be intriguing to compare the abundances of second-generation stars formed from Downloaded from https://academic.oup.com/pasj/article/73/3/609/6249457 by guest on 28 July 2021 low-energy supernovae with those of binary companion stars in close massive binaries.
There is a possibility of the formation of solar-metallicity stars in chemical evolution in the progenitors of the Galactic halo. Simulations of metal-enrichment in the Galactic halo by early generation stars predict metallicity of up to [Fe/H] ∼ −0.5 (Komiya et al. 2014;Sarmento et al. 2017;Côté et al. 2018). For example, Côté et al. (2018) found star-forming gas with solar-metallicities at the redshift of z = 7.29 by their simulations of the formation of the most massive galaxy realised by Wise et al. (2012). Sarmento et al. (2017) found a highly-polluted region with gas metallicity of Z = 0.1 Z close to the centre of its host halo by their simulations up to z = 5. These results may imply that the current chemical evolution models do not focus on the formation of the most metal-rich stars in the Galactic halo, i.e., we presume that solar-metallicity stars would not form or survive in the halo. If metal-rich stars are formed in the halo by consecutive star formation, their abundance patterns should be dominated by type II supernovae, not by type Ia supernovae. We will be able to identify the progenitors of such stars by larger [α/Fe] compared with disc stars with the same metallicity, if they exist.
It is difficult to follow the formation of close binary systems in a metal-free environment using hydrodynamical simulations. Since our scenario focuses on binaries with a separation of less than ∼1 au, we need simulations with very high resolution. All of the current simulations of first-star binaries deal with protostar binaries much wider than 1 au (see, e.g., Sugimura et al. 2020;Susa 2019). It is not clear if the formation of first-star binaries with small separations are supported by numerical simulations. However, it is interesting that higherresolution studies found more possible close binaries by decreasing the minimum spatial resolution of the simulations from 20 au (Stacy & Bromm 2013) to 5 au (Stacy et al. 2016).
We anticipate more applications of our scenarios to the origin of known peculiar stars, namely lithium-enhanced stars. There is an argument concerning the production of lithium in the hydrogen-rich envelope of massive stars by neutrino processes (Heger & Woosley 2010). If we simply apply the well-mixed lithium yield to the accretion on to the low-mass companion, the maximum surface lithium abundance can be A(Li) ∼ 3.0. This may account for some of the extremely lithium-enhanced stars as large as A (Li) ∼ 4.0 among the main-sequence stars (Li et al. 2018). Such lithium-enhanced stars are very rare, comprising 0.1 percent of the total population in almost the entire metallicity range (Kirby et al. 2016). This fact implies that the ejection of lithium-rich ejecta and the accretion on to the companion involve complicated physics.
Theoretical models are not able to explain lithium enhancement in the surface of low-mass stars. The transportation of 7 Be and the decay to 7 Li in the envelope require a high temperature at the bottom of the convective envelope, and hence only AGB stars with M 4 M can be A (Li) 3.0 by self-enrichment (Karakas & Lugaro 2016). Low-mass red giants are thought to suffer from extra mixing at the red giant branch bump (see, e.g., Charbonnel et al. 2020) and experience strong lithium-depletion. Even if the parameter for extra mixing is adjusted, it is too difficult to explain lithium-rich giants (Lattanzio et al. 2015). Nova outbursts produce a large amount of lithium (Starrfield et al. 2020). However, there are no known mechanisms which form low-mass stars in the Galactic halo from yields from novae. If there exist metal-rich and lithiumrich halo stars whose ingredients are dominated by nova yields, we may distinguish them by the abundance ratios of CNO isotopes (José & Hernanz 1998). Lithium enhancement in super-solar metallicity is also controversal. Theoretical models predict that lithium abundances are enhanced at super-solar metallicity, while the observations reveal a decreasing trend of lithium abundances. According to the models by Karakas and Lugaro (2016), AGB stars with M ≥ 4.75 M and Z = 0.03 evolve to lithium-rich stars up to A(Li) ∼ 5, while the observations show the decrease in A(Li) at [Fe/H] > 0 (Guiglion et al. 2016). Chemical evolution models for lithium in the Galactic disc can reproduce lithium abundances as high as A(Li) ≈ 3.5 (Prantzos et al. 2017). Therefore, the only possible contaminants are metalrich stars born in the Galactic disc that fly out into the halo.
An exploration of a different set of parameters and assumptions could be of interest. For instance, an assumption of a spherically symmetric explosion is not guaranteed. Aspherical explosions or jet explosions will change the amount of stripping and accretion (Tominaga 2009), while there is no convincing theory and observations for the geometry of supernova explosions. Still, there will be a way to produce metal-rich survivors with these models. A more eccentric orbit will have a variation in the stripping and accretion depending on the orbital phase at the explosion, although this effect can be incorporated into the binary separation. It is difficult, and beyond the scope of this study, to estimate the overall uncertainties and these new parameters. Various situations and models may better explain the diversity of the chemical compositions of known extremely metal-poor stars.
After the supernova event of the massive star, the binary system will be disrupted by the ejection. The remnant of the supernova will have a velocity of a few 100 km s −1 (Cordes et al. 1993). In the circumstance of the first binary star formation, it may be difficult for the remnants of the supernovae to be confined in the gravitational potential of the low-mass host halo. The low-mass companion star can remain in its host halo if the orbital velocity is lower or comparable to the velocity dispersion of the host halo.
Currently, there is no direct observational evidence of metal-free and metal-poor binary systems as we proposed because they are too old for massive primaries to survive. Therefore, the search for observational counterparts must be made for nearby, metal-rich stars. We searched for low-mass companions in OB star binaries (Moritani et al. 2018). The analysis of the 10 target stars is still ongoing, among which HD 164438 is reported to have an intermediate-mass star with the mass ratio of q = M 2 /M 1 = 0.1-0.2 (Mayer et al. 2017), where M 1 and M 2 are the initial masses of the primary and secondary, respectively. Although we have not yet established the population of binary systems in the solar vicinity, we can expect some fraction of stars that are metal-rich counterparts of the first-generation massive binaries. There are no reliable models or observations of binary star formation for the entire range of metallicity for OB stars, and hence we encourage the exploration of the binary star formation.
To explore the possibility of observational counterparts in the solar vicinity, we have compiled a catalogue of OB stars from the literature to find binaries where the masses of the primary star and the secondary star are 10-20 M and less than 1 M , respectively, with a separation less than 1 au. This criterion is translated into a mass ratio q < 0.1 and a binary period of a few days to a few months. We investigated more than 20 binary catalogues and surveys to compile the existing observational data. These data include spectroscopic binaries (Pourbaix et al. 2004), eclipse binaries (Malkov et al. 2006), astrometry binaries (Mason et al. 2001) and visual binaries (Sixth Catalog of Orbits of Visual Binary). 1 We have collected binary parameters for the stars classified as Galactic OB stars using the online catalogue (Skiff 2014). Figure 11 shows the distribution of mass ratio taken from our compilation of 332 OB stars with known binary mass ratios and periods (1 yr or shorter) from the observations. There is an observational bias in the distributions of both periods and mass ratios in our sample. In particular, binaries with small mass ratios are difficult to identify, and hence the number of such binaries is underestimated. In the figure, no corrections for the distribution are considered, as discussed in other studies (Moe & Di Stefano 2017).
Apart from the literature study, we have conducted a radial velocity monitoring to specify binary systems that are the metal-rich counterpart of those relevant to the massive binary scenario. Medium-resolution spectra have been obtained with the Medium And Low-dispersion 1 http://www.usno.navy.mil/USNO/astrometry/optical-IR-prod/wds/orb6 . Fig. 11. Mass ratio distribution function from the literature data for the binaries whose periods are 1 yr and shorter.
Long-slit Spectrograph (MALLS) (Ozaki & Tokimasa 2005) on the Nayuta 2.0 m telescope at the Nishi Harima Astronomical Observatory (NHAO). Also, high-resolution spectra has been obtained using two spectrographs, High Dispersion Echelle Spectrograph (HIDES) (Izumiura 1999;Kambe et al. 2013) on the 188 cm telescope at the Subaru Telescope Okayama Branch Office (OBO), and the Gunma Astronomical Observatory Echelle Spectrograph (GAOES) (Hashimoto et al. 2006) on the 1.5 m telescope at the Gunma Astronomical Observatory (GAO). The details of the observations and analyses are on-going and an initial report is given in a separate paper (Moritani et al. 2018).
Candidates for the survivors of the first stars exist in our Galaxy. A massive binary scenario presented here predicts the existence of metal-rich ([Fe/H] ∼ 0) and lithium-depleted, or perhaps lithium-rich [A (Li) ∼ 3], stars in the Galactic halo. We find 78 stars with [Fe/H] > −1 and six stars with [Fe/H] > −0.5 that have the space velocities of halo stars in the SAGA database, although none of them have reported lithium abundances (table 2). Among them, one star, BS 17587-021 is apparently a metal-rich halo star from its metallicity ([Fe/H] = 0.93, Lai et al. 2008) and its motion is different from the Galactic disc components, estimated using the data from the Gaia Data Release 2 (DR2) (Gaia Collaboration 2018). We need a more careful inspection of this star because only one iron line at 670.015 nm was used to determine its metallicity. There are also a number of solar-metallicity stars among halo stars kinematically selected to search for high-velocity stars ejected from the Galactic thick disc (Hawkins et al. 2015). The majority of such halo stars with −0.7 < [Fe/H] < −0.2 are argued to be the remnant of the thick disc component as a result of the last major merger event (Belokurov et al. 2020). These studies imply that the stars with [Fe/H] > −0.2 have different origins and have experienced a very efficient metal enrichment process, like the scenario  proposed here. However, it is unlikely that such an event produces even more metal-rich stars than the thick disc stars.
To identify the observational counterparts of the proposed scenario, we performed a cross-match of Galactic stars between the SAGA database and Gaia DR2 (Matsuno et al. 2019). Figure 12 shows the sample stars on the Toomre diagram that have space velocities larger than 180 km s −1 . The velocity data are calculated from the proper motion data taken from Gaia DR2 (Gaia Collaboration 2018). The boundary between the disc and halo stars are defined by the heliocentric velocity, V total = √ U 2 + V 2 + W 2 = 180 or 220 km s −1 as shown by the dashed lines in the figure, following the prescription in the previous studies (Nissen & Schuster 2010;Bonaca et al. 2017). The positions of the stars are calculated based on the right ascension and declination and the distance from parallax in Gaia DR2 where we imposed the criterion for the data accuracy by parallax_over_error >5. The total number of data with V total > 180 km s −1 is 1332. The massive binary scenario is a supplementary scenario to the existing three major scenarios for the origin of the known EMP stars, i.e., the scenario is consistent with the current framework. The three scenarios-(1) the binary mass transfer scenario, (2) the mixing and fallback in the first supernovae scenario, and (3) the fast-rotating massive star scenario-can be tested by surface lithium abundances, especially in scenarios (2) and (3). The degree of the contribution of a previous generation star to the abundance patterns of an observed star is measured by the dilution factor (Meynet et al. 2010), which is defined by the mass ratio of the ISM to the ejecta in the previous-generation star. The massive binary scenario can reproduce a variation in the dilution factor, by estimating the ratio of the envelope mass to the accreted mass. Therefore, the models of peculiar supernovae and possibly rotating massive stars may fit with the proposed scenario to explain all the abundance patterns, including lithium. In particular, the effect of winds from rotating massive stars is worth investigating because the winds are much slower than supernova ejecta and hence the companions can accrete more materials from the winds.
A connection with the binary mass transfer scenario is intriguing. One of the most iron-poor stars, HE 0107−5240, belongs to a binary system whose binary period is more than 10 yr (Arentsen et al. 2019), which is consistent with the theoretical prediction of the binary mass transfer scenario (Suda et al. 2004). If long-period binaries are responsible for CEMP stars at low metallicity, a variety of CEMP stars and other metal-poor stars are expected, depending on binary parameters. Also, we may imagine a more complex scenario for the origin of CEMP stars. If we try to apply the massive binary scenario to HE 0107−5240, we will need a triple system to complete the abundance pattern. The accretion of iron and lithium, initiated by the collision of a supernova explosion on to the observed stars, is followed by the accretion of carbon and possibly s-process elements, keeping the abundance of lithium and iron-group elements on the surface. Considering the third star to reproduce all the abundances, we need a fine-tuning of model parameters, while the formation of a triple system is not so rare, as seen from the observations of nearby main-sequence stars (Moe et al. 2019).
The formation of binary systems in the early universe will modify the current framework of the history of star formation and chemical evolution in the early universe. So far, only the first supernovae are thought to have contributed to the early chemical enrichment. There are many arguments for the formation of second-generation stars associated with the simulations of the first supernovae in the early universe. The major concerns are the metal mixing of supernova ejecta with the ISM (Whalen et al. 2008;Smith et al. 2015;Chiaki & Tominaga 2020), the formation of the low-mass, second-generation stars with dust cooling (Yoon et al. 2016;Chiaki et al. 2017), the initial mass function of the first-and second-generation stars (Hirano et al. 2014), and the chemical yields from the first supernovae (Heger & Woosley 2010;Tominaga et al. 2007). All of these topics will influence the proposed scenario in favor or disfavor of it.
Our proposed scenario also influences the argument of the initial mass function of the first-generation stars. Our scenario is consistent with most of the recent simulations on the formation of first-generation stars in which massive stars are dominant. It will be more important to consider the formation channel of low-mass companions around massive stars. It is in line with the argument that some low-mass first-generation stars are formed in the disc of the central massive stars (Susa et al. 2014). The formation of intermediate and massive stars in the early universe is also favored to explain the large fraction of CEMP stars among EMP stars under the binary mass transfer scenario (Komiya et al. 2007). It is proposed that massive stars dominate the EMP population as a byproduct of the need for many intermediate-mass stars to pollute low-mass companions.
The proposed scenario may provide an important clue to the early chemical enrichment of neutron-capture elements, like barium in EMP stars. We argue that the first massive stars are the first polluters with metals of lowmass companions, which may include the enrichment of neutron-capture elements. Observationally, the majority of observed EMP stars exhibit detectable barium absorption lines. This means that EMP stars will possess the evidence of neutron-capture processes such as the r-process and the s-process in any metallicity range. Recent observations and theories prefer neutron-star mergers as the site of the r-process (Abbott et al. 2017) in reproducing the early chemical enrichment of neutron-capture elements in dwarf galaxies in the Local Group (see, e.g., Ji et al. 2016;Roederer et al. 2016 for observations andHirai et al. 2015;Safarzadeh & Scannapieco 2017;Tarumi et al. 2020 for simulations). However, it is not sufficient to reproduce the even lower abundances of neutron-capture elements at the lowest metallicity range (Tsujimoto et al. 2017). Because the s-process is not expected in supernovae, we may need the production of r-process elements in supernovae as argued earlier (Wanajo 2013).

Conclusions
We have proposed a new scenario to find the survivors of the first stars in the local universe through the simulations of the collision between the ejecta of a first supernova and a low-mass star in a close binary system. The supernova ejecta do not have a significant effect on the stripping of the surface layers of low-mass stars at a distance of 0.1 au. The effect of the accretion of ejecta on to the companion star strongly depends on its binary separation. If the separation is small enough, namely 0.1 au, the companion star can be very metal-rich, although it depends on the chemical composition of the yield and how the supernova ejecta mix with the convective envelope of low-mass stars.
We have argued that survivors of the first stars can be found among metal-rich stars in the Galactic halo. This is due to the small radii in the first massive stars, which enables them to contaminate low-mass companion stars with supernova ejecta in binary systems. The existence of these stars in former massive binaries is in line with some theoretical models of the formation of the first stars. We also investigated the literature dealing with binary stars for OB stars in the Galactic disc to confirm the existence of short-period binaries consisting of massive stars and lowmass stars. Although nearby OB stars are metal-rich and the formation mechanism of binary systems must be different, we find potential candidates of the metal-rich counterparts of the progenitor binaries that match with our proposed scenario. We also looked for metal-rich halo stars in the sample of known Galactic stars with kinematic information taken Downloaded from https://academic.oup.com/pasj/article/73/3/609/6249457 by guest on 28 July 2021 from the Gaia DR2. There are several stars with [Fe/H] > −1 that have space velocities of a typical halo population.
The available data suggest that we may already have a number of the first stars in our Galaxy, the surfaces of which are contaminated by the ejecta of the first supernovae. The exploration of various possibilities of binary formation in the early universe, first supernovae and/or other sources of ejecta from massive stars, like a wind from a fast-rotating massive star, and dynamical process of the collision between the ejecta and stars, is encouraged. The detailed chemical abundances and kinematic information of metal-rich halo stars will provide a better constraint on the proposed models.

Funding
This work is supported by a Grant-in-Aid for Scientific Research (KAKENHI) (JP16K05287, JP15HP7004, JP16H02166, JP16H06341, JP19HP8019, JP20HP8012) from the Japan Society for the Promotion of Science. Figure 13 shows the time evolution of the mass of lithiumcontaining layers, which represents the unbound mass by the collision of supernova ejecta, for four models with different mass resolutions in their ejecta. The three best highresolution runs (H15W, H15Rc, and H15F) produce very similar results, and hence the H15F model is selected as a fiducial model. We attempted using the highest-resolution run (H15W), but it took too much time, and we stopped computation at ∼4 hr into the simulation time. The mass of lithium-containing layers for the H15W model agrees well with the other two high-resolution models, although the comparison can be made only for the early phase of the simulation. This confirms that anisotropic particle-mass distribution is a good prescription to simulate the collision of ejecta with a target object. In the lowest-resolution case, the unbound mass is larger compared to the other runs, which may come from the coarse sampling of the ejecta. Therefore, we adopted the run with N p,ej = 4800569 and m p, ej , i.e., 1 × 10 −7 M , as the fiducial resolution in this study.

Appendix 3. Dependence on the SPH scheme
Previous studies point out that the conventional formulation of SPH is not suitable to handle the contact discontinuity and is not able to detect some fluid instabilities (Agertz et al. 2007;Ritchie & Thomas 2001;Okamoto et al. 2003). To overcome the difficulties in treating the shock properly, many efforts have been made to improve the prescription and/or to invent a new scheme of SPH Price 2008;Read et al. 2010;Hosono et al. 2013;Hopkins 2013;Yamamoto et al. 2015). We compared the results using different SPH schemes, i.e., the conventional SPH and DISPH, under the same configuration. Figure 14 compares the conventional (standard) SPH (SSPH) and the DISPH with the same initial condition. The amount of unbound mass due to the stripping is less in the conventional SPH model than in the DISPH model. This can be interpreted as a consequence of introducing an unphysical surface tension at the contact surface between the star and the ejecta, which suppresses the stripping of the SPH particles.