Gravitational waves from inspiraling binary black holes (BBHs) provide insights into the lives and deaths of massive stars. Population synthesis allows us to model these binaries through isolated binary evolution, but its predictive power is limited by difficulties in varying the stellar models and their associated uncertainties. We present a new grid of stellar tracks computed with the open-source stellar evolution code MESA, spanning metallicities 10^-3≤ Z/Z_⊙≤ 7. We vary two stellar physics parameters: wind-driven mass loss and the convective boundary mixing (CBM) mechanism. We pair these models with the Method of Interpolation for Single Stellar Evolution (METISSE) and binary population synthesis code COSMIC to obtain synthetic populations of merging BBHs in the local Universe. We find a maximum in the primary mass spectrum near 10M_⊙ which in most model variations is composed of two sub-populations at ≈8M_⊙ and ≈13 M_⊙, with the higher-mass population dominated by BBHs whose progenitors underwent a mass ratio reversal (MRR). This population also suggests an anticorrelation between higher primary masses and mass ratio, as BBHs with m_1⪆10M_⊙ preferentially undergo MRR and prefer a final mass ratio of q≈0.7. However, the location and relative strength of these two sub-populations is sensitive to our assumed stellar physics: varying both the wind and CBM treatments can merge the MRR and non-MRR populations into a single peak near 9M_⊙. Variations in our stellar tracks, especially CBM, lead to a factor of ≈6 difference in the rate, primarily due to modulation of the common envelope formation channel.

The origin of the large-scale poloidal magnetic field required to power relativistic jets in collapsars remains uncertain. While such a field may be inherited during proto-neutron-star collapse, the efficiency of this process is unclear, motivating an in situ mechanism to generate poloidal fields out of the predominantly toroidal fields produced by stellar differential rotation. We present the first 3D general-relativistic magnetohydrodynamics collapsar simulations initialized with toroidal magnetic field profiles that closely follow those of precollapse stellar models. As the toroidal field in the disk becomes dynamically important, it seeds the dynamo, producing coherent poloidal magnetic loops that appear at $\sim { \mathcal O }(100)$ gravitational radii and are then advected inward along paths that may deviate from the disk midplane. The resulting poloidal fields thread the black hole (BH) and launch highly variable, wobbling relativistic jets on timescales of order seconds, with the onset depending on the initial magnetic field and the plasma circularization radius. Although the jets are highly variable and misaligned with the BH spin axis, they sustain ≳10 ^50 erg s ^−1 , comparable to that inferred for long gamma-ray bursts (LGRBs). We identify magnetic-flux inversions driven by the stochastic dynamo, leading to the formation of striped jets that could be imprinted in LGRB light curves. These results demonstrate that accretion-disk dynamos provide a robust pathway for jet production in collapsars across a broad range of progenitors.

SN 2024abup is a nearby broad-lined Type Ic supernova (SN Ic-bl) in NGC 0681 at a distance of 23.3 ±1.6 Mpc. As energetic explosions of massive stars, SNe Ic-bl are considered a plausible site for rapid-neutron capture nucleosynthesis (r-process) and chemical enrichment from short-lived progenitors. They may also contribute to dust production in the early Universe. We present JWST near- to mid-infrared (NIR+MIR) observations (1-14 micron) of SN Ic-bl 2024abup at +41 days after the V band maximum (+54 days after explosion), the first-ever JWST+MIR observation of a SN Ic-bl along with radio and optical data. Using the spectral synthesis code SUMO, we identify the observed broad IR line features in SN 2024abup and find significant contributions from C, O, Mg, and carbon monoxide (CO) – the earliest detection of molecules in a core-collapse SN so far. The spectrum shows continuum emission at wavelengths greater than 1.5 micron, which could be explained by dust – preexisting, newly formed, or a combination-heated by the SN. We do not find compelling evidence for infrared signatures of r-process elements, though our search is hampered by the presence of many broad and blended features from the non-r-process elements. These new observations indicate that SNe Ic-bl could be a contributor to early-universe dust production, and suggest that if r-process elements are produced, revealing their presence from spectra requires very high-quality data and models to disentangle blends.

The extreme Oe star AzV 493 is known to show unusual photometric and spectroscopic variability that suggest the presence of an unseen companion in a highly eccentric and long-period (7.3 or 14.6 yr) orbit. We obtained a Chandra/ACIS observation near the putative periastron for the 7.3 yr orbit to test for transient X-ray emission that would confirm its binary nature. Our data only place an upper limit to the X-ray luminosity of L-X < 2.5 & times; 10(33) erg s(-1) based on the 0.5-8 keV flux limit. Additionally, we obtained 4 new spectroscopic observations with the M2FS spectrograph at Magellan and 20 archive FLAMES/GIRAFFE and X-Shooter spectra from ESO/VLT to further constrain the possibility of radial velocity (RV) variation. Statistical analysis of the RV measurements yields inconclusive results regarding the existence of variations. We discuss possible mass limits for a potential companion, which may be a black hole, in the event that the variations are real. The violet-to-red Balmer ratio has also recently inverted, which may be a further indication of a companion.

We measure transverse proper-motion velocities of LMC Wolf-Rayet stars using Gaia DR3 astrometry. The combined velocity distribution of WNh, O If*/WN, and late-type WN (WNL) very massive stars (>100 M-circle dot; VMS) shows both slow, unejected objects (v(perpendicular to) < 10 km s(-1)) and stars dominated by fast runaway velocities (v(perpendicular to) > 24 km s(-1)). This supports expectations that VMS ages are comparable to the dynamical ejection timescale (similar to 1.5 Myr). These kinematics share similarities with those of lower-luminosity, classical WNh, O If*/WN, and WNL stars, as well as the SMC field OB stars, suggesting that dynamical ejections may also dominate these populations. In contrast, both single and binary early-type WN (WNE) stars are ejected populations that show single-peaked velocity distributions, suggesting a different ejection mechanism(s). We speculate that single WNE stars might result from explosive mergers onto the shell-burning layer, thereby stripping the H envelope. Binary WC stars appear to be faster (median v(perpendicular to) = 54 km s(-1)) and have higher luminosities than singles (median v(perpendicular to) = 38 km s(-1)), suggesting that single WC stars are not descendants of the binaries. Thus, the binaries are probably stripped by mass transfer, while the WC singles likely originate from another process. The high velocities of binary WC stars are consistent with some predictions that lower-mass clusters generate fast dynamical ejections. Single WC and WN3/O3 stars have ambiguous kinematics, but both show high v(perpendicular to) (median similar to 38 km s(-1)), possibly linked to their lower masses.

The catalog of Gravitational Wave (GW) events is rapidly growing, providing key insights into the evolution of massive binaries and compact object formation. However, a key challenge is to explain the origin of exceptional events such as GW190814, among the most asymmetric mass-ratio mergers to date (q≈ 0.1). We show that it shares an evolutionary pathway with the most unequal mass Galactic High Mass X-ray Binary (HMXB) 4U 1700-37/ HD 153919. We demonstrate this unique connection by utilizing a rich set of existing observational constraints for the HMXB and compute detailed binary evolution models to explain its formation history. We find that conservative mass transfer, along with a directed natal kick are essential to explain its current state. We show that this system is unlikely to form a GW source due to a failed Common Envelope (CE) phase in the future, in agreement with previous work. With additional models, we show that a similar pathway naturally forms GW190814-like events, provided the first phase of mass transfer remains conservative, and the first-born (lower mass) compact object receives a large natal kick (≳ 100 km/s) for the subsequent CE phase to be successful and form a asymmetric mass-ratio GW source. Anchored by the number of analogous Galactic HMXBs, we estimate rates for such GW events, which broadly agree with their observed rate. Our work demonstrates a unified formation pathway for highly asymmetric mass-ratio HMXBs and GW events. Moreover, it highlights the critical role of finding and characterizing local analogs in different evolutionary phases, and using them as a bridge to understand the origin of GW sources, especially the outliers like GW190814.

Interacting binary star systems play a critical role in many areas of astrophysics. One interesting example of a binary merger product are Thorne-Żytkow Objects (TŻOs), stars that look like red supergiants but contain neutron stars at their cores. TŻOs were theorized nearly five decades ago, and significant work has gone into understanding the physics of their formation, evolution, and stability. Several searches for TŻO candidates have also been carried out. Whether or not TŻOs could even exist or if they would be stable after formation has also been investigated. Understanding the existence and possible prevalence of TŻOs would have important effects on our understanding of binary evolution, stellar mergers, and inform binary population synthesis models. In this chapter, we review the formation channels, evolution and structure, final fates, and observable signatures of TŻOs, as well as candidates in the literature, from the inception of TŻO theory to recent progress in the field.

With upcoming data from Roman, Gaia DR4, and spectroscopic surveys, we will soon have an unprecedented dataset of Milky Way black holes (BHs) to constrain their formation and evolution. To prepare, we simulate the intrinsic population of Milky Way BHs with cogsworth, self-consistently accounting for their binary evolution and trajectories through the Galactic potential. We report the rate, demographics, and kinematics of these BHs, and their sensitivity to 32 variations in binary evolution, supernova physics, and Galactic potentials. In the fiducial model, ~$1.7 \times 10^8$ BHs have formed in the Milky Way (though this total spans an order of magnitude across our variations), where the vast majority (~91%) are currently isolated and ~3% have escaped the Galaxy. Most of the ~$10^7$ BHs in binaries have another BH or a white dwarf companion, but ~$10^5$ retain a luminous stellar companion. BHs are distributed more diffusely than visible stars, with a scale height around ~$2.5\times$ larger. BH masses correlate with present-day location: the most massive BHs are preferentially close to the Galactic plane. This correlation is especially strong for BH-star binaries, which separate into tight, low-mass post-common-envelope systems and wide, high-mass non-interacting ones. The BH mass distribution and kinematics are highly sensitive to the remnant mass prescription and natal kick model, so observations could constrain explodability criteria and BH kicks. Accounting for the time-evolution of the Galactic potential more than doubles the escape fraction and increases the bound population's scale height by ~20%, whilst neglecting binary interactions overestimates it by 30%.

Understanding the transfer of mass and angular momentum in binary interactions is crucial for modeling the evolution of any interacting binary after the first mass-transfer phase. Mass-transfer physics assumptions shape the predictions for later stages of binary evolution, such as the immediate progenitors of stripped-envelope supernovae and gravitational-wave mergers. We constrain the efficiency and stability of thermal-timescale mass transfer in massive binary evolution using the observed population of 62 massive interacting binaries on the main sequence ("Algols") in the Milky Way and the Large and Small Magellanic Clouds. We find that purely conservative or nonconservative mass transfer cannot explain the current mass ratio and orbital period of all massive Algols. Angular momentum conservation rules out conservative mass transfer in similar to 28% of massive Algols in the SMC. About three-quarters of all massive Algols are consistent with having undergone inefficient mass transfer (less than or similar to 50%), while the remaining systems, mostly residing in the LMC and the Milky Way, require mass transfer to have been more efficient than 25%. For our fiducial assumption on the extent of envelope stripping, the current sample of massive Algols does not require mass transfer to be efficient at the shortest orbital periods (similar to 2 days) at any metallicity. We find evidence that mass transfer on the main sequence needs to be stable for initial accretor-to-donor mass ratios as unequal as similar to 0.6. Unless biased by observational selection effects, the massive Algols in the SMC seem to have undergone less efficient mass transfer than those in the LMC and the Milky Way.

The origin of the large-scale poloidal magnetic field required to power relativistic jets in collapsars remains uncertain. While such a field may be inherited during PNS collapse, the efficiency of this process is unclear, motivating an in situ mechanism to generate poloidal fields out of the predominantly toroidal fields produced by stellar differential rotation. We present the first 3D general-relativistic magnetohydrodynamic collapsar simulations initialized with toroidal magnetic field profiles that closely follows those of pre-collapse stellar models. As the toroidal field in the disk becomes dynamically important, it seeds the dynamo, producing coherent poloidal magnetic loops that appear at ∼𝒪(100) gravitational radii and are then advected inward along paths that may deviate from the disk midplane. The resulting poloidal fields thread the black hole (BH) and launch highly variable, wobbling relativistic jets on timescales of order seconds, with the onset depending on the initial magnetic field and the plasma circularization radius. Although the jets are highly variable and misaligned with the BH spin axis, they sustain ≳ 10^50 erg s^-1, comparable to that inferred for long gamma-ray bursts (LGRB). We identify magnetic-flux inversions driven by the stochastic dynamo, leading to the formation of striped jets that could be imprinted in LGRB light curves. These results demonstrate that accretion-disk dynamos provide a robust pathway for jet production in collapsars across a broad range of progenitors.

Context. The blue supergiant (BSG) domain contains a large variety of stars whose past and future evolutionary paths are still highly uncertain. Since binary interaction plays a crucial role in the fate of massive stars, investigating the multiplicity among BSGs helps shed light on the fate of such objects. Aims. We aim to estimate the binary fraction of a large sample of BSGs in the Small Magellanic Cloud (SMC) within the Binarity at LOw Metallicity (BLOeM) survey. In total, we selected 262 targets with spectral types B0-B3 and luminosity classes I-II. Methods. This work is based on spectroscopic data collected by the FLAMES instrument, mounted on the Very Large Telescope, which gathered nine epochs over three months. Our spectroscopic analysis for each target includes the individual and peak-to-peak radial velocity measurements, an investigation of the line profile variability, and a periodogram analysis to search for possible short- and long-period binaries. Results. By applying a 20 km s(-1) threshold on the peak-to-peak radial velocities above which we would consider the star to be binary, the resulting observed spectroscopic binary fraction for our BSG sample is 23 +/- 3%. An independent analysis of line profile variability reveals 11 (plus 5 candidates) double-lined spectroscopic binaries and 32 (plus 41 candidates) single-lined spectroscopic binaries. Based on these results, we estimated the overall observed binary fraction in this sample to be 34 +/- 3%, which is close to the computed intrinsic binary fraction of 40 +/- 4%. In addition, we derived reliable orbital periods for 41 spectroscopic binaries and potential binary candidates, among which there are 17 eclipsing binaries, including 20 SB1 and SB2 systems with periods of less than 10 days. We reported a significant drop in the binary fraction of BSGs with spectral types later than B2 and effective temperatures less than 18 kK, which could indicate the end of the main sequence phase in this temperature regime. We found no metallicity dependence in the binary fraction of BSGs, compared to existing spectroscopic surveys of the Galaxy and Large Magellanic Cloud.

One of the main challenges in modeling massive stars to the onset of core collapse is the computational bottleneck of nucleosynthesis during advanced burning stages. The number of isotopes formed requires solving a large set of fully coupled stiff ordinary differential equations, making the simulations computationally intensive and prone to numerical instability. To overcome this barrier, we design a nuclear neural network (NNN) framework with multiple hidden layers to emulate nucleosynthesis calculations and conduct a proof of concept to evaluate its performance. The NNN takes the temperature, density, and composition of a burning region as input and predicts the resulting isotopic abundances along with the energy generation and loss rates. We generate training sets for initial conditions corresponding to oxygen core depletion and beyond using large nuclear reaction networks, and compare the predictions of the NNNs to results from a commonly used small net. We find that the NNNs improve the accuracy of the electron fraction by 280%–660%, the average atomic and mass numbers by 150%–360%, and the nuclear energy generation by 250%–750%, consistently outperforming the small network across all time steps. They also achieve significantly better predictions of neutrino losses on relatively short timescales, with improvements ranging from 100% to 1,000,000%. While further work is needed to enhance their accuracy and applicability to different stellar conditions, integrating NNN-trained models into stellar evolution codes is promising for facilitating the large-scale generation of core-collapse supernova progenitors with higher physical fidelity.

GW231123, the most massive binary black hole (BH) merger observed to date, involves component BHs with masses inside the pair-instability mass gap and unusually high spins. This challenges standard formation channels such as classical stellar evolution and hierarchical mergers. However, stellar rotation and magnetic fields, which have not been systematically incorporated in prior models, can strongly influence the BH properties. We present the first self-consistent simulations tracking a massive, low-metallicity helium star from helium core burning through collapse, BH formation, and post-BH formation accretion using 3D general-relativistic magnetohydrodynamic simulations. Starting from a 250 M _⊙ helium core, we show that collapse above the pair-instability mass gap, aided by rotation and magnetic fields, drives mass loss through disk winds and jet launching. This enables the formation of highly spinning BHs within the mass gap and reveals a BH spin–mass correlation. Strong magnetic fields extract angular momentum from the BH through magnetically driven outflows, which in turn suppress accretion, resulting in slowly spinning BHs within the mass gap. In contrast, stars with weak fields permit nearly complete collapse and spin-up of the BH to a ≈ 1. We show that massive low-metallicity stars with moderate magnetic fields naturally produce BHs whose masses and spins match those inferred for GW231123, and are also consistent with those of GW190521. The outflows launched during collapse may impart a BH kick, which can induce spin–orbit misalignment and widen the postcollapse orbit, delaying the merger. These outflows could further drive short-lived, high-luminosity jets comparable to the most energetic γ -ray bursts, offering a potential observational signature of such events in the early Universe.

We present nucleosynthesis and light-curve predictions for a new site of the rapid neutron capture process ( r -process) from magnetar giant flares (GFs). Motivated by observations indicating baryon ejecta from GFs, J. Cehula et al. proposed that mass ejection occurs after a shock is driven into the magnetar crust during the GF. We confirm using nuclear reaction network calculations that these ejecta synthesize moderate yields of third-peak r -process nuclei and more substantial yields of lighter r -nuclei, while leaving a sizable abundance of free neutrons in the outermost fastest expanding ejecta layers. The final r -process mass fraction and distribution are sensitive to the relative efficiencies of α -capture and n -capture freeze-outs. We use our nucleosynthesis output in a semianalytic model to predict the light curves of novae breves, the transients following GFs powered by radioactive decay. For a baryonic ejecta mass similar to that inferred of the 2004 Galactic GF from SGR 1806-20, we predict a peak UV/optical luminosity of ∼10 ^39 –10 ^40 erg s ^−1 at ∼10–15 minutes, rendering such events potentially detectable to several Mpc following a gamma-ray trigger by wide-field transient monitors such as ULTRASAT/UVEX. The peak luminosity and timescale of the transient increase with the GF strength due to the larger ejecta mass. Although GFs likely contribute 1%–10% of the total Galactic r -process budget, their short delay-times relative to star formation make them an attractive source to enrich the earliest generations of stars.

In close binary star systems, common envelope evolution (CEE) may occur after a previous phase of mass transfer. Some isolated formation channels for double neutron star binaries suggest that the donor of CEE was the accretor of a previous phase of stable mass transfer. Accretion should substantially alter the structure of the donor, particularly by steepening the density gradient at the core-envelope interface and rejuvenating the star. We study the CEE of a donor that was the accretor of a previous phase of stable mass transfer and has a rejuvenated structure. We perform 3D hydrodynamics simulations of the CEE of an 18 M _⊙ supergiant with a 1.4 M _⊙ companion using rejuvenated and non-rejuvenated 1D stellar models for the donor. We compare the two simulations to characterize the effect of the rejuvenation on the outcome of the common envelope phase and the shape of the ejecta. We find that accounting for a previous phase of mass transfer reduces the duration of the inspiral phase by a factor of two, likely due to the different structures in the outer layers of the donor. In the rejuvenated case, the simulations show more equatorially concentrated and asymmetric ejecta, though both cases display evidence for the formation of a pressure-supported thick circumbinary disk. During the dynamical inspiral phase, the impact of rejuvenation on the unbinding of the envelope is unclear; we find that rejuvenation decreases the amount of unbound mass by 20%–40% depending on the energy criterion used.

Common phases of mass transfer in massive stellar binaries are case A (during the donor’s main sequence) and case B (after the donor’s main sequence but before helium core depletion). Most stars see their radii significantly grow after the main sequence, making case B more common. However, very massive stars may already undergo significant expansion during the main sequence increasing the probability of case A mass transfer. We find that using convective boundary mixing informed by the width of the main sequence in 30 Doradus, case A mass transfer dominates for donor masses ≳75 M ⊙ . This is not the case without convective boundary mixing or in the stellar models commonly used in rapid binary population synthesis. Therefore, case A mass transfer may be more dominant than commonly assumed, with potential impact on rates of all post mass transfer binaries, from Wolf–Rayet-O-type binaries, to X-ray binaries and gravitational wave progenitors.

Rapidly rotating classical OBe stars have been proposed as the products of binary interactions, and the fraction of Be stars with compact companions implies that at least some are. However, to constrain the interaction physics spinning up the OBe stars, a large sample of homogeneously analysed OBe stars with well-determined binary characteristics and orbital parameters are required. We investigate the multiplicity properties of a sample of 18 Oe, 62 Be, and two Of?p stars observed within the BLOeM survey in the Small Magellanic Cloud. We analyse the first nine epochs of spectroscopic observations obtained over approximately three months in 2023. Radial velocities (RVs) of all stars are measured. Applying commonly-used binarity criteria we classify objects as binaries, binary candidates, and apparently single (RV stable) objects. We further inspect the spectra for double-lined spectroscopic binaries and cross-match with catalogues of X-ray sources and photometric binaries. We classify 14 OBe stars as binaries, and an additional 11 as binary candidates. The two Of?p stars are apparently single. Two more objects are most likely currently interacting binaries. Without those, the observed binary fraction for the OBe sample (78 stars) is f_OBe_obs=0.18+/-0.04 (f_obs_cand=0.32+/-0.05 including candidates). This fraction is less than half of that measured for OB stars in BLOeM. Combined with the lower fraction of SB2s, this suggests that OBe stars have indeed fundamentally different binary properties than OB stars. We find no evidence for OBe binaries with massive compact companions, in contrast to expectations from binary population synthesis. Our results support the binary scenario as an important formation channel for OBe stars, as post-interaction binaries may have been disrupted or the stripped companions of OBe stars are harder to detect.

Context. Common envelope evolution of a massive star and a neutron star companion has two possible outcomes: the formation of a short-period binary (a potential gravitational wave source progenitor) or the merger of the massive star with the neutron star. If the binary merges, a structure may form, comprised of a neutron star core surrounded by a large diffuse envelope, known as a ThorneZ(center dot)ytkow object (TZO). The predicted appearance of this hypothetical class of star is very similar to that of a red supergiant, making it difficult to identify them in observations. Aims. Our objective is to understand the properties of systems that are potential TZO progenitors; specifically, binary systems that enter a common envelope phase with a neutron star companion. We also aim to distinguish those that have been through a previous stable mass transfer phase, which can rejuvenate the accretor. Methods. We used the rapid population synthesis code COMPAS at solar metallicity, with the common envelope efficiency parameter set to unity, to determine the population demographics of TZOs. We used one-dimensional (1D) evolutionary TZO models from the literature to determine a fit for the TZO lifetime to estimate the current number of TZOs in the Milky Way, as well as to assess core disruption during the merger. Results. We explored the progenitors in the Hertzsprung-Russell diagram, calculated the formation rates, and investigated the kinematics of the progenitor stars. We find that the vast majority (approximate to 92%) of TZ(center dot)O progenitors in our population have experienced mass transfer and are rejuvenated prior to their formation event. For the Milky Way, we estimate a TZO formation rate of approximate to 4 x 10(4 )yr(-1), which results in approximate to 5 +/- 1 TZOs at present.

Kilonovae, the ultraviolet/optical/infrared counterparts to binary neutron star mergers, are an exceptionally rare class of transients. Optical follow-up campaigns are plagued by contaminating transients, which may mimic kilonovae but do not receive sufficient observations to measure the full photometric evolution. In this work, we present an analysis of the multiwavelength dataset of supernova (SN) 2025ulz, a proposed kilonova candidate following the low-significance detection of gravitational waves originating from the potential binary neutron star merger S250818k. Despite an early rapid decline in brightness, our multiwavelength observations of SN 2025ulz reveal that it is a type IIb SN. As part of this analysis, we demonstrate the capabilities of a novel quantitative scoring algorithm to determine the likelihood that a transient candidate is a kilonova, based primarily on its three-dimensional location and light-curve evolution. We also apply our scoring algorithm to other transient candidates in the localization volume of S250818k and find that, at all times after the discovery of SN 2025ulz, there are ≥4 candidates with a score comparable to SN 2025ulz, indicating that the kilonova search may have benefited from the additional follow-up of other candidates. During future kilonova searches, this type of scoring algorithm will be useful to rule out contaminating transients in real time, optimizing the use of valuable telescope resources.

Core-collapse supernova feedback models in hydrodynamical simulations typically assume that all stars evolve as single stars. However, the majority of massive stars are formed in binaries and multiple systems, where interactions with a companion can affect stars’ subsequent evolution and kinematics. We assess the impact of binary interactions on the timing and spatial distribution of core-collapse supernovae, using cogsworth simulations to evolve binary star populations, and their subsequent galactic orbits, within state-of-the-art hydrodynamical zoom-in galaxy simulations. We show that binary interactions: (a) displace supernovae, with ∼13% of all supernovae occurring more than 0.1 kpc from their parent cluster; and (b) produce delayed supernovae, such that ∼25% of all supernovae occur after the final supernova from a single star population. Delays are largest for low-mass merger products, which can explode more than 200 Myr after a star formation event. We characterize our results as a function of: (1) initial binary population distributions, (2) binary physics parameters and evolutionary pathways, (3) birth cluster dissolution assumptions, and (4) galaxy models (which vary metallicity, star formation history, gravitational potential, and simulation codes), and show that the overall timing and spatial distributions of supernovae are surprisingly insensitive to most of these variations. We provide metallicity-dependent analytic fits that can be substituted for single-star subgrid feedback prescriptions in hydrodynamical simulations, and discuss some of the possible implications for binary-driven feedback in galaxies, which may become particularly important at high redshift.