Simulations of stripped core-collapse supernovae in close binaries

We perform smoothed-particle hydrodynamical simulations of the explosion of a helium star in a close binary system, and study the effects of the explosion on the companion star as well as the effect of the presence of the companion on the supernova remnant. By simulating the mechanism of the supernova from just after core bounce until the remnant shell passes the stellar companion, we are able to separate the various effects leading to the final system parameters. In the final system, we measure the mass stripping and ablation from, and the velocity kick imparted to, the companion star, as well as the structure of the supernova shell. The presence of the companion star produces a conical cavity in the expanding supernova remnant, and loss of material from the companion causes the supernova remnant to be more metal-rich on one side and more hydrogen-rich (from the companion material) around the cavity. Following the removal of mass from the companion, we study its subsequent evolution and compare it with a single star not subjected to a supernova impact.


Introduction
There is substantial evidence that most massive stars evolve in binary systems (Duquennoy and Mayor, 1991;Rastegaev, 2010;Sana et al., 2012). Therefore, the presence of companion star is an important consideration in the theory and observation of supernovae (SNe) and supernova remnants (SNRs). In particular, while Type Ia (white-dwarf; WD) supernovae may have a companion which has deposited sufficient mass onto the WD to trigger a 'single-degenerate' explosion, many Type Ib/c (stripped core-collapse) supernovae may have close companions that have been at least partly responsible for the loss of mass from the progenitor (Bersten et al., 2014;Eldridge et al., 2015;Fremling et al., 2014;Kim et al., 2015;Kuncarayakti et al., 2015).
Observational searches for supernova companions have typically focused on Type Ia explosions. Possible companions have been a subject of scrutiny in order to determine the frequency of the two main suspected (single-degenerate or double-degenerate) explosion channels (Maoz et al., 2014). Hydrogen enrichment from a companion has been searched for in Type Ia SNRs, but so far there has been no evidence of hydrogen lines (Leonard, 2007;Lundqvist et al., 2015; * Correspondence: rimoldi@strw.leidenuniv.nl Leiden Observatory, Leiden University, Niels Bohrweg 2, 2333 CA Leiden, Netherlands Full list of author information is available at the end of the article Mattila et al., 2005). As noted in García-Senz et al. (2012), detection of H α lines may be difficult due to confusion with Fe and Co lines due to the mostly slow (< 10 3 km s −1 ) hydrogen mixing with iron-peak elements.
The presence of a supernova companion is difficult to directly detect if they are low-mass stars at very large distances, and so far definitive evidence of close companions to any SN progenitor, let alone those of Type Ib/c, has been lacking. Tycho G is probably the best example of a directly imaged, suspected companion, associated with the galactic Type Ia SN, Tycho (SN 1572;Ruiz-Lapuente et al., 2004), though recent observations put its status as a supernova companion into dispute (Kerzendorf et al., 2013;Xue and Schaefer, 2015). On the other hand, some direct searches for single-degenerate companions have ruled out giant/subgiant (evolved) stars (SN 2011fe andSNR 1006;González Hernández et al., 2012;Li et al., 2011) and even main-sequence companions (SNR 0509-67.5; Schaefer and Pagnotta, 2012).
The presence of a companion due to increased emission, and therefore modification of the standard light curve, from the ejecta interacting with the companion has also been ruled out in observations of Type Ia supernovae (Olling et al., 2015). However, a recently observed supernova (iPTF14atg; Cao et al., 2015) does show evidence of interaction with a companion arXiv:1510.02483v1 [astro-ph.SR] 8 Oct 2015 through the detection of an ultraviolet burst in the first several days. Though much of the focus of previous work has been on Type Ia explosions, the phenomena of companion interactions with single-degenerate Type Ia ejecta has parallels with core-collapse supernovae in binaries, and therefore this scenario still provides a useful context. Tauris and Takens (1998, hereafter, TT98) analytically investigated the consequences of a SN in a close, circularised binary, with the goal of finding the runaway velocities of the components of a binary disrupted by a Type Ib/c SN. This work was calibrated by estimations from early simulations of the effect of a SN shell impact on a star (Fryxell and Arnett, 1981) in order to determine the amount of mass lost and the change in velocity of the companion. Motivated by this problem, we perform simulations of SNe in binary systems with properties comparable to those used in TT98.
Simulations of supernovae have been performed at many scales, investigating regions of hundreds of kilometres around the nascent neutron star in the early supernova evolution (for a review of recent progress, see Janka, 2012) to effects of the impact on a companion, or the influence of the companion on the overall structure of the ejecta. The impact of Type Ia ejecta on companions has, in particular, been well investigated (Liu et al., 2012;Marietta et al., 2000;Pakmor et al., 2008;Pan et al., 2012). Hirai et al. (2014) investigated the fraction of mass stripped from a giant companion star due to a effect of a core-collapse (Type II) SN using a two-dimensional grid-based Eulerian code. Recently, Liu et al. (2015) also presented results on the consequences of a Type Ib/c supernova interacting with a binary companion using smoothed-particle hydrodynamics (SPH). These studies, however, have often focused on the companion star without following the explosion self-consistently from the moment of the supernova. As a consequence, the ejecta shell is initialised artificially, via analytic prescriptions, near the surface of the companion, without considering its earlier evolution. In addition, the response of the binary companion, and subsequent supernova remnant evolution, is analysed in these cases from a static configuration rather than placing the binary in an orbit.
For close binary orbits it is typically assumed that the binary has circularised by this point in its evolution, so that the eccentricity of the orbit can be set to zero (TT98). We follow the same assumption in this work. Moreover, despite these close separations and, therefore, short orbital periods, in theoretical work the binary period is taken to be much shorter than the timescale over which the ejected shell impacts the companion. This can be made more explicit (as in, for example, Colgate, 1970) by noting that the ejecta velocity must be larger than the escape velocity of the primary star, where M 1 and R 1 are the mass and radius of the primary. Since the distance a of the companion from the primary is larger than R 1 , and since the orbital velocity at that distance is then it must be that v ej > v orb . In practice, the typical ejecta velocities (10 3 ∼ 10 4 km s −1 ) are much larger than the orbital velocities (∼ 10 2 km s −1 ), hence the latter is typically ignored in analytic velocity calculations. However, matter in the ejecta in fact have a radially dependent velocity (approximately, in the homologous regime, v ej (R, t) ∝ R/t). Therefore, although the density also varies through the shell, during the latetime interactions of some lower-density, slower (and presumably high-metallicity) ejecta with the companion, we may no longer be justified in ignoring the orbital velocity. An additional important factor in the dynamics of supernovae in binaries is a possible kick imparted to the newly formed neutron star. This is likely due to a 'gravitational tugboat' effect from asymmetry in the ejecta surrounding the neutron star after the core bounce, and perhaps also high magnetic fields and the asymmetric emission of neutrinos from the protoneutron star (Kusenko and Segrè, 1996;Maruyama et al., 2011;Scheck et al., 2004Scheck et al., , 2006Wongwathanarat et al., 2013). For ultra-stripped supernova progenitors, which have very small ejecta masses, the shock expands very rapidly and the tugboat effect on the neutron star has been shown to be minimal (Suwa et al., 2015). For the range of hydrodynamic simulations in this paper we do not apply any additional kick to the neutron star.
To study this problem hydrodynamically, we simulate a supernova in an orbiting binary self-consistently from just after the moment of core bounce in the supernova. To this end, we first generate stellar structure models of the binary components using a onedimensional stellar evolution code, where we strip an evolved massive progenitor of the majority of its envelope. We then convert these stellar structures to threedimensional stars in an SPH code, and run simulations from the moment of the SN. We vary the mass of the primary star as well as the orbital separation independently. In particular, we are interested in investigating the dependence of the companion's removed mass and impact velocity on the initial orbital separation of the binary. We describe our numerical method in more detail in the following section.

Method
Throughout this work we use the Astrophysical Multipurpose Software Environment [1] (AMUSE; Pelupessy et al., 2013;Portegies Zwart et al., 2009 to perform our simulations. We first outline the technique used to generate the stellar models in our binary systems (Section 2.1). We then describe the set up of the initial hydrodynamical models from the stellar structure (Section 2.2). Finally, we describe the simulation of the SN in the binary, with some discussion on the initial convergence tests that were performed (Section 2.3).

Stellar models
In order to generate an SPH realisation of the binary, we require a stellar evolution code that can return the internal structure of the star. Two evolution codes in AMUSE fit this criterion: MESA (Paxton et al., 2011) and EVtwin (Eggleton, 2006(Eggleton, , 1971. We chose MESA to evolve the models to their final stellar structure, motivated by difficulties in previous work in using EVtwin to obtain solutions past the carbon flash in more massive stars (de Vries et al., 2014).
Due to interactions with the binary companion (and potentially also through stellar winds), much of the mass of the primary star is lost over its lifetime, resulting the helium star progenitor of a Type Ib/c supernova (for some observationally-motivated examples, see Kim et al., 2015). To obtain an estimated lifetime of the progenitor, we first use SSE, which is a fast predictor of stellar properties based on parametrised stellar evolution tracks (Hurley et al., 2000). With the intent of generating 3 M and 4 M Helium-star progenitors, we begin with a 12.9 M and 16.0 M zero-age models with metallicity Z = 0.02, [2] which are predicted by SSE to end with the required helium core masses. We do not model or speculate on the specific mechanisms of the mass loss from the primary star, but instead apply a constant mass loss (removed from the outer mass shells of the MESA structure model) until the final helium star mass is reached. Because the lifetime in SSE may be an overprediction compared to the actual lifetime reached in MESA, we apply this mass loss between 80% and 90% of the predicted SSE lifetime so as to not reach the end of the MESA evolution before all [1] www.amusecode.org [2] This 'canonical' value of the solar metallicity may be an overestimate; see Asplund et al. (2009) for a review. of the required mass is stripped. The stellar evolution is then continued until the final lifetime found in MESA.
For our helium star models, the very final stage of evolution involves a rapid expansion of the remaining, tenuous envelope. Due to interaction with the close binary companion, this small amount of material in the envelope is in fact expected to be lost from the system; for our progenitors, we use models just prior to this stage, at which the outer radius of the helium star is still compact. For consistency, we evolve our (Z = 0.02, M 2 = 1 M ) companion star to the same age as that used for the primary star. In practice, this means the companion is still on the early main sequence, and therefore the effect of this small duration of stellar evolution on the structure and composition of the companion is negligible.

Hydrodynamical model set-up
We model the hydrodynamics of the SN using the SPH code Gadget-2 (Springel, 2005), running in the AMUSE framework. The SPH formalism has been shown to be effective in three-dimensional simulations of stellar phenomena (for example, Pakmor et al., 2012). One reason is that computational resources are not expended on regions of vacuum or negligible density (for example, Church et al., 2009;Lai et al., 1993), which constitute a significant fraction of the simulation volume in the current problem. Modelling a binary in a vacuum is easily handled in SPH, without the need for (low density) background fields in grid codes, which can exhibit artificial shocks from motions of other bodies within this background.
The Lagrangian nature of SPH describes advection naturally, without suffering from complications of numerical diffusion found in Eulerian codes, and we do not have to restrict the simulation to a fixed volume, which is useful in the present problem of a rapidly expanding shell of gas. A benefit to running the simulation in three dimensions is the absence of any boundary effects, which can produce on-axis artefacts (Marietta et al., 2000) or preferential wave numbers in the formation of instabilities (Warren and Blondin, 2013). As with all hydrodynamical codes, the SPH method also has its drawbacks, and some of these are discussed further in the context of our convergence studies in Section 2.3.1.
The stellar models created in MESA are converted into SPH particles using the star to sph routine in AMUSE, in a similar method to that outlined in de Vries et al. (2014). The routine first extracts the one-dimensional hydrostatic structure of the star, represented as a function of mass coordinates, from the data generated by the stellar evolution code. The SPH particles are initialised in a homogeneous sphere constructed from a face-centred cubic lattice, and the radial positions of the particles are then adjusted so as to match the density profile from the evolution code. [3] The internal energies of the particles are then assigned from the temperature (and mean molecular weight) distribution from the stellar evolution code. We use equalmass particles throughout these simulations (unequalmass particles can cause additional complications such as spurious mixing; Rasio and Lombardi, 1999).
The primary star is configured with a purely gravitational core particle of 1.4 M to model the neutron star. The softening length is chosen to be equal to the smoothing length, such that, due to the compact support of the cubic spline, the smoothing kernel reaches zero at 2.8 . This equality is also maintained for the SPH particles to preserve equal resolution of the gravitational and pressure forces. The zero-kinetic-energy models are relaxed over 2.5 dynamical timescales of the gas using critical damping on the velocities of the particles, where at each step the magnitude of damping is reduced so that in the final step no constraint on the velocity is imposed (for a similar approach, see de Vries et al., 2014). This is required due to effects of mapping the one-dimensional stellar structure on to the particle grid, and differences in physics between the codes, such as the value of the adiabatic exponent.
We set up the binary models at different orbital separations, a, where the minimum separation is chosen to be greater than the limit of Roche-lobe overflow (RLOF) of the companion star, given by the Eggleton (1983) relation, where R 2 is the companion radius and q is the binary mass ratio M 2 /M 1 . Once both stars have been initialised to the SPH code, orbital velocities are determined for a circular orbit at the specified separation and applied to each star.

Simulation of the supernova explosion
The SN is initiated using the 'thermal bomb' technique (Hirai et al., 2014;Young and Fryer, 2007), which assumes the core bounce has just occurred, at which moment we inject energy into a shell of particles around the neutron star. As discussed in Young and Fryer [3] Randomisation of the angular orientation of the particles has the undesirable effect of the additional shot noise it generates in the initial density distribution; however, the further step of damped relaxation used here will ultimately result in a glass-like configuration.
(2007), thermal bomb approaches (along with alternative, piston-driven shocks) are not intended to embody the physical mechanism that drives the SN. Indeed, the actual processes by which the energy gain occurs near the proto-neutron star are still not fully understood, though recent observations and insights from three-dimensional simulations have shed some light on the role of instabilities, asymmetries and jets in driving this process (Bruenn et al., 2013;Couch and O'Connor, 2014;Couch and Ott, 2015;Hanke et al., 2013;Janka, 2012;Lopez et al., 2013;Milisavljevic et al., 2013). The boundary of energy injection is specified by radius (which is equivalent to a fixed enclosed mass) instead of particle number. This allows scaling of the problem over a range of SPH particle numbers while keeping fixed the mass fraction that receives the SN energy. The total thermal energy (a canonical 10 51 erg ≡ 1 foe [4] ) is distributed evenly amongst these N SN particles, so that the specific internal energy per particle is increased by (1 foe) / (N SN m SPH ).
We found that a careful investigation of the effect of the injection radius was necessary. Too small a radius (and therefore N SN ) results in large, uncontrolled asymmetries in the shock front that grow from intrinsic small-scale asymmetries in the initial particle distribution. On the other hand, too large a radius results in the internal energy of the SN being distributed amongst a large number of particles, lowering the specific internal energy and therefore reducing the overall temperature in the region and weakening the shock. We found that, for the helium star models used here, injecting the SN energy into a region R SN 0.05R generates a sufficiently spherical shock while still keeping N SN small as possible.
Two of the parameters we wish to predict, calibrated from our simulations, are the final velocities (formally, at infinity) of the runaway components of distributed binaries. For masses m relative to the neutron star mass (i.e. m ≡ M/M NS ), TT98 predict these values in terms of the following initial parameters • a: the pre-SN binary orbital separation • v: the pre-SN relative orbital velocities • w: the magnitude of the kick applied to the NS • θ and φ: the spherical polar angles of the NS kick vector with respect to the 'x'-axis aligned along the NS orbital velocity vector at the moment of the kick • v im : the change in velocity of the companion due to the impact of the SN shell [4] More recently, the 'Bethe' (B) has been proposed as an equivalent unit in honour of Hans Bethe's work on supernovae (Weinberg, 2006;Woosley and Heger, 2007).
• v ej : the velocity of the ejecta shell • m 2 , m 2f and m shell : the initial mass of the companion, the final mass of the companion after mass loss, and the mass of material in the shell, respectively (all relative to the neutron star mass) In the original work of Wheeler et al. (1975), during the supernova shell passage over the companion star, the effect of mass stripping is parametrised by the fraction of companion radius x = R/R 2 as a function of angle around the star. Above some critical fraction of the companion radius, x crit , a fraction F strip of the mass is stripped by the shell impact, and below it a fraction F ablate of mass is ablated.
In this framework, the magnitude of the impact velocity v im is theoretically predicted to be Here, we use the prediction by Wheeler et al. (1975) in the form adopted by Tauris (2015), which applies the substitution (F strip + F ablate ) = F → F * = (F strip + F ablate ) α for some α, as well as a free parameter η to account for the fact that this tends to overpredict the value of v im . Effectively, η represents the final change in momentum of the companion as a fraction of the incident momentum in the shell. The fraction of mass stripped and ablated from the companion, F strip and F ablate , is calculated in Wheeler et al. (1975) using a polytropic star of index n = 3.
As noted in Wheeler et al. (1975), corrections must be applied to this formula as it neglects the presence of a rarefaction wave back through the ejecta, geometrical effects of curvature in the shell (more important for small a), inhomogeneities in the ejecta and radiative losses behind the shock. Further phenomena can also modify the effective impact velocity, such as the deformation of the companion by the shock passage (altering its cross-sectional area), the formation of a bow shock in the ejecta (during which time the flow deflects around the companion star), and shock convergence on the far side of the star (causing the asymmetric emission of material from this side of the star).
One key element in the predictions of TT98 is the estimate of the mass stripped from the companion. This is based on early work from Wheeler et al. (1975) and simulations of a planar slab of material hitting a star Fryxell and Arnett (1981), which have a low resolution by today's standards. As noted above, results based on plane-parallel ejecta profiles need correcting for the fact that the shell in reality has some curvature. Higher resolution simulations such as the present one provide a test of these early estimates, which are one of the sources of uncertainty in the results of TT98.
The stripping of mass in our simulations is measured by calculating the specific energy for each particle, which is negative for bound particles. The amount of bound mass in the companion is time-dependent over the course of the simulation due to energy exchange between particles. The stabilisation of mass bound in the companion determines the end of our simulation, which in practice occurs within 10 dynamical timescales of the companion star (∼ 2 × 10 4 s).

Convergence test
One feature of SPH that demands caution is that the resolution is dependent on the local density, and therefore the method loses resolution in the lower-density, uppermost layers of the stars in our simulations. In the current problem, the mass stripped by the secondary is from these same layers. Therefore, a good test for the resolution of the simulations is to look for convergence in the quantity of mass stripped from the companion. During the SN, Richtmeyer-Meshkov (the impulsive form of Rayleigh-Taylor) instabilities (RMI) are expected to be present, which have been found to appear once reverse shocks form at the interfaces between discontinuities in the density gradient (Kane et al., 1999). Such discontinuities are present in Type Ib/c progenitors at the interface between the carbonoxygen boundary in the core and, if any substantial fraction of hydrogen remains in the envelope, also at the helium-hydrogen boundary. However, these discontinuities tend to be smoothed during the conversion from the 1D stellar model and subsequent relaxation of the SPH particles. Proper treatment of the RMI requires a prescription of artificial conductivity that is not included in the current SPH codes in AMUSE. This instability is expected to be a significant factor in the mixing of stellar material early in the evolution of supernova remnants (SNRs), and so any evaluation of the fate of the composition of the SN ejecta must take this into account.
Unless the growth of RMI is explicitly seeded by some structure at the density interface, these instabilities will grow from perturbations at the numerical level of the simulation and may not, in such cases, grow substantially (Kane et al., 1999). Therefore, there is the potential for instabilities to be dependent on numerical effects such as the resolution of the simulation. Additionally, during the stripping of mass from the companion star, the initial deceleration of the shell Results for a convergence study using the amount of mass stripped from the companion. The primary star mass was 3 M and the orbital separation was 4 R in all cases. The number of SPH particles used in each run is given in the legend.
impacting the companion may be Rayleigh-Taylor unstable, but also that the subsequent flow of the shell over the surface of the star may induce some shearing (Kelvin-Helmholtz) instabilities (KHI). Due to the smoothing of discontinuities after relaxation of the SPH models, a lack of artificial conductivity [5] in Gadget-2 and the only perturbations being from noise in our SPH distribution, we expect that instabilities will not be fully captured in our simulations. Since the overall capturing of instabilities is reduced, so too should be the effect of instabilities on our results. Figure 1 shows a test of varying the SPH particle number, N , based on the amount of mass lost from the companion star (evaluated using equation 5). For low N , there is noticeable noise in the bound mass determination over time, but for N ≥ 10 5 particles, this is no longer appreciable. As shown in Figure 1, we did not find any substantial difference in the results increasing N from 5 × 10 5 to 8 × 10 5 . Accounting for this, as well as available computational resources, our simulations were run with 5 × 10 5 particles.  Table 1 Simulation initial conditions and main results. The first three columns indicate the initial conditions, where M 1 is the mass of the primary (helium star) and before the SN, M ej is the total ejecta mass, and a is the initial orbital separation. The last two columns are the amount of mass stripped from the companion star and the (magnitude of the) impact velocity.

Results
After reviewing the initial conditions used for our simulations, we examine the early stages of the SN (Section 3.1). We then investigate the magnitude of mass lost from the companion as a function of the orbital separation (Section 3.2), as well as the velocity imparted to the companion and the fraction of imparted momentum compared to the incident shell (Section 3.3). Next, we examine the newly formed SNR for asymmetries in morphology and metallicity (Section 3.4). Finally, we consider the subsequent evolution of a star altered by a SN shell impact (Section 3.5). Table 1 shows the initial conditions used in our simulations. The choice of primary and companion masses is motivated by the binary parameters used in TT98 and Tauris (2015), while the minimum orbital separations are chosen to be outside the RLOF value (equation 3). The final two columns show the effects on the companion due to the shell impact, discussed in more detail in the remainder of this section.

Shock breakout
By approximately 20 s after the SN is initiated, the forward shock has broken out of the surface of the helium star, during which time the fraction of SPH particles bound to the 1.4 M neutron star drops smoothly to almost zero. We find at late times that there is some fall-back of a small amount of material, which remains bound to the neutron star. As we do not model here the complexities of the magnetic field of the new neutron star or any form of pulsar wind, it is possible that other mechanisms later expel some or all of the residual bound gas.
In Figure 2, it can be seen that there is a rapid conversion of energy from internal (thermal) energy from the moment of explosion to kinetic and potential energy as the shock passes through the star and the subsequent shell expands. By approximately 100 s following the SN explosion, very little of the original thermal energy remains in the gas as it has been almost entirely converted into kinetic energy in the expanding shell. Figure 3 shows the changes in density through the 3 M helium star from shortly after core bounce until around the time the forward shock reaches the outer layers of the star. A lower density cavity with a very shallow gradient is seen to lag behind the expanding ejecta shell, which has a steepening density gradient as a function of radius towards the leading edge of the ejecta. This variation in gradient is maintained over time, although the overall magnitude of the density drops during the expansion. We investigate the density and velocity distributions within the ejecta in more detail in Section 3.3.

Impact and mass loss from the companion
The passing shell first strips material from the outer layers of the companion. The compression of the companion along the direction of motion of the shock causes heating of the stellar material, which results in a subsequent mass loss through ablation. This ablation of material is found to be a slower form of mass loss than the initial stripping phase. The passage of the shock through the companion can be seen in the density slices of Figure 4. Aside from the mass stripping from the sides of the star as predicted in Wheeler et al. (1975), there is also mass loss from the far side after the shock has passed through the star. The black vectors in this figure show the velocities for a random sample of all the SPH particles that were originally in the companion which have subsequently become unbound. These vectors have had the orbital velocity vector of the companion subtracted, and they are then projected onto the orbital plane. Because each SPH particle has the same mass, these vectors also indicate the relative momentum of the unbound particles.
Panel (b) of Figure 4 shows that, once the shock passes through the centre of the companion, it converges at the far the side of the star as it accelerates down the density gradient (similar shock convergence is seen around other spherically symmetric density gradients, such as in Rimoldi et al., 2015). This increases the local pressure on this axis, resulting in expulsion of material from the far side of the star and can counter the effect of the outward kick imparted by the incident shell of material (see also Marietta et al., 2000).
In the last panel of Figure 4, the central density of the companion has dropped and it has noticeably expanded from the shock heating. During this later The energy is broken down into kinetic (E k ), potential (Ep) and internal (thermal) (Et). This example corresponds to a primary of mass 4 M and an orbital separation of 4.5 R period (final three panels) ablation occurs for material whose specific internal energy is greater than the square of the sound speed. Due to the shock heating, the companion is 'puffed-up' like a pre-main-sequence star, and its luminosity is expected to increase temporarily as it reverts to thermal equilibrium (Marietta et al., 2000). Finally, we find a quadrupole oscillation of the companion that is induced by the distortion from compression due to the shock This ringing subsides after about one dynamical timescale of the companion star. Figure 5 shows an example of the variation in companion mass due to the shell impact. The stripping of mass by the passing shell causes a rapid mass loss in the initial phase, followed by a more gradual mass loss due to the later ablation of shock-heated material. The proportion of mass lost drops rapidly even by moderate orbital separations.
In Figure 6, we show the amount of mass lost from the companion (as a fraction of its initial mass) as a function of orbital separation. The lost mass is found by subtracting the final bound mass at the last snapshot of each simulation (which occurs at 2 × 10 4 s) from the initial mass. We use the last time possible from the simulation as the final mass takes much of the total simulation time to reach its steady-state value. The dashed line shows the prediction from TT98, the dotted line shows the fit obtained from Type Ia simulations compiled by Tauris (2015), and the dotdashed green line is from the new work of Liu et al. (2015). A least-squares regression gives a fit to our data of 1.3 (R/R ) −2.6 for the M ej = 1.6 M data and 0.58 (R/R ) −2.2 for the M ej = 2.6 M data. We find comparable values of the lost mass, in particular the steeper power-law gradients, to those seen in Liu et al. (2015).
There are some differences between our initial conditions and those from the previous work shown in Figure 6. Compared with the calibration from Type Ia simulations, our explosion energies and ejecta masses are both slightly different. Additionally, the 1 M companion radius, R 2 , shrinks slightly after relaxation of the SPH models compared with the canonical 1 R . The predictions in TT98 and Tauris (2015) depend on these quantities in particular within the geometric parameter Ψ, used originally by Wheeler et al. (1975), defined as This parameter is used in the determination of x crit as well as F strip and F ablate in Wheeler et al. (1975) using tabulated data for an n = 3 polytrope. For our comparisons, we adjust these quantities (and therefore Ψ) in the TT98 and Tauris (2015) estimates to match the initial conditions of our simulations. Furthermore, the simulations in Liu et al. (2015) also use a slightly different companion mass and ejecta mass, and so their results are not completely equivalent to ours.  (2015) is not fully equivalent, as both the ejecta mass and companion mass (and therefore radius) are slightly different.

Momentum transfer and the velocity of the companion
When the orbital separation is very small, the impact of the ejecta causes not only a significant loss of mass from the companion star but also a large change in velocity. The largest change in velocity of the companion occurs during the transfer of momentum from the shell in the initial impact. However, as the end of the shell passes over the far side of the companion, there is an overpressure acting on this side of the star when the shock converges on this axis. This causes the companion to receive a slight change in momentum in the direction opposite to the shell motion (which has been suggested in other simulations such as Fryxell and Arnett, 1981;Marietta et al., 2000). In the theory of TT98, v im is defined to be an effective velocity that not only accounts for the momentum imparted to the companion by the passage of the shell but also the subsequent change in momentum due to (potentially asymmetric) mass loss.
We found that measuring the velocity of the companion with respect to the compact remnant is complicated by the difficulty to define the baryonic centres of the binary system with the ejecta that had not yet left the binary system, the oscillatory behaviour of the Markers and line styles correspond to those used in Figure 6. Again, as for Figure 6, the comparison with Liu (2015) is not fully equivalent due to the slightly different ejecta mass and companion parameters.
companion star as a result of the shell impact, and the Brownian motion of the neutron star due to the shot-noise of the limited resolution it its vicinity.
As an alternative technique, we set up a co-rotating frame of reference that matches the original circular orbital motion. Before the SN and up until the shell impact, there is no component of velocity of the companion perpendicular to this direction of motion. However, during and after the impact, the companion (as well as the mass unbound from it) gains a component of velocity, and therefore momentum, in the radial direction with respect to this frame. We use this to measure the momentum delivered to the companion and the material removed from the companion, as well as the effective kick velocities. Figure 7 shows the component of momentum in the radial direction for material unbound from the companion that was not unbound in the previous time step. This allows us to see clearly the burst out of the back of the star; we measure the impact velocity just after this peak. The right panel shows the breakdown of momenta in the radial direction for unbound and bound material originally from the companion (and their sum). This gives an alternative indication of η, where we see that although less than half of the total incident momentum in the shell is delivered to this material in total, only 30% of the momentum is delivered to the (bound material of the) companion star.
Our final impact velocity magnitudes are shown in Figure 8. A least-squares regression gives a fit to these data of 556 (R/R ) −1.4 for the M ej = 1.6 M data and 652 (R/R ) −1.4 for the M ej = 2.6 M data. The velocities for both ejecta mass conditions follow a similar gradient to earlier work presented in TT98 and Tauris (2015), although it is not quite as steep as the −1.9 power-law of (Liu et al., 2015). The overall scaling differs from previous work, however. The early estimate from TT98 of the impact velocities used a value of η = 0.5 (see equation 4), whereas fits in Tauris (2015) and Liu et al. (2015) sit closer to η = 0.2. Our impact velocities lie in between these values, and we consider one possible cause in the Section 3.3.1.
Finally, we also consider the effect of drag from the remaining material on the companion velocity, noting that there is still a non-negligible density of gas interior to the ejecta shell. For a conservative estimate of this drag force from the innermost ejecta, we neglect any outward velocity of this gas, and take a density of 10 −3 g cm 3 in this material. With these values, the drag force on the companion will be for v 2 = 300 km s −1 , and where we approximate the drag coefficient of the star with a solid sphere value of C drag ≈ 0.5. For a companion mass of M 2 = 1 M the acceleration associated with this drag is therefore only 10 −5 km s −2 . Although small, drag induced by the lower-velocity ejecta may alter the final effective impact velocity when integrated over a long timescale.

Ejecta profiles
We investigate in more detail our ejecta profiles as a potential cause of the discrepancy between our impact velocities and those of Liu et al. (2015). Previous work, such as that of Liu et al. (2015), has often initialised the ejecta with the assumption that it is in a homologous expansion by the time it impacts the companion, so that, for a given t, v ∝ R. The density and velocity profiles in this ejecta are constructed from broken power-law fits to analytic treatments of the shock through the progenitor. These treatments have, in particular, been based on the polytropic envelopes (or one-dimensional structure models) of supergiant stars, and the power-law fits are to the (small and large R) asymptotic limits of a varying density gradient in the ejecta (Chevalier and Soker, 1989;Matzner and Mc-Kee, 1999).
In Figure 9, we show the variation of ejecta density and velocity as a function of radius from the center of mass (by averaging the SPH particles over concentric shells) and compare with an analytic function from the equations used in Liu et al. (2015). We also show, in the bottom panel of Figure 9, the distribution of velocity over mass in the ejecta. It is clear from Figure 9 that in the ejecta from our helium star models we have a shallower density gradient through much of the outer regions compared with the power-law profiles. In this outer ejecta, the velocity and density are also higher in our models. As the kick velocity has a strong dependence on this high-velocity ejecta (Liu et al., 2015), this can explain increased impact velocities seen in our simulations.
Finally, we note that we examined the ejecta for large-scale asymmetries by determining the shellaveraged radial profiles of density and velocity in hemispheres corresponding to the directions toward and away from the companion star. We found that the values in either direction agreed to within a few per cent, and therefore do not produce a discernible difference on the logarithmic plots in Figure 9.

Properties of the larger-scale supernova remnant
The amount of accretion on the companion has previously been shown to decrease with increasing shell velocity (Fryxell and Arnett, 1981); therefore, the high ejecta velocities in Type Ib/c supernovae lead us to expect little pollution of the companion. Indeed, we find negligible pollution of the companion star with SN material. However, the converse-pollution of the SNR with material from the companion-can be appreciable. A few 10 −2 M of hydrogen-rich material may be stripped from the companion by the passing shell in our simulations, which may be detectable as an asymmetry in the metallicity of the SNR on the side of the companion. Figure 10 shows a 3D rendering of the SNR and companion at 10 3 s after the moment of the SN. At this point, a hole has been created in the passing shell due to the presence of the companion, which is seen to persist at later times. The hole in the ejecta caused by the companion is approximately 30 degrees in size for the minimum orbital separations.
Even if an ejecta hole cannot be detected morphologically, the presence of a hole in SNR ejecta may allow the inference of a companion from a burst of UV radiation generated during the impact with the companion, which can escape through the less optically thick region of the companion shadow cone (Kasen, 2010). The hole may persist to late stages of the SNR despite some amount of refilling due to the subsequent rarefaction wave along with hydrodynamic instabilities (García-Senz et al., 2012;Kasen, 2010).
Not only do we observe a hole in the SNR due to the companion star, but we also see an increase in the density in a ring surrounding the hole, as shown most clearly in Figure 10. As shell material impacts the outer part of the companion star, where material is stripped and swept up with the ejecta, this ring of gas is also compressed in contrast with the freely expanding ejecta that do not interact with the companion. Aside from augmentation of the early SN light curve, our results also suggest that ring-like enhancements in density of the SNR could indicate the presence of a companion star. Ring-like structures may easier to detect than a hole in the SNR as the enhancement in density may also be associated with an increase in radiative losses in the ring. Figure 11 illustrates that the metallicity of the SNR is highest in the innermost regions, where the ejecta represents material nearest the core of the supernova progenitor. A hole develops in this high-metallicity ejecta, at first primarily due to the shadow of the companion star (panel a). At later times (panel b), the ablation of companion material further reduces the metallicity of a large fraction of the inner part of SNR in the direction of the companion. The orbital motion of the companion star within the inner SNR during this longer period of ablation can also enlarge the region over which the gas is enriched with hydrogen.
3.5 Post-impact evolution of the companion Following the stripping and ablation of mass from the outer layers of the companion star, we use AMUSE to model the stellar evolution of the companion and compare with an unperturbed stellar model evolving from the main-sequence. As we associate composition with each SPH particle from the original stellar model, we can convert the final SPH state of the companion back to a 1D structure model by an inversion of the method to construct the SPH model outlined in Section 2.2. After the 1D model is loaded back into MESA, we continue the stellar evolution and compare the final age and HR diagram tracks to an undisturbed companion model.
A 1 M star with metallicity Z = 0.02 evolves to the through to a carbon-oxygen WD at 12.1 Gyr in MESA. On the other hand, the 1 M model which has lost 0.04 M of material from the SN impact reaches this stage at a later age of 14.0 Gyr. Both the primary and companion stellar models were evolved for the same length of time (that is, the time it took for the primary star to reach its helium star stage).
Although the final age of the stars is noticeably different, there is not a large evolutionary difference as represented on the HR diagram of Figure 12. The luminosity of the stripped star is somewhat lower during its evolution, but in general this discrepancy is only around 20% of the value seen in the undisturbed companion model. It may therefore be difficult to distinguish a companion that has lost part of its envelope due to a supernova from T eff and L alone. However, the stripping and contamination in the outer layers of the star might produce differences in chemical abundances that are spectrally distinguishable from the coeval stellar population (see also, for example, Pan et al., 2012).

Discussion and conclusions
For supernovae in close binaries, the impact of the ejecta shell can have non-negligible effects on the mass and velocity of the companion star. This is therefore an important phenomenon in considerations of the runaway velocities from supernova-disrupted binaries and, for example, their potential to contaminate searches for hypervelocity stars from other origins, such as the Hills mechanism with the supermassive black hole in the Galactic Centre (Hills, 1988;Yu and Tremaine, 2003). We have performed SPH simulations of SNe in close binaries to study the consequences of the shell impact on the companion. The overall hydrodynamic phenomena and trends we observe during these simulations are broadly consistent with previous studies of Type Ia (Liu et al., 2012;Marietta et al., 2000;Pakmor et al., 2008;Pan et al., 2012), Type Ib/c (Liu et al., 2015) and Type II (Hirai et al., 2014) SNe. In addition we find that the gradient in the impact velocity predicted by Wheeler et al. (1975) matches our results well, with some modification of the η parameter representing the total momentum received by the companion.
While this work was in preparation, Liu et al. (2015) presented work on the effect of a Type Ib/c supernova shell impacting onto a companion star, in order to derive the linear momentum transfer between the shell and star. The velocity induced onto the companion due to the shell impact in their work is a factor of 1.5 ∼ 2 lower than our results. Although it is not straightforward to separate the causes of these discrepancies, there are a number of differences between our calculations. Most notably, in our case, the shell is naturally formed from the supernova explosion mechanism, that originated in the exploding star, as opposed to the introduction of the supernova shell by an analytic description of the material that impacts the companion.
Using the calibration from our simulations, we return to the question of runaway velocities of the components of supernova-disrupted binaries considered in TT98 and Tauris (2015). We have created a python code that calculates the final speeds derived by TT98 in order to investigate the analytic predictions with our simulation results. In this Monte Carlo code, an impulsive increase in velocity, w, is imposed to the neutron star, randomly oriented from an isotropic distribution over a sphere. This is achieved by mapping from a uniform random distribution over t ∈ (0, 1] to 2πt for the angle φ, and from a uniform random distribution over u ∈ [0, 1] to cos −1 (2u − 1) for the angle θ. Figure 13 shows a comparison of the distributions of speeds with (red) and without (blue) the effect of applying v im and mass loss in the companion star.
From Figure 13, it is evident that, although adding an impact velocity to the companion (perpendicular to its orbital velocity) increases the minimum companion speed, it also in fact reduces the maximum companion speed. To clarify the discrepancies in the distributions that occur when adding v im , we consider the effect of NS kick angles on the final velocity of the companion star in disrupted binaries in Figure 14, analogous to Figure 4 in Tauris (2015). The white regions for high θ in each panel represent binaries that remain bound after the NS kick (and thus the runaway velocity is undefined). The grey regions represent NS kick angles for which the NS and companion star merge after the SN. It can be seen from the lower panel that the effect of applying an impact velocity to the companion star can stabilise the systems where the NS kick is counteraligned with the NS orbital velocity. In fact, the small  Figure 13 Comparison of Monte Carlo sampling over NS kick orientation for the case with no impact effects on the companion (blue) and with impact effects as calibrated from our simulations (red). Both cases show 10 4 samples of NS kick orientation with a uniform on distribution over the unit sphere. The magnitude of the NS kick velocity, w, is fixed at 800 km s −1 throughout. Mean values of the runaway velocities, as well as percentages of cases where the binary components remain bound and merged (and, in parentheses, merged cases that were calculated as bound). Neither bound nor merged cases appear in these distributions. region of parameter space giving large values of v 2 at φ = 0 and high θ is removed after adding v im (due to these systems now remaining bound). This explains the potentially counter-intuitive result that by adding an additional velocity to the companion star in fact reduces the maximum possible velocity of runaway stars.
The main results of our work are as follows: 1 We follow the supernova self-consistently from just after the core bounce in a helium star generated from a stellar evolution model. Exploding a SN in such a model produces an ejecta profile that is different from that used in previous work, which has employed an analytic function for the ejecta distribution derived from the theory of shocks travelling through a one-dimensional atmosphere. The progenitor model used here is still somewhat artificial in construction, with a constant mass loss parameter late in its evolution. As understanding of Type Ib/c progenitors improves (for some recent investigations, see Kim et al., 2015), future work would benefit from a more realistic progenitor model by modelling the mass loss processes in detail. 2 We have investigated the mass removed from the companion in the very close binary separations seen in Type Ib/c SNe, as well as the net change in momentum of the companion star due to the shell impact and later ablation of material. The change in momentum of the companion is used to calibrate theoretical predictions of the final velocity components in dissociated binaries, which is useful for studies of runaway stars from SNe, as seen in the recent work of Tauris (2015). 3 We investigated the morphology of the SNR shortly after the shell has passed the companion, as well as the pollution of the SNR with material stripped from the companion, which, for the case of Type Ib/c SNe, may be a detectable fraction of the total mass in the ejecta (several 10 −2 M out of ∼ 2M ). The metallicity of the SNR is found to be highest in the inner regions of the SNR, behind the expanding shell at the shock front, and in this region the ablation of hydrogen from the outer layers of the companion star can dilute the metallicity on the side of the SNR facing the companion, resulting in a strong asymmetry in metallicity in the orbital plane. 4 The companion star is additionally found to modify the morphology of the SNR in two distinct ways: as anticipated, a hole forms in the SNR on the side of the companion; also, an increase in the SNR density is seen in a ring around the hole, which may enhance the luminosity in SNR observations. 5 We have also followed the stellar evolution of the companion star after the removal of mass during the shell impact, and have shown that the (mainsequence) lifetime of a star that has suffered a supernova impact can be substantially lengthened while traversing a slightly different track on the HR diagram. v t (c) Figure 9 Blue lines show the profiles of (a) density, (b) velocity and (c) the mass distribution of velocity within the ejecta for our simulations (binned in radial shells form the center of mass of the ejecta). Black dashed lines show the power-law profiles used in Liu (2015). Both cases were calculated for a time of 90 s after the supernova. The transition velocity (and radius at which this occurs) in the Liu (2015) profiles are given as dotted lines.
(a) (b) Figure 10 A 3D rendering of the gas density in the system 10 3 s after the SN, viewed down the x-axis (a; the original axis of the binary) and y-axis (b). The companion is still distorted due to the impact, and has produced a hole in the expanding ejecta. We used the software Mayavi2 (Ramachandran and Varoquaux, 2011) for the visualisation. (a) (b) Figure 14 Distribution of NS kick angles from the x-axis (aligned with the pre-kick NS orbital vector) for the same parameters as the Monte Carlo runs shown in Figure 13 (but with 2 × 10 5 samples to better fill the space in angles). The span of θ over φ = 0 defines the orbital plane, where θ = 0 is for a kick aligned with the NS orbital velocity vector. Panel (a) shows the case of no impact effects on the companion (blue distribution in Figure 13) and panel (b) shows the case where impact effects calibrated from our simulations are included (red distribution in Figure 13). Colours represent the magnitude of the companion star from disrupted binaries (as a fraction of the maximum runaway velocity). Grey shows cases where the NS and companion star merge. The white regions for large θ are cases where the binary remains bound and so there is no runaway companion.