{'fraction_non_alphanumeric': 0.033401251637316255, 'fraction_numerical': 0.02672463979042352, 'mean_word_length': 4.51676033721397, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Nucleosynthetic Yields from Neutron Stars Accreting in Binary Common Envelopes

Massive-star binaries can undergo a phase where one of the two stars expands during its advanced evolutionary stage as a giant and envelops its companion, ejecting the hydrogen envelope and tightening its orbit. Such a common envelope phase is required to tighten the binary orbit in the formation of many of the observed X-ray binaries and merging compact binary systems. In the formation scenario for neutron star binaries, the system might pass through a phase where a neutron star spirals into the envelope of its giant star companion. These phases lead to mass accretion onto the neutron star. Accretion onto these common-envelope-phase neutron stars can eject matter that has undergone burning near to the neutron star surface. This paper presents nucleosynthetic yields of this ejected matter, using population synthesis models to study the importance of these nucleosynthetic yields in a galactic chemical evolution context. Depending on the extreme conditions in temperature and density found in the accreted material, both proton-rich and neutron-rich nucleosynthesis can be obtained, with efficient production of neutron rich isotopes of low Z material at the most extreme conditions, and proton rich isotopes, again at low Z, in lower density models. Final yields are found to be extremely sensitive to the physical modeling of the accretion phase. We show that neutron stars accreting in binary common envelopes might be a new relevant site for galactic chemical evolution, and therefore more comprehensive studies are needed to better constrain nucleosynthesis in these objects.


INTRODUCTION
A number of dedicated observing programs have shown that most massive stars are formed in systems with at least one binary companion (Kobulnicky & Fryer 2007;Kiminki et al. (CE) phase where the expanding star envelops its companion, causing the core of the expanding star and its companion to share a common envelope. The CE phase causes the binary to tighten its orbit and is postulated to explain many of the tight-orbit, massive-star binaries (Ivanova et al. 2013).
In the formation of a variety of massive star binaries including X-ray binaries and double compact object systems, the stellar system evolves through one, and often two, CE phases. In the first CE phase, the more massive star (primary) evolves off the main sequence, enveloping its companion. In some cases, the resultant tightening of the orbit produces a binary that is sufficiently close that, even after the subsequent collapse and explosion of the primary, the wind of the companion can accrete onto the neutron star. This is a common scenario behind the production of massive X-ray binaries. In some cases, when the companion star evolves off the main sequence, a second CE phase where a neutron star is enveloped by the companion occurs. This can tighten the orbit prior to the supernova explosion of the companion that produces binary pulsar systems and compact binaries that are believed to be the site of short-duration gamma-ray bursts (GRBs; Fryer et al. 1999b;Bloom et al. 1999), including merging neutron star systems (Dominik et al. 2012a(Dominik et al. , 2013(Dominik et al. , 2015 like the one recently detected by advanced LIGO/Virgo (Abbott et al. 2017). If the neutron star (NS) spirals into the core of the companion, it can produce a long-duration GRB, the so-called helium-merger model (Fryer & Woosley 1998).
The ultimate fate of the binary in this CE phase depends on the masses of the stars and the orbital separation at the onset of the phase. Many systems eject the hydrogen envelope, forming a binary consisting of a helium star and a NS. Others do not have sufficient orbital energy to eject the hydrogen envelope prior to merging with the helium core. These helium mergers were initially proposed to be a long-lived giant star powered by a central Eddingtonrate accreting neutron star known as a Thorne-Zytkow object (Thorne & Zytkow 1975). However, calculations including neutrino processes found that most of the energy released in the neutron star accretion would be radiated efficiently through neutrinos, allowing the neutron star to accrete at the Bondi-Hoyle rate, causing it to rapidly collapse to a black hole (Fryer et al. 1996;Fryer & Woosley 1998).
This helium-merger system, forming a black hole accreting system became one of the proposed black hole accretion disk (BHAD) GRB models (Popham et al. 1999;Zhang & Fryer 2001). Subsequent simulations have studied the potential of this system to produce ultra-long duration gammaray bursts or peculiar supernovae (Fryer et al. 2013;Soker & Gilkis 2018).
Material accreting onto neutron stars is not completely incorporated into the neutron star. If the material has enough angular momentum, it can form a disk that could ultimately drive a jet. This is believed to be rare in most CE scenarios (Murguia-Berthier et al. 2017). For neutron star systems, even if the material does not have a sufficient angular momentum to form a disk, some of the accreting material will be reheated and ejected (Fryer et al. 2006;Fryer 2009). During this accretion process, temperatures and densities become so high that both neutrino emission (that can alter the electron fraction) and nuclear burning can significantly alter the composition of the material. For the high accretion rates of supernova fallback, the reheated ejecta can burn into heavy r-process elements (Fryer et al. 2006). Fallback accretion rates range from a few times 10 −3 to 1 M s −1 . CE accretion rates are typically lower than these rates: ranging from 10 −4 −10 6 M yr −1 (note the former is in per second while the latter is per year).
In this paper, we will study the yields from these lower CE accretion rates. In section 2 we review NS accretion in CE, estimating accretion rates for a range of stellar models at different phases in the star's life. In section 3 we review the range of yields expected as a function of accretion rate from our single zone models. To determine the effect CE yields have on galactic chemical evolution, we must calculate the distribution of binaries and CE scenarios. By using these distributions and stellar models, we can estimate the accretion rates. In section 4.1, we use binary population systems to study yields from stars and stellar populations. We conclude in section 5 with a discussion of the role these yields play in galactic chemical evolution.

Estimating Mass Accretion
When a massive-star companion in a binary with a NS overfills its Roche lobe, its material accretes onto the neutron star. The accretion rate can be much faster than the NS can incorporate, ultimately developing into a CE phase. A number of assumptions are made in estimating this accretion rate and, especially for NSs, there seems to be some confusion on the validity of these assumptions. Here we review the basic physics assumptions and approximations used in this paper to estimate accretion rates. In astrophysics, the standard estimate for accretion onto a point source is the Bondi-Hoyle-Littleton solution (Hoyle & Lyttleton 1941;Bondi 1952). The Bondi radius (R B ) for a neutron star of mass (M NS ) can be determined by the radius that material of velocity v is bound to the NS, i.e.: where G is the gravitational constant. In the simplest case, v is set to the sound speed (c s ). But if the NS is moving with respect to the material, the relative motion (v m ) should also be included. One simple, and often standard, way to include both velocities is through a quadratic sum: v = v 2 m + c 2 s : The accretion rate is roughly the mass within this Bondi radius divided by the free-fall time at this radius. More accurately, this accretion rate is: where ρ is the density of the ambient medium and λ BHL is a non-dimensional parameter: λ BHL = 2/(3π) if we assume freefall. Calculations of Bondi-Hoyle accretion allow refinement of the value for λ BHL and determination of the accuracy of our solution to include the different velocity terms (Ruffert 1994a). In scenarios like our CE phase, there is both a velocity and density gradient across the Bondi radius, and these features drive instabilities in the accretion that can decrease the accretion rate (Ruffert & Arnett 1994;Ruffert 1994b;MacLeod & Ramirez-Ruiz 2014;MacLeod & Ramirez-Ruiz 2015;MacLeod et al. 2017;Murguia-Berthier et al. 2017). Bondi accretion also assumes that matter falling onto the neutron star accretes passively onto the neutron star. However, the gravitational potential energy released as matter accretes onto the neutron star is emitted in radiation and matter outflows that can significantly decrease the accretion rate below the Bondi rate. The Eddington limit is an extreme case of this radiative feedback that assumes all of the energy released is converted into radiation, and the momentum carried by this radiation exerts a force on the inflowing material. This radiation limits the amount of accretion onto an object. The radiative force at radius r is: where L rad is the radiative luminosity, σ is the cross section, and c is the speed of light. Setting this force equal to the gravitational force (F grav = GM NS m/r 2 where m is the mass of the accreting particle), we derive the Eddington luminosity: For accreting, fully ionized hydrogen, m is the proton mass and σ can be set to the Thompson cross-section. Although most studies use these assumptions to calculate the Eddington luminosity, for some scenarios, such as fallback accretion in supernova, the opacity per unit mass can be much higher (Fryer et al. 1999a), lowering the Eddington accretion rate. The Eddington limit on accretion assumes that the radiative luminosity is equal to the gravitational potential energy released: where r NS is the neutron star radius. This limit on the accretion rate assumes spherical symmetry, the radiation is not trapped in the flow and that all the accretion energy is released in radiation. For low accretion rates, many of these assumptions are valid. But as the accretion rate increases, these assumptions lose their validity.
Here we review each assumption individually. Determining whether the radiation is trapped in the flow is difficult without full calculations, but a first order estimate can be made by comparing the diffusive transport velocity (v diff ): where D is the size of the transport region (some fraction of the stellar radius) to the infall velocity, typically set to the free-fall velocity (v ff ): where M enc is the enclosed mass of the star at radius r. Using these approximations, it is found that, except at the beginning of the CE phase when the neutron star is in the outer layers of the hydrogen envelope when M < 10 −4 M yr −1 , the radiation is trapped in the flow (Chevalier 1993;Fryer et al. 1996). Recall that, for our nucleosynthesis models, we are only concerned with accretion rates above 1 M yr −1 where the radiation is truly trapped in the flow (v diffusion << v infall ) and the assumptions needed for the Eddington limit are not applicable.
What about the assumption that all of the energy is emitted in photons? If the radiation is trapped in the flow, the material will shock and settle onto the neutron star. By calculating the post-shock entropy of the material and assuming it piles onto the neutron star, we can derive the temperature and density properties of this accreted material. These estimates find that, unless the entropy is above 10, 000 k B per nucleon, the temperature and density conditions are such that most of the gravitational energy released will be converted into neutrinos and escape the star without impeding the inflow (Chevalier 1989;Fryer et al. 1996). In CE models, the entropy ranges between 10 − 100 k B per nucleon and, to date, no one has constructed a way to avoid rapid neutrino cooling in CE scenarios. For CE systems, only a fraction of the energy is emitted in photons and this assumption of the Eddington approximation is also never satisfied. For systems where the CE phase ends up with the NS merging with its companion's helium core, neutrino emission increases dramatically, allowing the neutron star to incorporate material at the high Bondi rates predicted for these dense conditions. This merger ultimately forms a rapidly accreting black hole that may produce ultralong gamma-ray bursts (Fryer & Woosley 1998).
Although it seems that Bondi-Hoyle accretion assumptions are most applicable to our problem, not all of the energy is converted to neutrinos and photons, some goes into kinetic energy that ejects a fraction of the accreting material (see Section 2.2). It is this ejecta that is the subject of our nucleosynthesis study. We will discuss this ejecta in more detail in Section 2.2, but it is important to understand that the ejecta may also alter the accretion. We will decrease the accretion rate in our Bondi-Hoyle solution to approximate this effect. In this project we use two bounds to match the range of efficiencies in simulations: λ BHL = 1/4, 1/40. The lower value is set to try to capture both asymmetric accretion and ejecta affects that lower the accretion rate.
The accretion during the CE phase depends on the structure of the star which, in turn, depends upon both the stellar mass and evolutionary stage. To estimate the NS accretion rates in the CE phase, we use a coarse grid of stellar models computed using the MESA stellar evolution code (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018, ranging from 8 to 25 solar masses with initial metallicities in the range 10 −4 ≤ Z ≤ 2×10 −2 . The massive star models (M ini ≥ 12 M ) are the ones from Ritter et al. (2018). The 8 M models were computed using the same input physics as in Jones et al. (2013). We refer the reader to those papers for a more thorough description of the stellar evolution calculations.
We study the structure of each of these models as the star evolves, focusing on periods of time when the star is expanding and a CE phase is likely to occur. Figure 1 shows the radius evolution of our stars as a function of time with the points showing the specific times used in our study. As the neutron star spirals into the star, the accretion rate onto it increases. The corresponding accretion rates (assuming λ BHL = 1/4) as a function of radius for the 12, 15 and 20 M stars at different evolutionary times is shown in Figure 2. In the bulk of the envelope, the accretion rate lies between 1 − 10 5 M yr −1 and we will focus on these rates, but if the neutron star spirals into the core, the accretion rate will be higher. Figure 1. Radius evolution of the stellar models used in this study. The x-axis is time remaining until core collapse (or, in the case of some of the 8 M models, envelope ejection). The points indicate at which times stellar structures were used from the models in order to calculate the accretion rates. 10 1 10 3 10 5 10 7 t fin − t

Ejection of Accreted Mass
The accreted material is explosively unstable and early calculations suggested that some of the infalling material would gain enough energy to be ejected (Fryer et al. 1996). Estimates of the convective timescale using the Brunt-Väisäla frequency (see, for example, Cox et al. 1983a;Fryer et al. 1996) suggest that the timescale that the material spends near the proto-neutron star surface is milliseconds in duration. This initial study focused on CE accretion scenarios, but most of the subsequent, more systematic, multidimensional work focused on the higher accretion rates seen in supernova fallback (Fryer et al. 2006;Fryer 2009). Although these studies focused on accretion rates above 10 4 M yr −1 , they showed the same features as the CE models studied in Fryer et al. (1996). The accreted material falls down toward the proto-neutron star surface. A fraction of this material is heated and accelerated to above escape velocities, ejecting it from the system. We designed a set of twelve trajectories based on these simulations, guided by the analytic models developed to understand these simulations. The accreting material accelerates nearly at free-fall until it falls within 10 km of the neutron star surface. The uncertainty in the flow lies in de- Figure 2. Accretion rate for a 1.4 M neutron star as a function of the position of the neutron star within the mass coordinate of its companion for the 12, 15 and 20 M , solar metallicity stars for a range of times. At early times, the envelope is compact and the accretion rate is higher. As the star expands, the accretion rate in the envelope decreases. In the bulk of the envelope, the accretion rate lies between 1 − 10 4 M yr −1 . Here we assume λ BHL = 1/4. The accretion rates may be an order of magnitude lower.
termining how quickly the flow is reversed and material is ejected. We study two extremes: a bounce scenario where the reverse is instantaneous, a convective scenario where the acceleration timescale is on par with the convective turnover timescale. For the latter, convective timescale, we can estimate the acceleration timescale from the Brunt=Väisäla freqency (t acceleration = 1/ω BV ) where where ρ and S are the density and entropy respectively (Cox et al. 1983b). For conditions in supernovae and fallback this timescale is on the range of 2 ms (Fryer & Young 2007). This timescale estimates the growth timescale of Rayleigh-Taylor instabilities, but this provides an approximate timescale for the acceleration timescale to reverse the shock. Fallback simulations suggest that the true answer lies between these two extremes.
We model these two extremes by two parameterized simulations: bounce and convective trajectories. The bounce trajectory assumes that material falls in adiabatically at free-fall until it reaches 20 km (10 km from the surface) and where it is assumed that the material bounce and is ejected at the escape velocity. For the convective trajectory, we assume that the material falls in at free-fall until it reaches a depth of 50 km, where we turn on a force that is strong enough to turn around the trajectory within the roughly 2 ms convective timescale.
The temperature evolution of these two paradigms is shown in Figure 3. The corresponding density evolution is shown in Figure 4. In the bounce trajectory, In our first model, sharp increase of the temperature evolution profile ( Figure 3) mimics a hard stop for the infalling material prior to this expulsion. In our the convective model, the acceleration begins sooner but is more gradual.
Although we expect some gain in entropy during the heating phase, we assume constant entropy evolution for our models. If we increase the entropy, we increase the temperature, and more heavy elements will be produced. Likewise, if we make the evolution even more gradual, the material will remain in a region of high nuclear burning longer, also producing more heavy elements. But our simulation-guided (Fryer et al. 1996(Fryer et al. , 2006Fryer 2009) toy models will provide a gauge of the importance of this ejecta in galactic chemical evolution.
Typically ∼10-25% of the accreted material is ejected along the angular momentum axis (Fryer et al. 2006). It is this material that is the focus of our nucleosynthetic studies. In this paper, we focus on accretion rates between 1 M yr −1 and 10 5 M yr −1 . The rate rarely exceeds our highest accretion rate in common envelop situations and lower rates (occurring in initial phases of the CE) do not produce much nuclear burning. The inner infall, nuclear burning and ejection phases are so rapid (< 1s) that we can assume the accretion rate is constant over any cycle. Table 1 lists the trajectories investigated in this paper, and shown in Figure 3. The trajectories labeled with ".C." are representative of the first model described above, while those labeled with ".D." are representative of the second delayed model.

NUCLEOSYNTHESIS CALCULATIONS FROM NEUTRON STAR ACCRETION
In this section we present the nucleosynthesis yields of the trajectories described in section 2.2 and listed in Table 1. The composition of the accreted material has a scaled solar isotopic distribution (Asplund et al. 2009), with metallicity Z m = 0.02. The nuclear network includes 5234 isotopes and 74313 reactions from the different nuclear physics compilations and rates available (see, e.g., Pignatari et al. 2016). The 3−α and 12 C(α,γ) 16 O by Fynbo et al. (2005) and Kunz et al. (2002) respectively, and the 14 N(p,γ) 15 O by Imbriani et al. (2005). The reaction rate 13 C(α,n) 16 O is taken from Heil et al. (2008), and the 22 Ne(α,n) 25 Mg and 22 Ne(α,γ) 26 Mg rates are from Jaeger et al. (2001) and Angulo et al. (1999). Experimental neutron capture reaction rates are taken when available from the KADoNIS compilation (Dillmann et al. 2006). For neutron capture rates not included in KADo-NIS, we refer to the JINA REACLIB database, V1.1 (Cyburt 2011). The weak rates for light and intermediate-mass species are provided by Oda et al. (1994) or Fuller et al. (1985), by Langanke & Martínez-Pinedo (2000) for mass number 45 < A < 65 and for the weak interaction between protons and neutrons. Finally, for heavy species with A > 65 we use Goriely (1999).
Besides the 14 N(p,γ) 15 O mentioned above, proton capture rates have different sources like the NACRE compilation (Angulo et al. 1999) and Iliadis et al. (2001) for isotopes in the mass region between 20 Ne and 40 Ca. Proton Trajectory Name Peak T 9 C Trajectory Name Peak T 9 D Accretion Rate (M yr −1 )  Table 1. List of trajectories investigated in this paper, along with peak temperatures in those trajectories. The first column includes those trajectories in Figure 3 which have undergone a sudden change in direction during the accretion, while those with a more gradual turnaround are listed in the third column. The accretion rates range from 1 to 10 5 M yr −1 as detailed in section 2. Temperatures have been clipped in these models at 10 GK, as reaction rates tables above this threshold are not available.
captures on isotopes heavier than 40 Ca are given by the JINA REACLIB database. The 13 N(p,γ) 14 O rate is taken from Caughlan & Fowler (1988), and proton captures on 27 Al from Champagne & Wiescher (1992). For our simulations we have used material accreted with solar composition. This approximation does not affect conclusions in this paper. Indeed, at the lower accretion rates seed nuclei might have large effects on the distribution of products, but their impact on the final yields are marginal. On the other hand, at the higher accretion rates investigated in this paper the material accreted enters in nuclear statistical equilibrium at extreme conditions (Table 1), so that the final yields are largely insensitive to the initial isotopic distribution (and to the electron fraction y e ). As we will see in the following sections, nucleosynthesis in these trajectories is dominating the total integrated ejecta for GCE.

Yields at Differing Accretion Rates
Production of isotopes at various accretion rates onto the surface of the compact object were investigated. Accretion rates varied from 1 M yr −1 , where the neutron star first comes into contact with the companion, to the maximum accretion rate of 10 5 M yr −1 , corresponding to the final stages of accretion in the system when the neutron star reaches the helium core.
Tables 2 and 3 summarize the results of this section and give a list of the isotopes with the highest overproduction, the increase in the abundance of an isotope relative to its initial abundance in the accreted material, at the accretion rates 10 4 M yr −1 and 10 5 M yr −1 for the trajectories investigated. Trajectories at lower accretion rates do not show the same extent of overproduction, and therefore contribute minimally to our GCE results. Table 4 shows the top 5 isotopes overproduced in each of the trajectories investigated. Those isotopes most produced in these trajectories are also given in the relevant sections. In Tables 6 and 7 the isotopic yields of stable isotopes (with contributions from decayed unstable isotopes) are provided for the trajectories listed in Table 1. Complete radiogenic contribution is considered for these abundances. Undecayed abundances are also provided in separate tables available on-line at CANFAR. The corresponding overproduction factors are given in Table 8 and  table 9.
Complete isotopic production factors and elemental abundances are shown in Figures 5 and 6, respectively, where results from the two sets of trajectories at accretion rates can be compared (see discussion here below). Figure 7 shows  Table 3. Table showing those isotopes overproduced in the given ranges for each of the accretion rates shown in the mod.D trajectories the integrated fluxes for each of the Mod.C accretion rates. Figures 8 and 9 show the final isotopic production factors, zoomed in the mass region 40 A 100, and the electron fraction Ye obtained. Our results show that the largest contributions to GCE must come from those trajectories at higher accretion rates. In these conditions, hydrogen is fully burned allowing high abundance overproductions at the iron group and beyond. Lower accretion rates contribute marginally to the enrichment of ejected material, both because the material is incompletely burnt and the amount of material ejected is much smaller. Proton rich material is highly overproduced in trajectory mod.C.ar1d4, and this is the source of most of the proton rich material ejected from the system. This trajectory reaches high enough temperatures that the accreted fuel burns efficiently, with a low enough peak density such that electron captures do not become dominant causing a shift in the Y e and peak production to more neutron rich isotopes.  Table 4. 5 most overproduced isotopes for each of the accretion rates investigated.
3.1.1 Accretion rate 1 M yr −1 : mod.C.ar1d0 and mod.D.ar1d0 The dominant reactions are proton captures from the H rich accreted material, followed by β + decays bringing the material back towards the valley of stability in the mod.C.ar1d0 case. Some (p,α) reactions are evident in this region, hindering the flow of material to heavier masses. Significant increases in the abundance of 7 Li can be seen for both mod.C.ar1d0 and mod.D.ar1d0 (top left panel of Figure 5) with both trajectories showing overproduction factors of ∼10 6 . This occurs due to efficient 3 He (α, γ) 7 Be, followed by β + decays. 15 N and 18 O are also increased significantly, by factors of about 150 and 250 for mod.C.ar1d0, and 280 and 510 for the delayed trajectory. Enhancements of light and intermediate-mass elements are comparable for the two trajectories, until the mass region A ∼ 50. Beyond this mass, the delayed trajectory does not show efficient production for the heavier elements, while in mod.C.ar1d0 final abundances are enhanced up to A ∼ 80. This is due to the higher peak temperature reached compared to the mod.D.ar1d0 model (1.303 GK compared to 0.728 GK), causing a build up of heavier nuclei mostly via proton capture reactions.
Panel a of Figure 6 shows that neither of the trajectories at this accretion rate allows any burning of material above Z=40. The conditions at this accretion rate are not extreme enough to allow proton captures on these heavier nuclei, or to cause photo-disintegration on the heavier elements present in the accreted material.
74 Se, 78 Kr and 84 Sr (the lightest of the p-nuclei) all have overproduction factors greater than 10 2 for mod.C.ar1d0 at this accretion rate, where as the delayed trajectory shows no production of these isotopes. Inspection of flux charts for mod.C.ar1d0 shows a proton capture path with β + decays in the proton-rich side of the valley of stability, like for the rp-process (Schatz et al. 2001). This becomes most evident above the N=20 neutron magic number (see panel (a) of Figure 7). As can be seen from the first columns of Tables 5 and 9, the H fuel accreted from the companion star remains largely unburnt under these conditions -only 0.2% of the initial abundance of H is burnt for both of these trajectories.
3.1.2 Accretion rate 10 M yr −1 : mod.C.ar1d1 and mod.D.ar1d1 Nucleosynthesis in mod.C.ar1d1 results in overproduction of a large number of isotopes up to A∼157, due to the higher peak temperatures experienced in this trajectory as compared with the mod.C.ar1d0 case, allowing for more efficient proton captures in the accreted material, whilst remaining below the threshold for activation of photodisintegration of heavier material seen in higher accretion rate models. The highest mass isotope to be overproduced under these conditions is 180 Ta (with a production factor of 5.4×10 2 ), and the isotopes with the greatest overproduction are 138 La and 31 P (3.3×10 3 and 1.0×10 3 times initial abundance respectively). 14 other isotopes -33 S, Ce all have production factors of between 10 2 and 10 3 . In mod.D.ar1d1 abundances are greatest at lower masses. 7 Li, 21 Ne, 23 Na, 33 S, 42 Ca and 84 Sr have overproduction factors of between 10 2 and 10 3 . The much reduced range of highly overproduced isotopes is again due to the lower temperatures experienced by the delayed trajectory on infall of the material. The peak temperature of mod.D.ar1d1 is similar to the peak temperature of mod.C.ar1d0 trajectory discussed before (see Table 1). Differences in the two abundance distributions are due to the temperature histories of the two trajectories, with the delayed trajectory exposed for longer to more extreme conditions. Production factors and abundances are shown in the top right panels of Figures 5 and 6, respectively. The mass fraction of the majority of stable isotopes above Z=40 remains unchanged, as was the case for the previous section. Peak temperatures from our simulations for these trajectories do not increase beyond ≈2 GK ( Figure 3) and burning occurs on a timescale of order milliseconds. Neither α nor proton captures have high enough probability at these temperatures to trigger complete burning of the accreted material.
Flux nucleosynthesis plots for this accretion rate are similar to those for the 1 M yr −1 case, and are shown in panel (b) of Figure 7.
3.1.3 Accretion rate 10 2 M yr −1 : mod.C.ar1d2 and mod.D.ar1d2 Abundance production for mod.C.ar1d2 is clustered at mass between A ∼ 25 and 100. 51 V and 61,62 Ni are overproduced by factors larger than 10 3 . Eight other isotopes have production factors greater than 10 2 : 42 Ca, 48 Ti, 52 Cr, 58,60 Ni, 64 Zn and 96,98 Ru. The flux nucleosynthesis plot (panel (c) of Figure 7) for this trajectory shows a relevant difference compared to lower accretion rates. Burning proceeds further from the valley of stability, due to the more efficient proton captures at high temperatures. The α-captures on intermediate-mass isotopes are activated. Heavier isotopes are destroyed by (γ,n) photodisintegration reactions, and a mild flux of (n,γ) reactions is also activated. As for lower accretion rates, a large proportion of the infalling material is not burnt (Table 6) with only ≈0.3% of the hydrogen fuel being consumed. The peak temperatures in this trajectory is 3.981 GK.
The delayed trajectory mod.D.ar1d2 shows an abundance distribution extremely different from mod.C.ar1d2.   Table 1 are shown. Abundances for delayed trajectories ("mod.D.") are shown in red, and in black for "mod.C." trajectories. Relative accretion rates are indicated in each panel.
Production extends to higher mass isotopes due to the longer time that the material spends at high temperature, and without yet relevant activation of photo-disintegration reactions. There are 22 isotopes enhanced by a factor between 10 2 and 10 3 , and among them 8 have masses above A = 150. 38 Ar, 41 K, 43,44 Ca, 45 Sc, 47,48 Ti and 136,138 Ce show overabundances between 10 3 and 10 4 , and 42 Ca greater than 10 4 . Temperatures are not high enough to initiate the photodisintegration reactions observed in mod.C.ar1d2.
It can be seen from the middle left panel of Figure 6 that there is a significant change in the mass fraction of elements with 50 < Z < 65, although the very heaviest elements (with Z > 65) which are investigated in this work remain largely unchanged in their abundances -The distributions of both of the C and D model are identical in the upper panels of Figure 6, implying that the distribution cannot be changed from the initial abundances. This is verified by comparison with Figure 5, where the production factors for all isotopes above A = 150 and most above A=100 are negligible. Some elements with 50 > Z > 65 have abundances increased by up to 3 orders of magnitude. However, because of their small intrinsic abundances and the small contribution to the final integrated ejecta, this effect is not important to GCE models.
3.1.4 Accretion rate 10 3 M yr −1 : mod.C.ar1d3 and mod.D.ar1d3 Material at this accretion rate reaches a peak temperature of 6.952 GK in the mod.C.ar1d3 trajectory and 4.090 GK in the mod.D.ar1d3 case, although for both trajectories the accreted material remains largely unburnt. The trajectory mod.C.ar1d3 reaches nuclear statistical equilibrium, and the destruction of heavier elements by photodisintegration as already seen for mod.C.ar1d2 now occurs in both trajectories. Isotopes are produced most efficiently in the mass region  Figure 6. Element abundance distributions in mass fraction with respect to atomic number are shown for the same cases in figure 5. Contribution from radioactive decay is considered. Abundances for delayed trajectories ("mod.D.") are shown in red, and in black for "mod.C." trajectories. Most of the element abundances in the first two panels remains unchanged from their initial distributions, while significant changes can be observed at higher accretion rates.

50
A 100, including species in the iron-group region and up to Ru. For both trajectories, a significant increase of greater than 10 3 in the abundances of 58,61,62 Ni is seen. In mod.C.ar1d3 overabundances greater than 10 2 are obtained for 5 other isotopes: 48 Ti, 42 Ca, 60 Ni, 64 Zn and 51 Va, and 9 in the delayed trajectory mod.D.ar1d3: 45 Sc, 42,44 Ca, 48,49 Ti, 52 Cr, 60 Ni, 80 Kr and 96 Ru. Elemental distributions of the two trajectories for this accretion rate are similar (middle left panel of Figure 6) up to Z = 30, however there are significant differences in the abundances of Mo, Nb and Ru, none of which are overproduced in mod.C.ar1d3, and for some isotopes of Kr and Sr where a difference in overproduction factors of over an order of magnitude can be observed.
3.1.5 Accretion rate 10 4 M yr −1 : mod.C.ar1d4 and mod.D.ar1d4 Among the trajectories described in this section, mod.C.ar1d4 is the first case to undergo complete burning of hydrogen, with this trajectory entering nuclear statistical equilibrium (NSE). Significant overproduction for a large number of isotopes is observed: 74 Se, 78 Kr, 84 Sr and 90,91 Zr have final abundances over 6 orders of magnitude higher than initial, 92 Mo is overproduced by a factor of 2.6×10 7 . Excepting the Zr isotopes, all of these nuclei are classically defined as products of the p-process (Arnould & Goriely 2003;Rauscher et al. 2013;Pignatari et al. 2016). While these are proton rich isotopes, for mod.C.ar1d4 we obtain a final neutron rich distribution (Ye = 0.46, table 6). This is mostly due to the neutron-rich abundance signature in the Ni region, with 58% of mass fraction as 62 Ni and only a minor contribution to 58 Ni. The zoomed isotopic distribution is shown in Figure 8, bottom left panel. The abundance pattern of heavy isotopes is similar to neutrino-driven winds ejecta with proton-rich composition (e.g., Fröhlich et al. 2006;Roberts et al. 2010;Arcones & Montes 2011). Partial burning of the hydrogen is observed for the delayed trajectory mod.D.ar1d4. Production for heavy isotopes is marginal compared to mod.C.ar1d4 (Figure 5), while intermediate-mass elements are made more efficiently (Figure 6). Isotopes 58,61 Ni both have overproduction factors greater than 10 4 . As also shown in the bottom left panel of Figure 9, p-process isotopes are not made in mod.D.ar1d4, with no significant production above A ∼ 80. As we have seen for trajectories in the previous section, material above this mass has been destroyed by photodisintegration. Panel (e) of Figure 7 shows the flux plot for mod.C.ar1d4. The material is in NSE, with the large amount of hydrogen allowing proton captures to occur and extending the abundance distribution to the proton rich side of the line of stability. The lower densities reached in this model as compared with the mod.C.ar1d5 trajectory ensure that electron capture reactions are not favorable enough to reduce the Y e .
The trajectory is in NSE with large fluxes in the iron group region, lower average and lower peak densities in this model lead to a more proton rich distribution of products than the mod.C.ar1d5 case (see panel f of Figure 7 for comparison). The flux plot is dominated by the reactions in NSE equilibrium around the iron group region during the NSE phase. In these simulations, the temperature and density freezout is extremely fast and do not allow to significantly modify the integrated fluxes after the trajectories exit from NSE.
3.1.6 Accretion rate 10 5 M yr −1 : mod.C.ar1d5 and mod.D.ar1d5 More neutron rich material is produced under these conditions with isotopic distributions skewed towards more neutron rich isotopes. This is due to the higher peak and average density conditions experienced in this model as compared with the mod.C.ar1d4 trajectory. The most overproduced isotope in trajectory mod.C.ar1d5 is 88 Sr, with an overproduction factor of 9.8×10 5 along with 87 Rb and 62,64 Ni, all with overproduction factors of greater than 10 5 . A similar distribution is observed in the mod.D.ar1d5 model, however 64 Ni is the most overproduced isotope in this case, with 87 Rb, 86 Kr, 88 Sr and 54 Cr being the other isotopes with a production factor greater than 10 5 . A large number of isotopes in both the mod.C.ar1d5 and mod.D.ar1d5 models are overproduced by factors of ∼ 10 4 -9 for the mod.C.ar1d5 case and 12 for the mod.D.ar1d5, all of which are clustered around the iron group region Temperatures in mod.C.ar1d5 have been clipped at 10 GK as reaction rate tables are not available beyond this value. The overall distribution and range of production is sim-ilar between mod.C.ar1d5 and mod.D.ar1d5, as can be seen in the bottom right panel of Figure 6 and Figures 8 and 9. As discussed by Fryer et al. (2006), fall-back trajectories can produce r-process abundances with mild neutron-rich conditions. In the scenario discussed here, within a realistic range of accretion rates, mod.C.ar1d5 and mod.D.ar1d5 both show neutron-rich Y e (0.447 and 0.443, respectively) and an abundance signature similar to the weak r-process (e.g., Seeger et al. 1965;Kratz et al. 1993;Arcones & Montes 2011;Wanajo 2013). Collapsars, as investigated by Siegel et al. (2018), with accretion rates from 0.3 to 30 times those investigated in these ar1d5 trajcectories have also been shown to have an increase in neutron rich material.
Panel (f) of Figure 7 shows the integrated flux plot for the mod.C.ar1d5 trajectory. It can be seen that the reaction pathways for this trajectory are through more neutron rich isotopes than in the mod.C.ar1d4 case. This is due to the higher density over the period of infall and ejection, making electron capture decays more efficient. In contrast to the mod.C.ar1d4 case (panel (e) of Figure 7), the material is more neutron rich, due to the higher peak density and increased time spent at higher densities in this trajectory, favoring electron capture reactions in the accreted material.

IMPLICATIONS FOR GALACTIC CHEMICAL EVOLUTION
In the previous section we have explored the large variety of nucleosynthesis that can be found in material infalling close to the neutron star at different accretion rates, and ejected back in the stellar host, and into the interstellar medium following evolution of the companion star. Both light and heavy elements can be made, showing proton-rich or neutron-rich isotopic patterns by slightly changing the trajectory conditions within a realistic parameter space. Interestingly, all of the different nucleosynthesis patterns shown could be produced in the same merging event, during the evolution of the merging system. The toy model presented in section 2 made to access the possible nucleosynthesis found in these systems is meant to be a first step to explore the production of elements in these systems, and motivate fully resolved hydrodynamics simulations. Based on the calculations present here, we can now also verify if these systems are relevant also for galactic chemical evolution (GCE).

Population synthesis study of CE phases
Before studying the role of these stellar objects in a galactic chemical evolution context, we need to determine the yields mass ejection in CE, and couple the accretion rates from Section 2 with the properties of the CE systems. To study the population of CE phases, we employ the StarTrack population synthesis code (Belczynski et al. 2002(Belczynski et al. , 2008 to generate a population of binary compact objects . The code is based on revised formulas from Hurley et al. (2000); updated with new wind mass loss prescriptions, calibrated tidal interactions, physical estimation of donor's binding energy in CE calculations and convection driven, neutrino enhanced supernova engines. A full description of these updates is given in Dominik et al. (2012b). The two most recent updates take into account measurements of initial parameter distributions  for massive O stars (Sana et al. 2012b) as well as a correction of a technical bug that has limited the formation of BH-BH binaries for high metallicity (e.g., Z = 0.02).
With this code, we modeled the binary interactions of 500,000 stars for 3 different metallicities (Table 5). Only a small fraction of the systems actually go through a CE phase. For each of these systems, we calculate the progenitor star mass, the radius at the onset of the CE phase, and the final separation. Using the progenitor mass and radius at the onset of the CE phase, we can determine which of our MESA progenitors to use and the time in its evolution. There is a slight inconsistency in our approach, because the formulae for the stellar radii versus time in the population synthesis models are not identical to our MESA models, but this discrepancy allows us to make a first approximation of the accretion rates. Once we determine the mass and time of the CE interaction, we can use the final separation to determine how deep the neutron star inspirals into its companion.
With our stellar structures and the separation evolution, we can follow the full range of accretion rates in the CE phase. During the CE phase, the NS is continuously accreting and ejecting mass, evolving through a range of accretion rates as the NS spirals down toward the stellar core. The accretion and ejection is rapid compared to the CE phase, and we can approximate this evolution as a series of phases with different accretion rates. For each binary system in our population synthesis calculation, we know the mass of the stellar component, the evolutionary timescale of the star at the initiation of CE phase (from the initial binary separation) and the range of accretion rates (based on the final binary separation and our stellar models). We then assume that 25% of the accreted mass is reheated and ejected based on the results of accretion simulations (Fryer et al. 2006) to get the rate of mass ejection.
To get the total amount of mass ejected, we integrate the mass ejection rate over the inward spiral of the NS. Based on CE simulations (Ricker & Taam 2008, 2012Passy et al. 2012;Ivanova et al. 2013), we assume that the duration of a typical CE phase persists for roughly 3 times the orbital period (P onset ) at the onset of the CE phase. In our calculations, we will assume that the time spent at each radius (t r )  is 3 times this orbital period (P r ) at radius r, for each bin i, normalized by the orbital period of each radius: In this paper, we consider 3 options for the mass accretion and ejection: λ BHL = 1/4, 25% mass ejecta, λ BHL = 1/40, 25% mass ejecta, and λ BHL = 1/40, 10% mass ejecta. Recall that λ BHL is a non-dimensional parameter relating to the efficiency of accretion in our system (see section 2.1). Figure 11 shows the distribution of accretion rates for 10 of the close binaries in our population synthesis calculation with our model using λ BHL = 1/40, 25% mass ejecta.
Using our full population of binary systems, we can estimate the average mass ejected per binary system as a function of metallicity and our value for λ BHL (Figure 12). The fraction of merging systems in our models are 34%, 61% and 75% for 1, 0.1 and 0.01 solar metalicities repectively. On average, a binary system ejects a few hundredths of a solar mass of material (in practice, some systems can eject nearly a solar mass of material while many systems eject very little mass).
Our MESA stellar models consist of a course grid in mass and time, which presents challenges in the outer layers of the companion star in matching with our population synthesis models, extracted from Hurley et al. (2000). This is not a concern however at higher accretion rates near the core of the companion, where the majority of burning in these models is observed, so long as the core masses of our stars are similar. As we focus on higher accretion rates for our yields, inconsistencies introduced by this are minimal in our models.

Galactic Chemical Evolution
To estimate the contribution of CE mass ejection events in a galactic chemical evolution context, we use the OMEGA code described in Côté et al. (2017). This is a classical onezone open-box model (e.g., Tinsley 1980) that is part of the open-source NuPyCEE package 1 . We adopt the same default Milky Way setup as in Côté et al. (2018). Our input parameters are tuned to reproduce the star formation rate, the gas fraction, the gas inflow rate, and the Type Ia and core-collapse supernova rates currently observed in the Milky Way (see Table 1 in Kubryk et al. 2015), within the observational errors. Our model is also calibrated to reach solar metallicity (Z m = 0.014, Asplund et al. 2009) when the Figure 10. Initial versus final separation for all the binary systems undergoing a CE phase from our population synthesis simulations, given at 3 metallicities. Figure 11. Mass ejected versus accretion rate for 10 sample binary systems in our population synthesis calculations. As the neutron star spirals into its companion, the accretion rate increases. Throughout this inspiral phase, mass is ejected and the conditions of this ejecta evolve as the neutron star spirals into deeper and deeper stellar layers. The peak accretion rate depends upon the structure of the stellar core and the depth of the inspiral. Figure 12. Average mass ejected per binary system (based on the population synthesis models) as a function of metallicity (solar, 1/10th solar and 1/100th solar) for two extremes in our assumptions for accretion rate and mass ejected ("high" denotes λ BHL = 1/4, 25% mass ejecta, low denotes λ BHL = 1/40, 10% mass ejecta). Some binaries will not eject much mass at all, others will eject over a solar mass.
Sun forms, which is 4.6 Gyr before the end of the simulation (Connelly et al. 2017). The evolution of metallicity in our model is generated using NuGrid yields (Ritter et al. 2017) for massive and low-mass stars, and the yields of Iwamoto et al. (1999) for Type Ia supernovae.
To include the contribution of CE events, we use the NuPyCEE delayed-extra source implementation 2 which allows to include additional enrichment sources based on input metallicity-dependent delay-time distribution (DTD) functions and yields. Because the high NS accretion rate CE events studied here occur in systems involving two massive stars, the mass ejection rate is assumed to follow the lifetime of massive stars. In practical terms, for each stellar population formed in our model, all CE events occur between 5 and 40 Myr following the formation of the progenitor stars. This first-order implementation will be improved in follow-up studies. For the chemical composition of CE events ejecta, we convolve our nucleosynthesis calculations (see Section 3.1) with the metallicity-dependent mass ejection rates inferred from our population synthesis analysis (Section 4.1). An example of the resulting yields is shown in Figure 13. The complete set of yields used in our chemical evolution calculations is available on-line at CANFAR. Figure 14 shows the contribution of CE events to the solar isotopic composition predicted by our models, using our three mass accretion and ejection options described in Section 4.1: λ BHL = 1/4, 25% mass ejecta, λ BHL = 1/40, 25% mass ejecta, and λ BHL = 1/40, 10% mass ejecta. We also ] mod.C mod.D Figure 13. Yields ejected by CE events on average per massive star binary (MSB) at Z = 0.02, representing the convolution of our population synthesis models (Section 4.1) with the yields calculated for different accretion rates (see Section 3.1), using the trajectories "mod.C" (black line) and "mod.D" (red line). This is for λ BHL = 1/4 and 25% mass ejecta (see Section 4.1). The complete set of yields used in our galactic chemical evolution calculations for different metallicities and mass accretion and ejection options is available on-line at CANFAR.
varied the binary fraction of massive stars between 25 % and 100 %. Overall, the different mass accretion and ejection options cause more variations in our predictions than the binary fraction. Those figures show that in some cases, CE events could significantly contribute to the chemical evolution of some iron-peak and first-peak neutron-capture isotopes in the Galaxy. In fact, some isotopes are overproduced in our models relative to the solar composition. Our result here show only the contributions from common envelope events to the solar composition, as direct comparison to isotopic abundances of other contributing events depend sensitively on the models chosen. As mentioned in the previous sections, the nucleosynthesis of CE events is sensitive to the physical conditions (i.e., temperature and density). For example, when using the yields from the trajectories "mod.C.", the p-isotopes 74 Se, 78 Kr, 84 Sr, and 92 Mo, as well as 90,91 Zr, are always overestimated by more than an order of magnitude. However, when using the yields from the delayed trajectories "mod.D.", none of these isotopes are significantly produced. Instead, the isotopes contributing the most to the solar compositions are rather in that case concentrated on the neutron-rich side, such as 64 Ni, 70 Zn, and 86 Kr.
We note that the nucleosynthesis has been calculated assuming an initial metallicity Z = 0.014, and that we used those yields for all metallicities. The purpose of our chemical evolution calculations is to verify whether or not CE events could be important for galactic chemical evolution. Our results should be seen as a first-order approximation and a motivation for future work.

FINAL DISCUSSION AND CONCLUSIONS
In this work we have explored the nucleosynthesis produced by neutron stars accreting in binary common envelopes. A realistic range of accretion rate conditions is explored (between 1 M yr −1 and 10 5 M yr −1 ), for two sets of trajec-tories. In the first set, material is assumed to be near freefall, and suddenly ejected with acceleration comparable (and in opposite direction) to the gravitational acceleration. In the second set, material is ejected with more gradual acceleration, resulting in lower temperature and density peaks compared to the first set. A large variety of nucleosynthesis patterns were obtained. Heavier elements are produced with increasing accretion rate, due to the higher temperature and density peaks. In particular, weak interactions become extremely important in defining the final composition at accretion rates of 10 4 -10 5 M yr −1 , and leading to a proton-rich or neutron-rich nucleosynthesis pattern for heavy elements between Fe and Ru.
We test for the first time the impact of CE events in a galactic chemical evolution context. We find that accreting neutron stars could contribute in a non-negligible way to the solar composition for some isotopes. In particular, using yields from the first set of trajectories, we overproduce many p-isotopes such as 74 Se, 78 Kr, 84 Sr, and 92 Mo, among others. Using the second set of trajectories where there is a gradual change in trajectory, these events do not contribute to the solar abundances of these isotopes. The yields from these events is therefore highly dependent on the specific physical conditions experienced in the CE event, and the conditions in which the nucleosynthesis takes place is therefore crucial for galactic chemical evolution.
To summarize, we have shown that neutron stars accreting in binary common envelopes are potentially an important (unaccounted for) nucleosynthesis site for the chemical evolution of the Galaxy. Due to the simple approximations made in this first study for the nucleosynthesis trajectories, present sets of yields are still highly uncertain. But these results are a first important step, demanding more detailed simulations in the future.  Figure 14. Contribution of CE events on the solar isotopic (mass number) composition of the Sun, using the galactic chemical evolution framework described in Section 4.2 and the yields from the trajectories "mod.C" (top) and "mod.D" (bottom). Different colors represent our different mass accretion and ejection options (see Section 4.1), where λ BHL defines the magnitude of the accretion rates (see Equation 3) and the percentage represent the fraction of accreted mass that is reheated and ejected. For each color, the band represents the range of solutions assuming different binary fractions from 25 % to 100 % for massive stars. The dotted horizontal line shows a contribution of 100 %. Anything above this line implies an overestimation relative to the solar composition. Isotopes of interest which are overproduced are labeled.  Physical Review Letters, 96, 142502 Fryer C. L., 2009, ApJ, 699, 409 Fryer C. L., Woosley S. E., 1998, ApJ, 502, L9 Fryer C. L., Young P. A., 2007, ApJ, 659, 1438Fryer C. L., Benz W., Herant M., 1996, ApJ, 460, 801 Fryer C. L., Colgate S. A., Pinto P. A., 1999a, ApJ, 511, 885 Fryer C. L., Woosley S. E., Hartmann D. H., 1999b Table 7. The same as Table 6, but for trajectories mod.D.ar1d0, mod.D.ar1d1, mod.D.ar1d2, mod.D.ar1d3 mod.D.ar1d4 and mod.D.ar1d5 respectively.