Explosions in Roche-lobe Distorted Stars: Relativistic Bullets in Binaries

State-of-the-art surveys reveal that most massive stars in the Universe evolve in close binaries. Massive stars in such systems are expected to develop aspherical envelopes due to tidal interactions and/or rotational effects. Recently, it was shown that point explosions in oblate stars can produce relativistic equatorial ring-like outflows. Moreover, since stripped-envelope stars in binaries can expand enough to fill their Roche lobes anew, it is likely that these stars die with a greater degree of asphericity than the oblate spheroid geometry previously studied. We investigate the effects of this asymmetry by studying the gas dynamics of axisymmetric point explosions in stars in various stages of filling their Roche lobes. We find that point explosions in these pear-shaped stars produce transrelativistic ejecta that coalesce into bullets pointed both toward and away from the binary companion. We present this result and comment on key morphological differences between core-collapse explosions in spherical versus distorted stars in binary systems, effects on gravitational wave sources, and observational signatures that could be used to glean these explosion geometries from current and future surveys.


INTRODUCTION
Massive stars predominantly form in binary and multiple systems and the majority of them will interact and exchange mass with at least one companion star (Sana et al. 2012;Moe & Di Stefano 2017).Because these stars are so massive, they will die as core-collapse supernovae (CCSNe), explosions that synthesize many of the heavy elements and drive outflows from galaxies (e.g.Fowler & Hoyle 1964;Woosley et al. 2002;Efstathiou 2000).
The shape of a star in a binary system can deviate significantly from spherical symmetry.This is particularly true when a star expands to fill its Roche lobe.This forces the star to become "pear-shaped" due to the combined effect of the gravitational pull of the companion and the centrifugal effect from the motion around the common center of mass.In this paper, we explore what happens if a supernova explosion occurs in such a highly asymmetric, Roche-lobe-deformed star.
Explosions in non-spherical stars are of interest because of the effects that an oblique shock can have when breaking out of the stellar surface.The non-radial motions can dramatically alter the dynamics and emission of shock emergence (Matzner et al. 2013).In a previ-ous paper, DuPont et al. (2022) numerically explored the effect of point explosions in rapidly rotating, oblate progenitors and showed how these produce relativistic equatorial ring-like outflows.Here we build upon this work and extend this to the pear-shaped stars in binary systems.
The fraction of core-collapse supernovae progenitors that are significantly Roche-lobe deformed is not very well known, but binary evolutionary simulations suggest that this may be common for lower mass progenitors of stripped-envelope supernovae (type Ib/c and type IIb, Filippenko 1997).These are produced in binary systems where the progenitor loses most of its H-rich envelope during a first mass transfer phase (Kippenhahn & Weigert 1967;Götberg et al. 2017).The resulting stripped star first contracts as a hot helium star.Drout et al. (2023) and Götberg et al. (2023) recently presented the discovery and analysis of such stars, which so far had been missing.
During their later evolution, helium stars in the mass range 2 − 4M ⊙ can expand significantly to become helium giants (e.g., Divine 1965;Habets 1986;Dewi et al. 2002;Dewi & Pols 2003a;Yoon et al. 2010;Sravan et al. 2019).If they retain a layer of hydrogen, the expansion is even more drastic (Laplace et al. 2020;Klencki et al. 2022).This implies that these low-mass stripped progenitors are (close to) filling their Roche lobes at the onset of core collapse and therefore deformed.
Of special interest are supernova explosions in Rochelobe deformed stripped stars where the companion star is a neutron star.These systems may lead to ultrastripped SNe (Tauris et al. 2015).This is the last stage in the pathway towards the formation of double neutron stars and gravitational wave sources (Tauris et al. 2017).In this work, we present a suite of two-dimensional axisymmetric simulations of point explosions in a 3M ⊙ pre-supernova helium star at various stages of filling its Roche lobe, described in Section 2. We study how the shock propagates and breaks out from the surface.We report on the formation of a trans-relativistic ejecta pointed both away and towards the companion in Section 3. We discuss the astrophysical implications and observational consequences in Section 4. Section 5 provides a summary with conclusions.

Governing Equations
The governing hydrodynamical equations in the special relativistic regime describe the conservation of baryon number, and conservation of stress-energy, where we work in the fluid frame so that T µν = ρhu µ u ν +pη µν .In these equations, h = 1+ϵ+p/ρ is the total specific enthalpy, ϵ is the specific internal energy, ρ is proper mass density, η µν is the Minkowski metric, p is pressure, is Lorentz factor, and β is the fluid velocity in units of c.The signature of the metric is (− + ++).The system of equations is closed by assuming an adiabatic equation of state p = (γ − 1)ρϵ, where γ = 4/3 is the adiabatic gas index.We work in units where c = 1.The point explosions are launched in a 2D axisymmetric spherical-polar mesh in the (r−ϕ) plane.All quantities are made dimensionless by combinations of R ⊙ , M ⊙ , and c.The code used to run the simulations, entitled SIMBI, was written by this study's first author, and is an open-source, second-order, Godunov-type scheme.SIMBI is is parallelized for multi-CPU frameworks enabled by OpenMP, and supports GPU-acceleration on both AMD and Nvidia platforms (DuPont 2023).

Initial Conditions
Using the MESA (Paxton et al. 2011, 2013, 2015, 2018, 2019) code, Farmer et al. (2021) evolved a series of binary stars with a focus on the primary star post helium burning.The authors followed the evolution of said primary stars just before core collapse.The binary system had various initial masses, but constant initial mass ratio q 0 = M 2init /M init = 0.8, and each star was nonrotating.We use these stars as the progenitors in this study.From knowing the initial mass, M init , final mass, M 1 , and initial mass ratio, q 0 , we can arrive at the final mass ratio of the system pre-supernova assuming the mass exchange was conservative, This parameter is crucial for constructing the equipotential surfaces of the Roche lobe.The separation of the stars can be calculated by inverting the Eggleton (1983) formula to arrive at where R eq is the volume equivalent radius defined such that the volume of the Roche lobe is V = 4/3πR 3 eq .For simplicity, we choose a single primary star from the grid of binary simulations for our entire analysis.This primary had an initial mass of M init = 11M ⊙ and presupernova mass of M 1 = 3M ⊙ , which gives a final mass ratio of q = 5.21.Its volume equivalent radius is roughly 5.4R ⊙ .
In this analysis, we assume a steady wind environment, where Ṁ is the stellar mass loss rate, v wind is the velocity of the wind, and R * is the outer radius of the distorted envelope.A common convention for modulating stellar wind magnitudes is the mass-loading parameter A = A * × 5 × 10 11 g cm −1 , which assumes Ṁ = 10 −5 M ⊙ yr −1 and v wind = 10 8 cm s −1 .We fix A * = 1 in our simulations.
The MESA density profiles are computed in 1D using the volume equivalent radius, R eq .Since the equipotential surfaces of the Roche lobe are themselves isobars, we can invoke hydrostatic equilibrium (−∇P/∇Φ = ρ) to map the density ρ(R eq ) to the surfaces of constant potential, Φ.We also assume synchronous rotation for the system such that the dimensionless potential has the form (Kopal 1959): where we set θ = π/2 so that r − ϕ lies in the orbital plane of the binary.The formulae above allow us to transform the 1D density profiles from ρ(R eq ) to multidimensional profiles ρ(Φ), conserving the stellar mass in the process.To map the density contours for a Roche lobe filling factor, F, we use the functional form of the modified potential described in Leahy & Leahy (2015): where Φ L1 is the potential at the L1 Lagrange point of the system.We model the supernova explosion as a "thermal bomb" in a small region of the grid, where δr is the radial extent of the overpressure, and H(r) is the Heaviside step function.Since we are constrained geometrically by the detailed evolution of the binary, we are only afforded the liberty of varying E/M and the Roche lobe filling factor, F. Results are presented with 2048 radial zones per decade on a logarithmically-spaced spherical-polar mesh with ∆r = r sin(π/2)∆ϕ.The radial extent is r min = 10 −3 R ⊙ to r max = 10R ⊙ and the azimuthal extent covers half the sphere, ϕ ∈ [0, π].This results in a grid of 8192 radial zones by 2795 angular zones.In other words, the explosion is resolved in 361 radial zones, and the star is resolved within 7644 radial zones.We set δr = 1.5r min .

RESULTS
Figure 1 depicts the shock wave evolution from a 10 51 erg explosion set off inside the 3M ⊙ progenitor.The first snapshot of the explosion shows that at first, the supernova shock wave propagates outward quasispherically while still interior to the stellar mantle in panel (a).This is because the iso-density contours near the stellar core are quasi-spherical due to gravity, so there is no preferential density gradient in either direction that would initially accelerate and deform the supernova shock (Sakurai 1960;Matzner & McKee 1999).However, because of the geometry of the Roche lobe filling star, the shock waves will break out at the stellar poles first, i.e., at θ = 0 and not ϕ = 0, since the ratio of the distance to one of the stellar poles to the distance to the back of the star is R z /R back ≈ 0.85.The result of this asynchronous breakout is shown in panel (b) of Figure 1 where the shock continues its outward expansion, but the oblique breakout causes the shocks to wrap around towards the nose of the Roche lobe as well as the back of the progenitor.The wrap-around effect is caused by the naked blast wave realizing a lateral pressure gradient in the φ direction resulting in the blast wave forming a pattern wave along the stellar surface, pointing towards impending collision points along the center-ofmass axis of the binary (see e.g., Matzner et al. 2013;Afsariardchi & Matzner 2018;Scully et al. 2023, and references therein).Finally, in panel (c) we see the aftermath of the lobe shocks colliding toroidally at a point located along the center-of-mass line of the binary system.The shocked lobes coalesce at this point, forming a hot, pressurized blob of ejecta whose secondary acceleration is caused by the steep pressure gradient directed along the binary center-of-mass axis.This ultimately results in a violent spray of ejecta that form collimated relativistic bullets pointing both away and towards the companion.The front bullet has a mass of 6 × 10 27 g and the back bullet has a mass of 5 × 10 25 g, which are of order the masses of the Earth and Moon, respectively.
To quantify some of the observable properties of the relativistic ejecta, a useful measure is the cumulative kinetic energy distribution above some four-velocity, E k (> Γβ). Figure 2 plots E k (> Γβ) for all geometries and explosion energies considered in this work.There is non-monotonic behavior in the Roche lobe stars where the geometries F ≥ 0.98 extend the high-velocity tail beyond the ballistic velocities shown for their spherical counterparts.The F = 0.95 stars, however, produce no material in such high-velocity ejecta and are slower overall than the spherical counterparts.Although the F ≥ 0.98 explosions extend the high-velocity tail of the kinetic energy distribution by factors of only a few, the angular distribution of this fast tail shows a clear split in the azimuthal angle.
The solid angle of the bullets at the different velocity cutoffs is from where E iso is the isotropic-equivalent energy weighted by the energy per solid angle.With Equation 9, we compute E k (>{0.5, 1})| ϕ=0 • = {2 × 10 44 erg, 5 × 10 41 erg}, and E k (>0.5)| ϕ=180 • = 6 × 10 43 erg as the total kinetic energies pointed towards and away from the companion for the 10 51 erg explosion.For the 10 52 erg explosion, the total energies are E k (>{0.5, 1})| ϕ=0 • = {4 × 10 45 erg, 2 × 10 45 erg}.Lastly, Figure 3 shows the separation between 10 51 erg explosions in the canonical spherical progenitor and a star that has filled its Roche lobe.At freeze out of the acceleration in both cases, the Roche lobe filling star produces trans-relativistic ejecta that gets focused along the binary separation axis while there is no such fast ejecta in the perfectly spherical supernova explosion.This distinctive asymmetry in the velocity and density profiles might lend itself to various hints at the supernova signatures observed and are further discussed in the next sections.We remind the reader that these results are idealized to 2D axisymmetry, which may enhance the effect of the shock converging at the nose and backside of the Roche-lobe.Exploration of the 3D nature of this problem will be needed and we plan to undertake this in a future study.

A relativistic bullet aimed at secondary and implications for shock breakout
The late-stage morphology of explosion in the Rochelobe distorted star implies that the companion lies along the path of one of the relativistic bullets.For the companion star to dodge the relativistic bullet, its sectorial orbit would need to sweep an angle greater than the sum of the solid angle of the companion with that of the bullet which leads to the companion radial constraint: where R C is the radius of the companion, d is the binary separation, ∆ϕ orb ≡ tan −1 (V orb /β bullet c) is the sectorial width of the companion orbit in the time it takes the relativistic bullet to reach the companion's original position, V orb ≈ (GM/d) 1/2 is the orbital velocity of the companion, β bullet is the velocity of the ballistic bullet, and ∆ϕ bullet is the half-opening angle of the bullet.With the bullet moving at the speed of light and having ∆ϕ bullet ≃ 0.01, Inequality 10 requires V orb > 0.01c in order for the argument in the tangent function to be non-negative.Since typical orbital velocities in binary systems don not exceed a few hundred kilometers per second, this means that the bullet will hit any companion regardless of its type whether it is a neutron star (NS), a black hole (BH), or a main sequence (MS) star.We reserve discussion of compact companions for later sections since their predicted rates are lower than those of MS companions (see e.g., Zapartas et al. 2017;Renzo et al. 2019).
While the energy budget of the directed bullet is not enough to unbind typical envelopes for massive companion stars, the deposition of energy could potentially affect the companion's evolution through induced mixing effects, seismic activity, or the collision could release a burst of radiation visible from the shadow cone (Kasen 2010;Ogata et al. 2021).Furthermore, the companion star will be shock heated, expand, and radiate, which Figure 3. Upper : The cumulative kinetic energy per unit solid angle for the 10 51 erg explosion as a function of azimuthal angle, ϕ, in the star that has filled its Roche lobe, i.e., F = 1.The inset marks the region aimed at the companion where the fastest ejecta is focused, carrying total energy E k (>1) = Ω bullet dE k /dΩ = 5×10 41 erg towards the companion in a bullet of angular width ∼ 1 • .There is also a counter bullet pointed away the companion carrying away total energy ∼ 10 44 erg in ejecta with Γβ > 0.5.Lower : The cumulative kinetic energy per unit solid angle for the 10 52 erg explosion where most of the focusing happens towards the companion with E k (>1) = Ω bullet dE k /dΩ = 2 × 10 45 erg.will produce late-time radiation signatures at the supernova location (see e.g., Maund et al. 2004;Folatelli et al. 2014;Maund et al. 2016;Ryder et al. 2018;Maund 2019;Sun et al. 2020, for up-to-date observations of strippedenvelope SNe with companion detections), but since the energy imparted by the bullet will be small, the companion's inflation timescale will be shorter lived than values expected for spherically symmetric SN hitting their companions.We also note that if the explosion energy is large enough, the ejecta-ejecta collision along the binary axis could produce a flash of photons before the SN light is seen as discussed in Scully et al. (2023), and the combined energies of the bullets in the 10 52 erg explosion simulated in this work are ∼ 10 46 erg, roughly the released energy of XT 080901 (Soderberg et al. 2008).Moreover, the stripped-envelope progenitors studied in this work are possible candidates for rapidly evolving SNe now better known as a subclass of fast blue optical transients (FBOTs; Drout et al. 2014;Kleiser & Kasen 2014;Kleiser et al. 2018b,a), though the luminosity peak will shift since the breakout flash is partially obscured by the obliquity of the breakout (Matzner et al. 2013).All of these together -the ejecta-ejecta flash, short shock cooling timescale of the companion, and non-standard light curve morphology -might provide hints towards classification of asymmetric SNe in binaries versus their spherically symmetric counterparts.

Polarization signatures
In principle, polarization measurements can be made if there is asymmetry in the supernova progenitor itself or if the explosion engine was asymmetric at the outset (Höflich 1991;Yamada & Sato 1991;Höflich et al. 1996).In the case for the Roche lobe bullets discussed in this work, the asymmetric outflow will likely give rise to measurable non-zero net polarization, assuming that it is dominated by Thompson scattering on free electrons.To get a sense of the lifetime of the asymmetry in the flow, we estimate the deceleration of the backwards bullet based on results from our simulations and assumptions about the stellar environment.The front bullet will be quickly destroyed by interactions with the companion, but the backwards bullet will decelerate after a time M bullet 4πA /βc, which, when scaled for a few reference parameters from the 10 51 erg explosion, gives and we have chosen A * = 10 −3 such that asymmetry of the outflow lasts ∼ week which is roughly the earliest cadence at which we can currently observe SNe post explosion (see Wang & Wheeler 2008;Patat 2017, and references therin).Equation 11 is two-pronged.That is to say, the stellar environment surrounding the Roche lobe explosion would have to be quite rarefied or the mass of the bullet (and therefore the explosion energy) would have to be larger in order to measure polarization signatures for these types of explosions.At least one of these conditions can be satisfied based on: (i) the large upper limits placed on the explosion energies of corecollapse events; (ii) the stellar environment surrounding stripped-envelope stars being uncertain.Furthermore, if the asymmetry lasts much longer than we estimate, it is possible that bullets in binaries could produce corecollapse SN remnants with "ears" (Grichener & Soker 2017).

Asymmetric engines
Although this work is in the context of spherically symmetric explosions in aspherical progenitors, there lies a plethora of observational evidence that the resultant NS star formed from core-collapse can have a natal kick from an asymmetric explosion mechanism (see e.g., Hobbs et al. 2005a;Scheck et al. 2006;Wongwathanarat et al. 2013;Janka 2017).An asymmetric explosion in an asymmetric progenitor will likely greatly change the distribution of energy over solid angle shown in Figure 3 since the resultant bullet velocities strongly depend on the time and the angle at which they collide along the binary axis.The advantage of the point explosion is the exterior lobe shocks colliding simultaneously with the interior SN shock that breaks out at the Roche lobe nose.This advantage of simultaneous collisions depends on the explosion energy and F as shown in Figure 2. When the explosion itself is asymmetric, the dynamics of the collision outside of the Roche lobe filling star will change in a non-trivial way and should be studied in a future work.

Fast Radio Bursts
In the scenario in which the binary companion is a neutron star (NS), the relativistic bullet aimed toward the companion will encounter the strong magnetic field of the NS magnetosphere.In this case, as the plasma bullet hits the magnetic field and is decelerated, some of its kinetic energy can be emitted in the form of radio waves.The mildly relativistic bullet has energy ∼ 10 44 erg and initial radius of ∼ 3 × 10 7 cm.Assuming the bullet is stopped almost instantaneously by the strong NS magnetic field, then the time for the reverse shock to traverse the bullet is roughly 10 −3 s, of order the duration of a fast radio burst (FRBs; Lorimer et al. 2007).We use the plasma frequency, ν, as an estimate for the radio wave frequency, where n e = ρ/m p is the electron number density in the bullet, m p is the proton mass, e is the elementary charge, m e the electron mass, r e is the classical electron radius, and we've accounted for the ratio of solid angles of the bullet with that of the magnetosphere of the NS where R SO is the stand-off radius with value: where B d is the surface dipole field of the NS and we follow the convention that quantities Q x = Q/10 x with cgs units implied.From this, we find This implies that the bullet meets the energy, duration, and frequency requirement for an FRB if it crashes into an NS magnetosphere, assuming an efficiency 10 −5 −10  (Zapartas et al. 2017;Renzo et al. 2019), we predict that supernova explosions in close binaries with a magnetized NS companion are a viable source for nonrepeating FRBs via this mechanism.In this scenario, the FRB would be followed by the rising light curve of a stripped-envelope supernova (SN Type Ib/c) several days later.

Relevance for the formation of Double NS and GW sources
Explosions of Roche-Lobe distorted stars are expected in the main binary formation scenarios for double compact objects.This is especially true in the main formation pathway for double neutron stars (e.g.Dewi & Pols 2003b;Tauris et al. 2017), where both progenitors may be significantly distorted.Our findings show potentially observable signatures that could be used to search for systems evolving through this formation channel, see for example De et al. (2018).
The first explosion is typically coming from the star that was once the most massive star in the system.Depending on the separation and the amount of swelling of the progenitor, this may already be an example of a Roche-distorted star with an main sequence companion.Observable consequences may result from ejecta hitting the companion, which is most likely still a main sequence star in this phase (Sect.4.1) and the aspherical outflows which may give rise to polarization signatures (Sect.4.2).
Most systems unbind during the first explosion, due the supernova kick resulting from asymmetries in the explosion (e.g.Hobbs et al. 2005b).Although we show in this work that additional asymmetries may result from the deformation of the star, we find that these will not play any significant role for the momentum imparted onto the system and the newly born neutron star.
Systems that remain bound can evolve further as a binary system.The second star now evolves to fill its Roche-lobe.It first loses its hydrogen-rich envelope in a highly non conservative event during which the orbit shrinks.The result is a compact helium star orbiting a neutron star in a tight orbit.The helium star is expected to swell once its start burning helium in a shell around its core (often referred to as case BB mass transfer).It will fill the Roche lobe anew now feeding gas to the neutron star, a stage that may be observable as an X-ray binary.Eventually the helium star will complete the advanced burning stages it will explode too, likely while still filling its Roche-lobe.At this stage we expect a heavily Roche-lobe distorted star with a neutron star companion, possibly with the observational consequence of a fast radio burst (Sect.4.4) accompanied by a supernova.

SUMMARY AND CONCLUSIONS
This work has demonstrated the following: (1) pointlike explosions in pear-shaped stripped-envelope stars cause shocks to converge to a point, producing relativistic bullets pointed away and towards their companions even for the typically non-relativistic 10 51 erg explosions; (2) such explosions only impart ∼ 10 44 − 10 45 erg onto their companions, drastically reducing the shock cooling timescales for main sequence companions assumed for spherical SN in binaries; (3) these explosions can obscure the breakout flash and therefore the luminosity peak as seen for many rapidly evolving SNe; (4) pointlike explosions in pear-shaped stripped-envelope stars might be candidates for Fast Radio Bursts associated with a core-collape SNe if the relativistic bullet interacts with the magnetic field of a neutron star companion.In short, as the vampire companion star depletes the primary of its envelope, at the end of the primary's life it responds by blasting the vampire star with a relativistic bullet.
This work was supported in part by NASA ATP grant 80NSSC22K0822.MD and SdM thank Rüdiger Pakmor, Mathieu Renzo, Shazrene Mohamed, Nando Patat, Rob Farmer, Stephen Justham, Eva Laplace, and Norbert Langer for useful discussions and comments during the MPA-Kavli summer school 2023, and MD acknowledges useful discussions with Andrei Gruzinov and support from LSSTC Data Science Fellowship Program, which is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation, and the Moore Foundation.(Paxton et al. 2011, 2013, 2015, 2018, 2019), astropy (Astropy Collaboration et al. 2013)

Figure 1 .
Figure1.The dynamical evolution of a 10 51 erg explosion within a Roche-lobe filling stripped star (F = 1, shown as the grey dashed line with volume-equivalent radius ∼ 3.8 × 10 11 cm).Logarithmic pressure is shown on the right half and four-velocity on the left half.In panel (a), the shock wave evolves quasi-spherically while still inside the stellar mantle.In panel (b), the supernova shock has broken out of the stellar poles first and begins wrapping around the star towards collision points at the Roche-lobe nose and the back of the progenitor.In panel (c), the lobe shocks collide and force material along the center-of-mass axis of the binary, forming two relativistic bullets directed towards and away from the companion with masses 6 × 10 27 g and 5 × 10 25 g, of order the masses of the Earth and Moon, respectively.

Figure 2 .
Figure2.The ballistic cumulative kinetic energy distribution for the 10 51 erg and 10 52 erg explosions at 162 s and 58 s, respectively, for all geometries considered.The progenitors with Roche-lobe filling factors F ≥ 0.98 extend the highvelocity tail of the ejecta by factors of a few when compared with their spherical counterparts.In general, the F = 0.95 produces the least amount of relativistic material.

Figure 4 .
Figure 4. Evolution of the 10 51 erg supernova a few minutes post explosion in the stripped star (a) and in the same stripped star that has filled its Roche lobe (b).Lab frame density is plotted on the right side and four-velocity is plotted on the left.The grey-dashed curves mark the original outer surface of the 3M⊙ progenitor with volume-equivalent radius 3.8 × 10 11 cm.The outflow velocity for the spherically symmetric star peaks at around Γβ = 0.5 while the fastest ejecta is focused along the binary axis for the Roche lobe filling star and has velocity Γβ ≳ 1.
CMasher (van der Velden 2020), SIMBI (DuPont 2023), MESA −4 for converting the kinetic energy into the isotropic radio emission of 10 39 erg − 10 40 erg quoted for many FRBs (see Cordes & Chatterjee 2019, for a detailed review of the current state of FRB science).Based on recently calculated rates of SNe in the universe (∼ 1 s −1 Li et al. 2011) and FRB rates (∼ 0.1 s −1 Thornton et al. 2013) being in rough agreement with the rate of NS companions existing in close binaries with massive stars