The population synthesis of Wolf-Rayet stars involving binary merger channels

Wolf-Rayet stars (WRs) are very important massive stars. However, their origin and the observed binary fraction within the entire WR population are still debated. We investigate some possible merger channels for the formation of WRs, including main sequence (MS)/ Hertzsprung Gap (HG) + MS, He + HG/ Giant Branch (GB). We find that many products produced via binary merger can evolve into WRs, the MS/ HG + MS merger channel can explain WRs with luminosities higher than $\sim 10^{5.4}$\,L$_{\odot}$, while the He + HG/ GB merger channel can explain low-luminosity WRs in the range of $10^{4.7}$\,L$_{\odot}$\,$\sim$\,$10^{5.5}$\,L$_{\odot}$. In the population synthesis analysis of WRs, we assume an initial binary fraction ($f_{\rm ini,bin}$) of 50\% and 100\% for massive stars. We also assume that MS/ HG + MS merger products are non-rotating or rapidly rotating ($\omega/\omega_{\rm crit}=0.8$). In different cases, the calculated single fractions of WRs range from $22.2\%$ to $60.6\%$ in the Milky Way (MW) and from $8.3\%$ to $70.9\%$ in the Large Magellanic Cloud (LMC). The current observations fall within the range of our calculations. When the merger product of MS/HG+MS rotates rapidly, we estimate that there are approximately 1015 to 1396 WRs in the MW and 128 to 204 WRs in the LMC. Our model also roughly reproduces the observed single-peak luminosity distribution of WRs in the MW. However, the weak bimodal luminosity distribution observed in the LMC is not reproduced in our model. We assess that this may be due to the model underestimating the mass-loss rate in the LMC. In conclusion, we consider that the binary merger is significant formation channel for WR formation, and can explain the observed high fraction of the single WRs in the total population.


INTRODUCTION
Wolf-Rayet stars (WRs) are typically massive heliumburning stars that have lost their hydrogen envelopes.It is characterized by broad and intense emission lines (Crowther 2007).The strong stellar winds of WRs contribute to the nucleosynthesis of elements and drive the chemical evolution of galaxies (Heger et al. 2003;Langer 2012).Moreover, the study of WRs is crucial for gaining insights into various important objects, such as X-ray binaries and black holes (van den Heuvel et al. 2017).
The origin of WRs is still poorly understood, despite their significance and the progress made in studying them.Some calculations using single-star evolutionary models can ex-Corresponding author: Chunhua Zhu, Guoliang Lü chunhuazhu@sina.cn,guolianglv@xao.ac.cnZhuoWenli2024@163.com plain a significant portion of the properties of highly luminous WRs through high mass-loss rate caused by strong stellar winds and/ or rotation induced mixing (Langer et al. 1994;Meynet & Maeder 2005;Cui et al. 2018).However, the mass-loss rate ( Ṁloss ) of massive stars remains highly uncertain (Mauron & Josselin 2011).Additionally, observations indicate that most massive early O/B-type stars do not exhibit very high rotation (Ramírez-Agudelo et al. 2013, 2015, 2017).
Theoretically, binary evolutionary models provide a more robust explanation for the formation of WRs.In these models, the donor evolves into WRs by losing most of its H-rich envelope during the mass transfer (MT) process (Shenar et al. 2019).Pauli et al. (2022) also speculate that 80% of the WR population in the Large Magellanic Cloud (LMC) consists of binary systems.However, observationally, many WRs do not have identified companions, and the fraction of single WRs to the entire WR population is estimated to be over 60% in the Milky Way (MW) (Vanbeveren & Conti 1980;Vanbeveren et al. 1998;van der Hucht 2001;Dsilva et al. 2020Dsilva et al. , 2022)), and around 70% for the LMC (Bartzakos et al. 2001;Foellmi et al. 2003a,b;Schnurr et al. 2008;Neugent & Massey 2014;Hainich et al. 2014;Shenar et al. 2019).Additionally, the optically thick wind effect in H-free WRs leads to observed effective temperatures that are much lower than model calculations, posing challenges for theoretical models in explaining H-free WRs (Gräfener et al. 2012;Lu et al. 2023).
It is well known that more than 50% of massive stars are born in binary systems (Sana et al. 2013).Additionally, over 30% of the massive stars that exchange mass with a companion will merge (Sana et al. 2012).Therefore, it is interesting to investigate whether massive binary mergers are the main formation scenario for WRs and what properties these WRs possess.In this paper, we study the WR population resulting from binary merger channels.Our paper is structured as follows: in Section 2, we describe the modeling of different merging channels, the evolutionary trajectories, and the properties of the resulting WRs.In Section 3, we simulate the WR population using a synthetic population method.Finally, in Section 4, we present our conclusions.

WRS PRODUCED VIA BINARY MERGER
In this section, we provide details of the modelling of the different merging channels and compare the evolutionary trajectories of the merger products from these channels with the observed positions of WRs in the Hertzsprung-Russell diagram (HRD).

Model
Binary mergers are complex and challenging to simulate accurately.To simplify the process, we utilize the rapid binary code of Massive Objects in Binary Stellar Evolution (MOBSE) (Hurley et al. 2000(Hurley et al. , 2002;;Giacobbo et al. 2018) to calculate the binary evolution leading up to the merger event.During this evolution, the dynamical MT between the binary components can result in a common envelope evolution (CEE), ultimately leading to the formation of a close binary or direct merger into a new star (eg., Han et al. 2020).
In the MOBSE code, the evolution and structure of stars are described by a series of fitting formulae, which are based on calculations using stellar models from Pols et al. (1998) (Hurley et al. 2000(Hurley et al. , 2002;;Kiel & Hurley 2006;Giacobbo et al. 2018).The outcomes of binary system interactions are mainly determined by parameters for common envelope evolution (CEE) and the mass transfer (MT) model.During the CE phase, MOBSE uses a standard energy prescription for treatment (Hurley et al. 2002;Ivanova et al. 2013;Giacobbo et al. 2018).The efficiency for the orbital energy used to eject the envelope is set to α CE = 1, while the envelope structure parameter is λ = 0.1 (Xu & Li 2010).Additionally, the MT model and any unmentioned parameter settings align with the default model (values) outlined in Giacobbo et al. (2018).
Based on the studies by Hurley et al. (2002) and Giacobbo et al. (2018), we consider the following merger scenarios for the formation and evolution of WRs: merger between a main sequence (MS)/Hertzsprung Gap (HG) star and a MS star, and merger between a helium star (He) and a HG/ Giant Branch (GB).The structure and evolution of the newly formed star resulting from the binary merger are simulated using the Modules for Experiments in Stellar Astrophysics (MESA) code (version 10398) (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018)).The subsequent subsections provide detailed explanations of these merger channels.In addition, we use MESA to cross-check the MOBSE calculations.We try to keep the two codes as consistent as possible.If there is an inconsistency between the results of MOBSE and MESA calculations, we use the results calculated by MESA.In MESA, based on Stancliffe & Eldridge (2009), Shao & Li (2021), Pauli et al. (2022), andDutta &Klencki (2023), for binary systems with very short initial orbital periods, the donor undergoes very efficient mass stripping while the gainer rapidly expands and reaches critical rotation because it accretes too much material within a short timescale.The gainer that reaches critical rotation is usually considered not to continue accreting material, and excess material is lost from the vicinity of the gainer (Stancliffe & Eldridge 2009;Shao & Li 2021).In our simulations, we find that most binary systems with initial orbit periods of about ≤ 4 days, the gainer expands significantly and reaches critical rotation during MT.Only a few binary systems have the gainer rotation speeds that may stay below the critical value due to tidal synchronization (Ritter 1988).Following van den Heuvel (1994) and He et al. (2024), we use a parameter β = 0.75 (β is the fraction of mass lost from the vicinity of the gainer) for this case.For binary systems with initial orbital periods of about > 4 days, we find that the mass stripping efficiency of the donor is relatively weak, and conservative MT can achieve MT occurring in case A/B.

MS/HG + MS merger model
In the case of massive binaries with small orbital periods and/ or large eccentricity, the merger occurs when the two parent stars are brought into deeper contact.This is caused by the continued expansion of the primary as it evolves into the MS/ HG phase, as well as the sharp orbital contraction resulting from the loss of angular momentum (Podsiadlowski et al. 1992;Wellstein et al. 2001;de Mink et al. 2014).The merger product has a higher mass and luminosity compared to the two parent stars (the primary and secondary before the merger) (Glebbeek et al. 2013), which in turn drives stronger stellar winds (Vink et al. 2001).This implies the possibility of the formation of WRs in the subsequent evolution of these merger products.
In the MOBSE code, the two stars in the binary system will merge if the condition (R 1 +R 2 ) > a(1−e) (Hurley et al. 2002;Eldridge et al. 2017;Giacobbo et al. 2018;Riley et al. 2022) is met during the expansion of the MS/ HG phase.Here, a represents the binary separation, e is the eccentricity, and R 1 and R 2 are the radii of the two components.Following the approach of Tout et al. (1997) and Hurley et al. (2002), the material of the two merging components is assumed to be completely mixed.The hydrogen in the core of the merger product will be replenished, resulting in the rejuvenation of the merged star to approximately a zero-age main sequence (ZAMS) star (Hellings 1983;Vanbeveren et al. 1998;Menon et al. 2021).Due to the relatively large binding energy of the MS/ HG stars (Taam & Sandquist 2000;Ivanova & Taam 2004;Dominik et al. 2012), no mass loss is assumed during the collision.In other words, the mass of the merger product is equal to the sum of the masses of the two parent stars, i.e., The post-merge evolution of the merger product is performed using the MESA code.For the simulations, we adopt typical metallicities of 0.014 for the MW (Ekström et al. 2012) and 0.006 for the LMC (Eggenberger et al. 2021), as WRs are mainly found in the MW and the LMC.In the MESA simulations, convection is treated using the Ledoux criterion and a mixing length parameter of 1.5 (Langer 1991).A semi-convective parameter of 1 is also used (Langer 1991).The overshoot area is calculated using a step-function approach, with an overshooting parameter set to 0.335 (Brott et al. 2011).Additionally, a thermohaline mixing parameter of 1 is employed (Kippenhahn et al. 1980).To avoid computational convergence failure of the model near the Eddington limit, we enable the option MLT++ (Paxton et al. 2013).
Mass loss is calculated with reference to Brott et al. (2011) andPauli et al. (2022).When the temperature is below the bistability jump temperature, the maximum value of the massloss rate from Vink et al. (2001) or Nieuwenhuijzen & de Jager (1990) is used.When the temperature is higher than the bi-stability jump temperature, the mass-loss rate for MS stars is calculated using the prescription of Vink et al. (2001).For the surface H abundance (X H ) is between 10 −5 and 0.4, the mass-loss rate from Nugis & Lamers (2000) is used.For the H-free WR stage (X H < 10 −5 ), the mass-loss rate from Yoon (2017) is employed.For X H is between 0.4 and 0.7, the wind mass-loss rate is linearly interpolated between the prescriptions of Vink et al. (2001) and Nugis & Lamers (2000).Since the work of Fullerton et al. (2006), Puls et al. (2006), Krtička & Kubát (2017), Beasor et al. (2020), andBjörklund et al. (2021) suggests that the mass-loss rate of massive stars at different stages is not as high as expected, we adopt a similar approach to Yoon et al. (2006) and Brott et al. (2011), reducing the mass-loss rate from Nugis & Lamers (2000) by a factor of 10.It is worth noting that we use the same stellar wind scheme described above for both the MW and LMC.Considering that the MW has a higher mass-loss rate, while the LMC has a lower one, our stellar wind model may be reasonable for the LMC, while it may be underestimated for the MW.
In the post-merge evolution, the effects of angular momentum conservation, off-center collisions (Costa et al. 2022), and the formation of an excretion disk around the merged product (Fryer & Heger 2005) are taken into consideration.These complexities can result in the merged product either rotating or not rotating.To quantify the effect of rotation on the merger product, a ratio of the rotational angular velocity to the critical rotational angular velocity, ω/ω crit , is used.In the simulations, values of ω/ω crit = 0 and 0.8 are adopted.The efficiency parameter for rotational mixing is set to 1/30 (Brott et al. 2011), and the inhibition parameter for the chemical gradient is set to 0.1 (Brott et al. 2011).We also considered various instabilities caused by rotation, such as dynamical shear instability, secular shear instability, Eddington-Sweet circulation, and the Goldreich-Schubert-Fricke instability (Heger et al. 2000).Furthermore, the massloss rate of a rotating massive star is increased according to the equation Ṁ = , where β = 0.43 (Langer 1998).

He + HG/ GB merger model
The binary systems including a He and a HG or GB star generally have undergone a stable MT or a CEE.Its companion can fill its Roche lobe at a subsequent evolutionary stage (HG/GB).In MOBSE, the He star and the HG or GB star experience a CEE if q = M donor M gainer > q crit .For a HG star as the donor: q crit,HG = 4 (Tout et al. 1997); for a GB star as the donor: (Hjellming & Webbink 1987), where M 2,c is the core mass of donor.Based on a standard CEE model, He star and its companion can merge into a new star if E bind > α CE ∆E orb , where ∆E orb is the change of orbital energy of the binary from CEE begin to end, E bind is the binding energy of the donor envelope, and α CE is the efficiency for the orbital energy used to expel the envelope (Hurley et al. 2002;Ivanova et al. 2013).Based on Eggleton (1983), Han et al. (1995), Tout et al. (1997), Xu & Li (2010) and Farrell et al. (2020), the mass lost (∆m) of the donor due to the orbital energy changes (∆E orb ) during CEE can be approximately calculated by where M 2,env is the envelope mass of the donor, u is the specific internal energy of the gas, the thermal ejection efficiency parameter α th is 1, α CE is 1, the core boundary is defined as the central Hpoor region of X H < 10 −4 .
For the structure of the merged product, following Lombardi et al. (2002) and Gaburov et al. (2008), the He star with higher density (lower entropy) sinks into the core of the donor as the core of the merged product.According to Zhang & Jeffery (2013) and Zhang et al. (2020), the merging process is approximated by He stars accreting the average mass fraction of the donor core initially, followed by accreting the average mass fraction of the remaining envelope (M env −∆m), with an accretion rate ( Ṁacc ) of 1 M ⊙ /yr.That is, the new born star has a core mass combining He star and the donor core, and a H-rich envelope mass of M env − ∆m.

Evolutionary Tracks of WRs produced via binary merger
Our study focuses on the evolution of the merger products and explores whether these products from different merger channels can evolve into WRs.We aim to determine which subtypes of WRs can be explained by these different merger channels.
Based on previous studies by Smith & Maeder (1991), Meynet & Maeder (2003) and Georgy et al. (2012), the different stages or subtypes of stars in our evolutionary model are defined as follows: He stars (log(T eff /K) > 4.0 and X H < 0.3), WRs (log(T eff /K) > 4.0, X H < 0.3, and the mass-loss rate > 10 −6 M ⊙ /yr) (Crowther 2007).H-rich WNs (10 −5 < X H < 0.3), H-free WNs (X H < 10 −5 and X N > X C ), and WCs (X H < 10 −5 and X N < X C ). Aguilera-Dena et al. ( 2022) point out that WOs may be special WCs, so possible WOs phase is not considered in the model.When considering the observed sample of WNs, we follow the distinction made by Hainich et al. (2014) and Pauli et al. (2022), where WNs with X H > 0 are classified as H-rich WNs and those with X H = 0 are classified as H-free WNs.In addition, there are some H-rich WNs with X H > 0.3 in the observed WRs of Hainich et al. (2014) and Hamann et al. (2019).For example, X H of BAT99-49 and BAT99-111, their X H may be higher than 0.5.However, Hainich et al. (2014) and Pauli et al. (2022) considered that these stars are unlikely to be WRs but are instead more likely to be central hydrogen-burning stars.Considering that the X H values of these samples exceed the range designated for WRs in the studies by Smith & Maeder (1991), Meynet & Maeder (2003), Georgy et al. (2012), and Yusof et al. (2013), we exclude these observed samples from our analysis.Furthermore, some WCs observed in the LMC lack information about their effective temperature.Consequently, this subset of the sample is excluded from our HRD analysis.

Evolution of MS/ HG + MS mergers
For the MS/MS+HG merger product, based on Schneider et al. (2016), hydrogen can mix into the center of the merger product during the merger process, serving as fuel.The core mass and luminosity of the merger product are larger and higher than those of its two parent stars (the primary and secondary before the merger).In our model, we assume that the merger product rejuvenates to become a ZAMS star.In our model, we assume that the merger product is rejuvenated to a ZAMS star.
Fig. 1 presents the evolutionary stage of the merger product as a function of its initial mass and the normalized time during the post-MS evolution.Rotation lowers the lower limit of the initial mass of the merger product for the formation of H-rich WNs, H-free WNs, and WCs, attributed to the increased mass-loss rate of merger products due to rotation (Langer 1998).In the MW, without rotation, a merger product with mass higher than approximately 30 M ⊙ can evolve into the H-rich WNs phase.However, the lower mass limits in models with a rapid rotation of 0.8 ω crit are around 19 M ⊙ .Without rotation, a merger product with mass of about 84 M ⊙ can enter the H-rich WNs phase at the terminal-age main sequence (TAMS).However, in models with a rapid rotation of 0.8 ω crit , this mass requirement drops to around 34 M ⊙ .In the LMC, without rotation, the merger product needs approximately 55 M ⊙ to evolve into the H-rich WNs phase.With a rapid rotational velocity of 0.8 ω crit , this lower limit decreases to about 24 M ⊙ due to rotation induced chemically homogeneous evolution (CHE).CHE enables the merger product to enter the H-rich WNs phase at TAMS with mass greater than about 30 M ⊙ .Our results are consistent with Pauli et al. (2022), the H-free WNs phase appears to act as a transition phase between the H-rich WNs and WCs phase, and it typically has the shortest timescale throughout the WRs phase.
Fig. 2 illustrates the evolutionary tracks of the merger products and compares them with the observed WRs in the HRD.The formation of larger mass products through MS/HG + MS mergers can result in stronger stellar winds, leading to the formation of WRs.When compared with observations, we find that the merger products can approximately explain only the high-luminosity WRs.For the H-rich WNs, the lower limit of the luminosity for ω/ω crit = 0 (and ω/ω crit = 0.8, respectively) is approximately over 10 5.6 L ⊙ (10 5.4 L ⊙ for the MW and 10 6.1 L ⊙ (10 5.5 L ⊙ ) for the LMC.For H-free WNs, only luminosities higher than about 10 5.8 L ⊙ (10 5.5 L ⊙ ) for the MW and 10 6.1 L ⊙ (10 5.6 L ⊙ ) for the LMC can be explained.For WCs, only luminosities higher than about 10 5.7 L ⊙ (10 5.4 L ⊙ ) for the MW and 10 6.4 L ⊙ (10 5.7 L ⊙ ) for the LMC.In the LMC, models without rotation have a higher lower limit for H-rich WNs luminosity than those in Pauli et al. (2022), which is 10 6.0 L ⊙ .The main reason is that our models have a lower mass-loss rate.In addition, since the option MLT++ eliminates the effect of envelope inflation for massive stars near the Eddington limit, it leads to a much larger effective temperature for our H-free WRs model than most observed H-free WRs.

Evolution of He + HG/ GB mergers
During the He + HG/ GB merger, the ejection of the Hrich envelope outside the system suggests that these mergers do not necessarily require high stellar wind to become WRs.However, not all He + HG/ GB mergers can form WRs. If ∆E orb during the CEE is too large or too small, it can result in either successful ejection of the envelope or the retention of a significant amount of the H-rich envelope.
Fig. 3 shows a typical example of He + HG/ GB merger.Initially, the components with Z=0.014 are on the ZAMS, the primary's and the secondary's masses are 11 M ⊙ and 9.5 M ⊙ respectively, the initial orbital period is 2 d and the MT is highly non-conservative, with a MT efficiency of 1−β = 0.25.After two stable MT, the primary evolves into a 2.3 M ⊙ He star, and the secondary still stays in MS phase but its mass increases to 11.4 M ⊙ .Approximately 6.8 M ⊙ of mass is lost outside the binary system (referred to as isotropic reemission).When the secondary evolves into HG, it fills its Roche lobe and CEE occurs.At this time, M 1 =2.2 M ⊙ , M 2 =11.3 M ⊙ , the binary separation is 110 R ⊙ , which results in about 6.4 M ⊙ mass of the envelope of the secondary being lost from the system during CEE.It means that the He star and its companion merger into a new star with mass of 7.1 M ⊙ .
However, the initial parameter space of the He + HG/ GB merger channel is very small, specifically, the initial primary's masses is 11 ∼ 16 M ⊙ , and the initial mass ratio q ini is 0.85 ∼ 0.95.This is because He stars have a very short lifetime, and most of them cannot survive until the HG/ GB stars undergo reverse MT.Thus, the mass of the merger prod-uct for this channel is roughly 7 ∼ 15 M ⊙ .Fig. 4 shows the evolutionary tracks of the products formed via He + HG/ GB merger and compares them with the observed WRs in the HRD.After the merger, the products enter the central He-burning phase, and they evolve towards the bluer region of the HRD.Compared with observed WRs, the He + HG/ GB merger products can explain H-rich WNs with luminosities between 10 4.7 L ⊙ ∼ 10 5.4 L ⊙ for the MW and between 10 5.0 L ⊙ ∼ 10 5.5 L ⊙ for the LMC.Due to the low mass of these merger products, the stellar wind is not strong enough to completely remove the H-rich envelope, preventing these merger products from evolving into the H-free WRs phase.

WR POPULATION
Due to most WRs are affected by optically thick stellar winds, the observed effective temperature and radius of WRs are not always reliable.Therefore, comparison with the observed position of WRs in the HRD alone are often insufficient.
To address this limitation, our study utilizes population synthesis to simulate the entire population of WRs.Through this approach, we investigate both single and binary WRs systems, enabling a comprehensive analysis of their characteristics.

Population synthesis method
Considering that more than 50% ∼ 90% of massive stars are in binary or multiple systems (Sana et al. 2012(Sana et al. , 2014;  2019).Observations of WCs for MW were obtained from Sander et al. (2019).Observations of H-rich WNs and H-free WNs for LMC were obtained from Hainich et al. (2014).Observations of WCs for LMC were obtained from Crowther et al. (2002).It should be noted that we do not consider observations with X H > 0.3.
In order to investigate the binary fraction of WRs, Dsilva et al. ( 2023) consider that a WRs is single if its radial velocity change (∆RV) satisfies ∆RV < 50 km/s (the observational inclination is assumed to be 90 • ), or else this WRs is considered as WRs binary.In our simulations, there are some WR binaries with very wide orbits.Their ∆RV should be smaller than 50 km/s, and these WRs can be treated as single WRs on the observations.Following Dsilva et al. (2023), we assume that WRs are single if these WRs in binary systems at perihelion have ∆RV < 50 km/s (In our grid, the orbital period is approximately greater than 10 3.0 ∼10 3.5 days), otherwise these binaries are called as WRs binaries.
To distinguish them from single WRs mentioned above, WRs formed through binary mergers are referred to as merger WRs.Therefore, in this work, based on the formation scenario, isolated WRs are divided into single and merger WRs.
Using MOBSE, we evolve 10 7 massive binary systems, and obtain a large grids for WRs candidates via single star and binary merger channels.Considering that MOBSE does not provide detailed stellar structure, we use MESA to evolve, cross-check, and linearly interpolate connections within the grids of different channels in the MOBSE calculations at a certain step size.For the single-star channel and the MS/ HG + MS merger channel, the mass range (initial primary or initial isolated single-star or the merger product) is roughly 20 M ⊙ ∼ 100 M ⊙ .For He + HG/ GB merger channel, the model grid covers roughly the initial primary mass in the range of 11 M ⊙ ∼ 16 M ⊙ , with the mass ratio in the range of 0.85 ∼ 0.95, and the initial orbital period is in the range of 1.5 ∼ 2.3 days.In addition, for the single star channel, we assume that all models are non-rotating.For the MS/ HG +MS merger products both ω = 0 and ω = 0.8ω crit are assumed.For the He + HG/ GB merger channel, we assume that the angular momentum is dissipated through the CEE phase, and thus the merger product of He+HG/GB is non-rotating.The above channels contain a total of 206 MESA models.Simultaneously, some massive binary systems experience MT and can eventually evolve into WRs + O-type binary systems.We also use MESA to evolve, cross-check, and linearly interpolate connections within the grids of binary nonmerger channels in the MOBSE calculations at a certain step size.In MESA, we refer to Dutta & Klencki (2023) and treat the secondary as a point mass, mainly focusing on the structural information of the primary forming WRs.We adopt the scheme from Ritter (1988) to calculate the MT.The accretion efficiency is described in Section 2.1.In this grid, the initial mass range for the primary is from 20 M ⊙ ∼ 100 M ⊙ .The range for the initial mass ratio is from 0.3 ∼ 0.95.The initial orbital periods range from 5 ∼ 1000 days.Through a three-dimensional interpolation of primary mass, mass ratio and orbital period, we can obtain a detail structure and evolution of primary.The binary channel grid in this work includes 2214 MESA models.Fig. 5 shows some of the possible results of the formation of WRs by binary channel.

Properties of WR population
Combining MOBSE and MESA, we obtain WR population via the method of population synthesis.This method is used in many literatures (Lü et al. 2012(Lü et al. , 2017;;Yu et al. 2019;Han et al. 2020;Zhu et al. 2021;Guo et al. 2022;Zhu et al. 2023).

WRs expected and the single/ binary fraction of WRs
Estimating WRs' number is crucial for understanding the overall population of WRs and for conducting effective searches.Following the method in Dutta & Klencki (2023), we calculate the WRs' numbers in the MW and the LMC.We assume that both the MW and LMC have a constant star formation rate (SFR).Based on Chomiuk & Povich (2011) and Harris & Zaritsky (2009), SFRs are 1.9 M ⊙ yr −1 in the MW and 0.4 M ⊙ yr −1 in the LMC, respectively.The WRs' number can be approximately estimated via where ⟨M⟩ is the average mass of initial mass function; f WR is the fraction of systems evolving into WRs; ⟨t tot ⟩ is the average lifetime of such systems; ⟨t WR ⟩ is the average lifetime during the WRs phase.Moreover, WRs are predominantly found as single stars in the observations, whereas most of them were expected to be in binaries according to previous theoretical models.Involving binary merger, we find that single WRs may dominate in the WR population.
Fig. 6 shows the calculated number of WRs and the fractions of single/ binary WRs for different cases, in which single WRs include WRs in binary systems with an orbital period longer than about 1000 days (∆RV < 50 km/s), evolved from massive stars born in isolated systems (when f ini,bin = 50%) and isolated WRs produced via binary merger.For the merger scenario our models predict 704 ∼ 1396 WRs and single fraction 22.2% ∼ 60.6% for the MW, 47 ∼ 204 WRs and single fraction 8.3% ∼ 70.9% for the LMC.In the MW, the number of WRs has been estimated, ranging about 6500 (van der Hucht 2001) and 1200 (Maeder & Lequeux 1982;Crowther 2015).Our results align with the previous estimates.So far, the known WRs in the MW and LMC are about 642 (Crowther 2015) and 154 (Neugent et al. 2018), respectively.It means that there is still a significant fraction of unknown WRs in the MW.This may be attributed to the Galactic disk and a large amount of dust obscuring many WRs.Furthermore, our simulations give upper and lower limits on the true number of WRs for the MW and LMC.Clearly, the actual number of WRs depends on f ini,bin , the rotation of MS/ HG + MS merger products, and the SFR.Considering the observation precision, the calculated single fraction of WRs is very consistent with the current observation value for the case of a high rotation rate (ω = 0.8 ω crit ) of the MS/HG + MS merger product.
On the other hand, the decrease in f ini,bin and the increase in rotation of the MS/ HG + MS merger products always favours an increase in the single fraction of WRs.This is because the initial increase in the number of isolated massive stars and the rapid rotation increase the possibility of WRs formation by enhancing the mass-loss rate and extending the lifetime of WRs through enlarging the mass of the burning core (Sibony et al. 2023).In the LMC, the very small contribution of isolated massive stars to the WR population, and the rotation induced CHE are key to the formation of WRs.Additionally, a decrease in f ini,bin will always lead to a significant decrease in the calculated number of WRs.This is because the lifetimes of WRs formed through the single star channel are much shorter than those formed through binary channels.
Without considering mergers, the population of WRs obtained is mostly dominated by binary stars, and the calculated number of WRs in the MW is much lower than the estimates by Maeder & Lequeux (1982), Crowther (2015) and van der Hucht (2001) and also significantly lower than the observed number of WRs in the LMC.It is worth pointing out that in the LMC for f ini,bin = 50% our calculated binary fraction of WRs is essentially in agreement with Pauli et al. (2022).

Luminosity distribution of WRs
It is important to explore the luminosity distribution of WRs.Figs. 7 and 8 show the luminosity distributions of different subclasses of WRs in the MW and LMC predicted to be obtained after normalisation of the observed samples.
In Fig. 7, within the MW, the luminosities of the observed WRs from different subclasses exhibit a roughly unimodal distribution.Among all subclasses, a decrease in f ini,bin significantly increases the contribution from the single-star channel, thereby augmenting the number of WRs in the highluminosity region (> 10 5.8 L ⊙ ).Rapidly rotating MS/HG + MS merger products form fainter WRs than in non-rotating scenarios due to an increased mass-loss rate, which leads to a shift of the predicted peak of the distribution toward lower luminosities.Upon comparison, we assess that the scenario where merger products are rapidly rotating roughly aligns with the observed distribution of H-rich WNs and H-free WNs.However, in the case of WCs, even with rapidly rotating MS/ HG + MS merger products, the number of lowluminosity WCs (< 10 5.4 L ⊙ ) in our simulation is still much lower than the observed estimate.We evaluate that appropriately enhancing the mixing efficiency of the merger products can facilitate the formation of low-luminosity WCs.
In Fig. 8, the luminosities of the observed WRs from different subclasses in the LMC generally exhibit a slightly bimodal distribution.Unlike in the MW, f ini,bin has a weaker influence on the luminosity distribution in the LMC, which is attributed to the significantly lower contribution of the  and 4, but for WRs which are formed by the primaries through the binary channel.In all the subgraphs, the evolution trajectory from bottom to top corresponds to initial masses of the primary of 20 M ⊙ , 34 M ⊙ , 50 M ⊙ , 65 M ⊙ , 80 M ⊙ and 100 M ⊙ respectively.Moreover, for all these evolutionary trajectories, the initial mass ratio is consistently set at q ini = 0.6.single-star channel to WRs in the LMC.Among all subclasses, rapidly rotating MS/ HG + MS merger products shift the predicted peak further to lower luminosities due to CHE.Considering the observed high fraction of single WRs in the LMC, we only compare scenarios where MS/HG + MS merger products are rapidly rotating.For H-rich WNs, H-free WNs, and WCs, the luminosity distributions simulated in our models roughly match the first peak of the observed distributions.However, the second peak observed in the samples (∼ 10 6.4 L ⊙ for H-rich WNs and H-free WNs, and ∼ 10 6.8 L ⊙ for WCs) does not appear in our simulations.We evaluate that appropriately increasing the mass-loss rate in the single-star channel could potentially improve this discrepancy, as the single-star channel's contribution to WRs luminosities primarily fall between 10 6.0 L ⊙ and 10 6.6 L ⊙ .

CONCLUSIONS
Considering binary interaction and using MOBSE and MESA code, we investigate the WRs formation via massive star merger.We find that the massive-star merger can successfully produce WRs.The MS/ HG + MS merger, producing a massive star with high rotation and high mass-loss rate, can evolve into the WRs with luminosity higher than 10 5.4 L ⊙ .The He + HG/ GB merger, ejecting most of the H-rich envelopes, can produce the WRs with low luminosity WRs of 10 4.7 L ⊙ ∼ 10 5.5 L ⊙ .
Combining MOBSE and MESA codes and using the method of population synthesis, we simulate WR population.We estimate that there are approximately 704∼1396 WRs in the MW and 47 ∼ 204 WRs in the LMC.Here, the contribution of WRs produced via binary merger to total WRs is about 12.6% ∼ 48.5% for the MW and 5.5% ∼ 69.1% for the LMC.In addition, the fraction of single WRs in whole WRs is about 51.9% ∼ 60.6% for the MW and 69.8% ∼ 70.9% for the LMC in the models with rapid rotation after binary merger, which is consistent with the observed estimate.These indicate that the binary merger is significant formation scenario for WRs formation, and can explain the high fraction of the single WRs in the total population.

Figure 1 .
Figure 1.MS/HG + MS merger products exhibit different WRs phases during post-MS evolution.These phases, including the pre WRs phase (black), H-rich WNs phase (blue), H-free WNs phase (green), and WCs phase (red), can be observed as a function of the initial mass of the merger products, the normalized time to He burning, and variations in metallicities and rotations.

Figure 2 .
Figure 2. Positions of WRs observed in the MW and LMC, and evolutionary trajectories of MS/HG + MS merger products with different masses.In the MW and LMC, the evolutionary trajectories from bottom to top correspond to initial masses of 20,24,30,40,56,80,100 M ⊙ for the merger products.The different types/stages are marked with colours.Black: pre WRs; blue: H-rich WNs; green: H-free WNs; red: WCs.Observations of H-rich WNs and H-free WNs for MW were obtained from Hamann et al. (2019).Observations of WCs for MW were obtained fromSander et al. (2019).Observations of H-rich WNs and H-free WNs for LMC were obtained fromHainich et al. (2014).Observations of WCs for LMC were obtained fromCrowther et al. (2002).It should be noted that we do not consider observations with X H > 0.3.

Figure 3 .
Figure 3. Evolutionary tracks in HRD for WRs produced via the He + HG/GB merger channel.Initially, the components with Z = 0.014 are on the ZAMS, the primaries and the secondaries masses are 11 M ⊙ and 9.5 M ⊙ respectively, and the initial orbital period is 2 d.Markers of different colors and shapes indicate the beginning and end of the Roche lobe overflow (RLOF), respectively.It is worth noting that the accretion efficiency of MT is 1 − β = 0.25.

Figure 5 .
Figure 5. Similar to Figs. 2and 4, but for WRs which are formed by the primaries through the binary channel.In all the subgraphs, the evolution trajectory from bottom to top corresponds to initial masses of the primary of 20 M ⊙ , 34 M ⊙ , 50 M ⊙ , 65 M ⊙ , 80 M ⊙ and 100 M ⊙ respectively.Moreover, for all these evolutionary trajectories, the initial mass ratio is consistently set at q ini = 0.6.

Figure 6 .
Figure6.The fractions of single and binary WRs, as well as the total number of WRs calculated in different scenarios, are determined through population synthesis.We assume different rotations for the MS/HG + MS merger product, either ω = 0 or ω = 0.8ω crit , and consider initial binary fractions of f ini,bin = 50% and 100% for massive stars.The single channel pertains to WRs in binary systems with orbital periods exceeding approximately 10 3.0 ∼ 10 3.5 days (∆RV < 50 km/s) and WRs evolving from massive stars in isolated systems.The merger channel relates to single WRs resulting from binary mergers.The binary channel involves WRs binaries with orbital periods shorter than about 10 3.0 ∼ 10 3.5 days (∆RV > 50 km/s).These channels are represented by distinct stripes.Observations of the single and binary fractions of WRs in the MW and LMC come fromDsilva et al. (2022)  andShenar et al. (2019), which are shown as dashed lines.The total number of currently observed WRs in the MW and LMC is sourced fromCrowther (2015),Massey et al. (2014), andNeugent et al. (2018).

Figure 8 .
Figure 8. Similar to Fig. 7, but for WRs in the LMC.The luminosity distribution of different subclasses of WRs obtained after normalisation with the observed sample.Observations of H-rich WNs and H-free WNs in the LMC were obtained from Hainich et al. (2014).Observations of WCs in the LMC were obtained from Bartzakos et al. (2001).Additionally, four of the 26 WCs from Bartzakos et al. (2001) were excluded based on the work of Neugent et al. (2018).