Explosions Driven by the Coalescence of a Compact Object with the Core of a Massive-Star Companion Inside a Common Envelope: Circumstellar Properties, Light Curves, and Population Statistics

TL;DR

This paper models explosions caused by compact objects merging with the cores of massive stars during common envelope phases, revealing their circumstellar properties, light curves, and population statistics, especially highlighting black hole mergers.

Contribution

It introduces hydrodynamic models of binary coalescence to predict circumstellar structures and light curves of merger-driven explosions, emphasizing the role of black holes and CSM interaction.

Findings

Merger-driven explosions are brightened by CSM interaction by up to three magnitudes.

Black hole mergers produce the brightest explosions with $M_V ext{~} -18$ to $-19$.

Black hole coalescence events constitute about 50% of merger-driven explosions.

Abstract

We model explosions driven by the coalescence of a black hole or neutron star with the core of its massive-star companion. Upon entering a common envelope phase, a compact object may spiral all the way to the core. The concurrent release of energy is likely to be deposited into the surrounding common envelope, powering a merger-driven explosion. We use hydrodynamic models of binary coalescence to model the common envelope density distribution at the time of coalescence. We find toroidal profiles of material, concentrated in the binary's equatorial plane and extending to many times the massive star's original radius. We use the spherically-averaged properties of this circumstellar material (CSM) to estimate the emergent light curves that result from the interaction between the blast wave and the CSM. We find that typical merger-driven explosions are brightened by up to three magnitudes by CSM interaction. From population synthesis models we discover that the brightest merger-driven explosions, $M_V \sim -18$ to $-19$, are those involving black holes because they have the most massive and extended CSM. Black hole coalescence events are also common; they represent about 50% of all merger-driven explosions and approximately 0.3% of the core-collapse rate. Merger-driven explosions offer a window into the highly-uncertain physics of common envelope interactions in binary systems by probing the properties of systems that merge rather than eject their envelopes.

Equations

$$M_{\mathrm{BH},\mathrm{final}}=M_2+M_{\mathrm{acc}}=M_2+f_{\mathrm{acc}}{M}_{\mathrm{core}}.$$

$$\Delta E_{\mathrm{mec}}\sim 3\times 10^{50}\left(\frac{f_{\mathrm{mec}}}{10^{-3}}\right)\left(\frac{\eta }{0.1}\right)\left(\frac{M_{\mathrm{acc}}}{2M_{\odot }}\right)\mathrm{erg}.$$

$$\Delta E_{\mathrm{BZ}}\sim 10^{52}{a}^2{\left(\frac{M_{\mathrm{BH}}}{3M_{\odot }}\right)}^2{\left(\frac{B}{10^{15}\mathrm{G}}\right)}^2\left(\frac{t_{\mathrm{acc}}}{10^2\mathrm{s}}\right)\mathrm{erg},$$

$$\Delta E\approx \frac{d{M}_{\mathrm{CSM}}}{M_{\mathrm{tot}}}E,$$

$$\frac{d{E}_{\mathrm{rad}}}{d{M}_{\mathrm{CSM}}}\approx \frac{E}{M_{\mathrm{tot}}}\frac{r}{R_{\mathrm{ph}}},$$

$$\frac{E_{\mathrm{rad}}}{E}$$

$$\approx \frac{4\pi }{R_{\mathrm{ph}}{M}_{\mathrm{tot}}}{\int }_{R_{\ast }}^{R_{\mathrm{ph}}}\rho r^{3}dr.$$

$$\rho (r)={\rho }_{\mathrm{ph}}{\left(\frac{r}{R_{\mathrm{ph}}}\right)}^{-n},$$

$$M_{\mathrm{CSM}}\approx 4\pi {\rho }_{\mathrm{ph}}{R}_{\mathrm{ph}}^3\times \{\begin{array}{ll}1& n=2,\\ \ln (R_{\mathrm{ph}}/R_{\ast })& n=3,\\ R_{\mathrm{ph}}/R_{\ast }& n=4,\end{array}$$

$$\frac{E_{\mathrm{rad}}}{E}$$

$$\approx \frac{4\pi {\rho }_{\mathrm{ph}}{R}_{\mathrm{ph}}^3}{2M_{\mathrm{tot}}},$$

$$\approx \frac{1}{2}\frac{M_{\mathrm{CSM}}}{M_{\mathrm{tot}}}(\mathrm{for}n=2),$$

$$\frac{E_{\mathrm{rad}}}{E}$$

$$\approx \frac{4\pi {\rho }_{\mathrm{ph}}{R}_{\mathrm{ph}}^3}{M_{\mathrm{tot}}},$$

$$\approx \frac{1}{\ln \left(\frac{R_{\mathrm{ph}}}{R_{\ast }}\right)}\frac{M_{\mathrm{CSM}}}{M_{\mathrm{tot}}}(\mathrm{for}n=3),$$

$$\frac{E_{\mathrm{rad}}}{E}$$

$$\approx \frac{4\pi {\rho }_{\mathrm{ph}}{R}_{\mathrm{ph}}^3}{M_{\mathrm{tot}}}\ln \left(\frac{R_{\mathrm{ph}}}{R_{\ast }}\right),$$

$$\approx \left(\frac{R_{\ast }}{R_{\mathrm{ph}}}\right)\ln \left(\frac{R_{\mathrm{ph}}}{R_{\ast }}\right)\frac{M_{\mathrm{CSM}}}{M_{\mathrm{tot}}}(\mathrm{for}n=4).$$

$$M_{V,\mathrm{peak}}\approx -18.9$$

$$-0.229{\log }_{10}(M_{\ast }/30M_{\odot })$$

$$-1.41{\log }_{10}(E/10^{51}\mathrm{erg})$$

$$-2.20{\log }_{10}(q/0.1)$$

$$M_V\approx -11.42$$

$$+1.25{\log }_{10}(M_{\ast }/M_{\odot })$$

$$-2.08{\log }_{10}(E/10^{50}\mathrm{erg}).$$

Full text

binaries: close, methods: numerical

††software: Athena++, Stone et al. (in preparation) http://princetonuniversity.github.io/athena, Astropy (Astropy Collaboration et al., 2013), SNEC (Morozova et al., 2015a), COMPAS (Stevenson et al., 2017; Barrett et al., 2018; Vigna-Gómez et al., 2018)

1 Introduction

Binary and multiple systems are ubiquitous among massive stars. Of these systems, a large fraction are at separations so close that they will interact over the stars’ lifetimes (Sana et al., 2012; de Mink et al., 2014; Moe & Di Stefano, 2017). As these multiple-star systems evolve to leave behind compact object stellar remnants, the stage is set for interactions between the evolving stars and the stellar remnants. In some cases, a phase of escalating, unstable mass transfer from a massive star donor onto a compact-object companion leads to a common envelope phase (Paczynski, 1976), in which the compact object is immersed within the envelope of the massive star and spirals closer to the massive star’s core (e.g. Taam et al., 1978; Iben & Livio, 1993; Terman et al., 1995; Taam & Sandquist, 2000; Armitage & Livio, 2000; Ivanova et al., 2013; De Marco & Izzard, 2017).

Common envelope interactions can lead to either the ejection of the shared, gaseous envelope and a surviving, binary pair or to the merger of the donor star core with the companion. The distinction between these cases for a given binary remains highly uncertain, but is of significant interest given its importance (e.g. Belczynski et al., 2002; Kalogera et al., 2004, 2007; Belczynski et al., 2008; Andrews et al., 2015; Tauris et al., 2017; Vigna-Gómez et al., 2018), for example, to estimating rates of compact object mergers and associated gravitational-wave transients (Abbott et al., 2016a, b; Abbott, 2017a, b, c; The LIGO Scientific Collaboration & the Virgo Collaboration, 2018a, b). In the case of massive stars interacting with lower-mass, compact object companions, the theoretical expectation is that only donor stars with the most-extended, weakly bound hydrogen envelopes are susceptible to ejection, while the remainder of systems are likely to merge (Kruckow et al., 2016).

What happens when a compact object merges with the helium core of a massive star? At least two possible outcomes have been suggested. Perhaps a neutron star embedded in a stellar envelope could burn stably, forming a Thorne-Zytkow object (Thorne & Zytkow, 1977; Podsiadlowski et al., 1995; Podsiadlowski, 2007; Levesque et al., 2014; Moriya, 2018). Alternatively, it is possible for the angular momentum of the merger to lead to the formation of a disk around the compact object (be it a neutron star or a black hole), from which material forms a rapidly accreting neutrino-cooled disk (Houck & Chevalier, 1991; Chevalier, 1993, 1996; Fryer et al., 1996; Fryer & Woosley, 1998; Popham et al., 1999; MacFadyen & Woosley, 1999; Fryer et al., 1999; MacFadyen et al., 2001; Zhang & Fryer, 2001; Lee & Ramirez-Ruiz, 2006; Chen & Beloborodov, 2007; Barkov & Komissarov, 2011; Chevalier, 2012; Fryer et al., 2014; Song & Liu, 2019). Accretion of the surrounding core material liberates on the order of $\eta M_{\odot}c^{2}\sim\eta 10^{54}$ erg, where $\eta$ is an efficiency factor of order 0.1 (Frank et al., 2002). Much of this energy emerges in neutrinos, which stream freely away from the accretion object (Chevalier, 1993; Fryer et al., 1996; Zhang & Fryer, 2001). A fraction, however, emerges as Poynting flux or mechanical energy (either a disk wind or collimated outflow) and can feed back on the surroundings, powering a blast wave with energy similar to a supernova (Fryer & Woosley, 1998; Zhang & Fryer, 2001; Chevalier, 2012; Dexter & Kasen, 2013; Fryer et al., 2014; Song & Liu, 2019; Gilkis et al., 2019; Soker, 2019; Soker et al., 2019).

Chevalier (2012) observed that, in the case that the merger was initiated by a preceding common envelope phase, the blast wave necessarily interacts with the dense surrounding medium of the common envelope ejecta (Soker & Gilkis, 2018; Soker, 2019; Gilkis et al., 2019; Soker et al., 2019). The resulting transient is described as a ”common envelope jets supernova” by Soker et al. (2019). One of the key points that Chevalier (2012) and Soker et al. (2019) mention is that the distribution of common envelope ejecta is crucial in shaping the observed light curve (see Kleiser et al., 2018, for a similar discussion in the context of rapidly-fading supernovae). Here, we pursue this line of examination.

Light curves from supernova explosions of single stars within a dense circumstellar material (CSM) have been considered as the origin of the variations in type II supernova for many years (Chugai, 1997, 2001; Chugai & Danziger, 1994; Chugai et al., 2004; Smith & McCray, 2007; Chevalier & Irwin, 2011; Moriya et al., 2013a; Ofek et al., 2013b; Dessart et al., 2015; Morozova et al., 2017), and even a population of objects that appear to transition from type I to type II (Margutti et al., 2017). Specifically type IIn, with their narrow line features, are known to be interacting with a slow moving CSM (Kiewe et al., 2012; Taddia et al., 2013). But the origin of the material around the pre-SN star is not yet clear. Ideas include late stage stellar winds or small outbursts of gas (Quataert & Shiode, 2012; Shiode & Quataert, 2014; Fuller, 2017) or formations of a disk-like torus (Andrews & Smith, 2018; McDowell et al., 2018). With this paper we add to the calculations by Chevalier (2012) the light curves expected from an engine-driven explosion inside a merger ejecta profile. We show how the atypical CSM distribution leads to light curves powered in part by CSM interaction, resembling type IIn and IIL and lacking a plateau phase (e.g. Das & Ray, 2017; Morozova & Stone, 2018; Eldridge et al., 2018).

We briefly review the engine-driven explosion model in 2. We then self-consistently model the circum-merger density distribution from the common envelope phase using a three-dimensional hydrodynamic simulation in Section 3. To explore the impact of this CSM on the resultant engine-driven explosions, we produce a number of spherically-symmetric radiative transfer models of blast waves interacting with the (spherically-averaged) common-envelope ejecta profiles in Section 4. Next, we map the expected populations of compact object–core mergers and their resultant transients in Section 5. Finally, in Section 6 we study the imprint of CSM on merger-driven explosions, compare to known supernovae, and discuss possible identification strategies. In Section 7 we conclude.

2 Merger-driven explosions

2.1 Inspiral, Merger, and Central Engine

Following its inspiral through the common envelope, a compact object can tidally disrupt and merge with the helium core of the massive star (Zhang & Fryer, 2001; Chevalier, 2012; Soker & Gilkis, 2018; Soker et al., 2019). The disrupted material forms an accretion disk surrounding the now-central compact object. The local densities of more than $10^{2}$ g cm*-3* imply that accretion occuring on a dynamical timescale is very rapid of the order of $10^{-3}$ to $10^{-1}M_{\odot}$ s*-1* (Fryer & Woosley, 1998; Zhang & Fryer, 2001), making these conditions very similar to those of a classical collapsar scenario of rapid accretion onto a newly formed black hole (MacFadyen & Woosley, 1999; MacFadyen et al., 2001; Siegel et al., 2018). Neutrinos mediate the bulk of the accretion luminosity at these accretion rates (MacFadyen & Woosley, 1999), so accretion can occur onto either neutron stars or black holes under these conditions (Fryer et al., 1996). At higher-still accretion rates, neutrino energy deposition can overturn the accretion flow (e.g. Pejcha & Thompson, 2012). This rapid accretion, over approximately the dynamical time of the core (Zhang & Fryer, 2001), $10^{3}$ s, will lead a neutron star companion object to quickly collapse to a black hole.

These conditions of hypercritical accretion from the core onto an embedded black hole set the stage for an explosion powered by this central, accreting engine. In the context of collapsing helium stars stripped of their hydrogen envelopes, the result is a long gamma ray burst (GRB), in which a relativistic jet, perhaps powered by the coupling between the magnetic field in the accretion disk coupling to the rotational energy of the black hole (via the Blandford & Znajek (1977) process (e.g. Barkov & Komissarov, 2008)) tunnels out of the helium star and is accelerated to high Lorentz factor (MacFadyen & Woosley, 1999; MacFadyen et al., 2001).

In the context of a helium core surrounded by an extended hydrogen envelope, the beamed power of the jet is not expected to be able to tunnel out of the envelope under most circumstances. This is because the jet head must displace the ambient stellar material and, therefore, expands at a rate which balances the ram pressure of this interaction with the momentum flux of the jet (e.g. Matzner, 2003). Following the estimates of Quataert & Kasen (2012)’s Section 2.2, a typical jet whose working surface expands at a few percent the speed of light (as it displaces the surrounding stellar envelope) reaches only approximately $0.03c\times 10^{3}{\rm s}\sim 10^{12}$ cm before the core-accretion event ends and the jet shuts off. Lacking the pressure to continue driving its expansion, the jet is choked by the surrounding gas distribution, and expands laterally, distributing its power more isotropically into the envelope material and powering an outburst (e.g. Murguia-Berthier et al., 2014; Senno et al., 2016). We note that Quataert & Kasen (2012); Woosley & Heger (2012); Dexter & Kasen (2013) discuss an alternate scenario in which the jet might escape if long-lived accretion from the envelope material persists over a duration of hundreds of days. The typical energies range from $10^{50}-10^{52}$ erg (Fryer & Woosley, 1998; Zhang & Fryer, 2001), as we discuss in Section 2.2.

An alternative process of stable burning has been suggested when the engulfed object is a neutron star. In these cases a Thorne-Zytkow Object is said to form – a star with a neutron-star core, which might stably burn hydrogen for up to $10^{5}$ yr (Thorne & Zytkow, 1977; Podsiadlowski et al., 1995; Podsiadlowski, 2007), potentially showing unique surface features (Levesque et al., 2014). If the core of the helium star is near its original density when neutron star enters it, the flow convergence rate is high enough that it seems very difficult to avoid a neutrino-cooled accretion state (Chevalier, 1993; Fryer & Woosley, 1998). However, it is interesting to note that because the neutron star enters the core from the outside in (rather than inside-out as in a collapsar) there is a possibility that feedback from lower accretion-rate common envelope inspiral (e.g. Fryer et al., 1996) would prevent the object from ever reaching the hypercritical, neutrino-cooled accreting branch (see, for example, the trajectories of infall and accretion in MacLeod & Ramirez-Ruiz, 2015; Holgado et al., 2018). Thus, more work is needed to conclusively distinguish between these alternatives.

2.2 Model Adopted: Engine Mass and Energetics

Before merger the donor star consists of a core of mass $M_{\rm core}$ and an envelope of mass $M_{\rm env}$. The total mass of the star is $M_{\ast}=M_{\rm core}+M_{\rm env}$. The compact-object companion that ends up in the center of the donor’s core has mass $M_{2}=qM_{\ast}$. When the two stars merge, $M_{2}$ accretes a fraction, $f_{\rm acc}$, of $M_{\rm core}$. This either causes it to grow if it is already a black hole, or, if it is a neutron star, collapse to a black hole, then grow to a final mass,

[TABLE]

The released energy from the accretion is, therefore, $\Delta E_{\rm acc}=\eta M_{\rm acc}c^{2}$, where $\eta\sim 0.1$. Much of this energy is radiated by neutrino emission. However, some energy emerges mechanically, from a magnetohydrodynamic disk wind (e.g. Feng et al., 2018). The mechanical power is an uncertain fraction of the accretion energy, such that $\Delta E_{\rm mec}=f_{\rm mec}\Delta E_{\rm acc}$. Dexter & Kasen (2013) estimate $f_{\rm mec}\sim 10^{-3}$ given an inflow-outflow model in which the accreting material decreases as a power law with radius, implying

[TABLE]

The energy that emerges in a jet via the Blandford & Znajek (1977) process can be estimated as (Zhang & Fryer, 2001),

[TABLE]

where $B$ is the magnetic field in the disk, $t_{\rm acc}$ is the timescale of accretion, and $a$ is the black hole’s dimensionless spin. Because $M_{\rm acc}\gtrsim M_{\rm BH}$, values of order unity for $a$ are expected regardless of the initial spin state. In what follows, we will assume that the compact object accretes the entire core mass ($f_{\rm acc}=1.0$) and we explore explosion energies between $3\times 10^{50}$ erg and $10^{52}$ erg. We further assume that regardless of the precise mechanism or energetics, the energy is shared roughly spherically with the hydrogen envelope (Chevalier, 2012; Soker & Gilkis, 2018; Soker, 2019).

3 Circumstellar Material Expelled During Binary Coalescence

As a basis for our analysis, we model binary coalescence and the circumbinary ejecta that results from the merger of two example mass ratio binaries. Here we describe our numerical method, the unstable mass exchange that leads the binary to merge, and the CSM mass distribution that this runaway mass transfer creates. The distribution of CSM at large scales, $r\sim 10^{15}$ cm, is of particular importance for the outburst light curves. We therefore focus our numerical modeling on the early phases of runaway, unstable mass exchange that expels this largest-scale CSM.

3.1 Hydrodynamic Models of Binary Coalescence

Our models are simulated within the Eulerian hydrodynamic code Athena++ (Stone, J.M., in preparation), and are based on the methodology described in MacLeod et al. (2018a, b, c). We use a spherical polar coordinate system surrounding the donor star in a binary pair. We model the interaction of this donor star with a softened point mass representing an unresolved companion object. The domain extends over the full three-dimensional solid angle from 0.1 to 100 times the donor-star’s original radius.

The donor star is modeled by polytropic envelope, with structural index $\Gamma=1.35$. The donor has a core mass of 25% its total mass. The gas in the simulation domain follows an ideal-gas equation of state, with index $\gamma=1.35$. These choices are intended to approximately represent a convective, isentropic envelope of a massive star in which radiation pressure is important in the equation of state (in which case $\Gamma=\gamma\rightarrow 4/3$, e.g. MacLeod et al., 2017a; Murguia-Berthier et al., 2017).

We initialize the calculation at a separation slightly smaller than the analytic Roche limit separation, where the donor star overflows its Roche lobe (Eggleton, 1983), and halt the calculation when the companion star has plunged to 10% of the donor’s original radius – the inner boundary of our computational domain. The binary is initialized in a circular orbit and the donor star is initially rotating as a solid body with rotational frequency matching the orbital frequency at the Roche limit separation.

The calculations themselves are carried out in dimensionless units where the donor’s mass, radius, and gravitational constant are all set to unity. They may, therefore, be rescaled to a physical binary of any mass or size. Below, we report on two models which have secondary to donor star mass ratios of $q\equiv M_{2}/M_{\ast}=0.1$ and $q=0.3$. As we will see in Section 5, $q\sim 0.1$ is a very common mass ratio, while $q=0.3$ is near the upper end of the range of events that result in mergers.

3.2 Unstable Mass Transfer Leading to Binary Merger

Mass transfer is unstable in our model binary system in that it runs away to ever increasing rates and drives the binary toward merger. This process begins with Roche lobe overflow of the donor star into the vicinity of its companion. In general, mass transfer proceeds unstably when the loss of material from the donor causes the donor star to increasingly overflow its Roche lobe – either because it grows in radius, or because its Roche lobe shrinks. In binary systems such conditions are often realized in binary pairs where a more massive donor star transfers mass onto a less massive accretor, causing the binary separation to shrink.

In Figure 1, we show the binary system separation as a function of time in our model system, and snapshots of the gas density distribution in the orbital plane. In these figures, we have rescaled our dimensionless simulations to a fiducial donor star mass of $30M_{\odot}$ and radius of $1000R_{\odot}$. In the upper panel, time is zeroed at the time at which the companion plunges within the original donor-star’s radius, $t_{1}$, where $a(t_{1})=R_{\ast}$. Over the preceding 50 years, the binary separation continuously shrinks, at first gradually, but with increasing rapidity (MacLeod et al., 2018a). After the companion object plunges within the donor’s envelope it spirals to the inner boundary of our computational domain (at 10% the donor’s radius) within about 5 orbital cycles, or two years, approximately the orbital period at the donor star’s surface.

Mass loss from the donor star at the expense of orbital energy drives this rapid decrease in binary separation and the pair’s coalescence. Turning our attention to the lower panels of Figure 1, we note that as the donor star overflows its Roche lobe, material is pulled, primarily from the vicinity of the $L_{1}$ Lagrange point, toward the companion object. As the orbital separation decreases, from $1500R_{\odot}$ to $1000R_{\odot}$ to $500R_{\odot}$, the breadth and intensity of this mass transfer stream increase dramatically. MacLeod et al. (2018a) studied the dynamics of this runaway, unstable Roche lobe overflow in detail, and found that the mass loss rate from the donor increases by orders of magnitude over this period. However, MacLeod et al. (2018a) also show that the analytic model of Paczyński & Sienkiewicz (1972) coupled to a point-mass binary orbit evolution model captures the key features of these stages once the specific angular momentum of the ejecta has been measured (e.g. Huang, 1963).

These high mass exchange rates quickly exceed the Eddington limit mass accretion rate that material can accrete onto a compact object companion, and much of the material pulled from the donor is lost to the circumbinary environment (as seen in the snapshots of Figure 1, though in the case of the simulation this is because accretion onto the companion object is not modeled, see MacLeod et al., 2018b). Much of this mass loss occurs near the $L_{2}$ Lagrange point, near the lower-mass companion object (Shu et al., 1979; Pejcha et al., 2016a, b; Metzger & Pejcha, 2017; MacLeod et al., 2018b). By the final panel, where the separation is $100R_{\odot}$, the core of the original donor and companion object are mutually immersed in a significantly extended common envelope that originated from the donor star (Paczynski, 1976). Once immersed, some binary systems deposit enough energy into their environments to expel this envelope. Others do not, and the companion object merges with the core – powering the sort of engine-driven explosions discussed in Section 2.

Because orbital tightening and coalescence of the binary system is a direct result of angular momentum loss to ejected material, the amount of ejecta relates directly to the binary properties. First using semi-analytic scalings (MacLeod et al., 2017b), then hydrodynamic simulation results (MacLeod et al., 2018b), we have found that the expelled mass at the onset of coalescence (defined as mass at $r>R_{\ast}$ at $t=t_{1}$) is always on the order of 25% the mass of the merging companion object (Section 4.2 of MacLeod et al., 2018b). In the calculation shown in Figure 1, which has a mass ratio $q=0.1$, at $t_{1}$ the ejecta mass (measured as the mass at radius greater than the donor’s original radius) is 16% the companion’s mass, or approximately $0.49M_{\odot}$. At the termination of our calculation, when the separation has decreased by a further factor of ten, the ejecta mass has increased to roughly 150% of the companion object’s mass, or $4.45M_{\odot}$. By comparison, our calculation with $q=0.3$ expels nearly identical percentages of mass relative to the more massive black hole: $1.44M_{\odot}$ at a separation equal to the donor’s radius and $13.9M_{\odot}$ in our final snapshot (separation 10% the donor’s radius). In the following, we analyze the distribution of this material in the circumstellar environment.

3.3 Resultant Circumstellar Distribution

Next, we analyze the three-dimensional distribution of debris expelled by the merger episode. To do so, we analyze the final snapshot of our hydrodynamic simulation, when the separation has tightened to one-tenth the donor’s original radius, or $100R_{\odot}$ in our fiducial, $30M_{\odot}$, $1000R_{\odot}$ model. Because the binary separation is tightening extremely rapidly at this phase, material ejected subsequently in the merger does not affect the largest-scale gas distribution prior to the compact object’s merger with the core. This, therefore, is the CSM that any explosive outburst will interact with as it expands, particularly when we consider the crucial scales of interest of $10^{14}$ to $10^{15}$ cm that lie near the photosphere of the explosive transient.

In Figure 2, we show the large-scale density distribution out to 50 times the initial donor radius ($5\times 10^{4}R_{\odot}$, or approximately $3.5\times 10^{15}$ cm). The panels show a slice through the orbital, $x-y$ plane, and an aziumthal average, plotted in $z-R$, perpendicular to the orbital plane. Figure 2 shows that a thick, extended circumbinary torus of expelled material from the donor’s envelope has formed around the merging pair of stars. This torus is roughly azimuthally symmetric, but has distinct structure in polar angle, with relatively evacuated poles and dense equator, representative of the fact that material is flung away from the merging binary in the equatorial plane. Pejcha et al. (2016a, b, 2017) analyzed the thermodynamics of similar outflows and show that heating, arising continuously from internal shocks, and radiative diffusion and cooling very likely regulate the torus scale height. Thus, the precise scale height observed in Figure 2, modeled under the simplification of an ideal-gas equation of state, would likely be modified by the inclusion of more detailed physics.

MacLeod et al. (2018b) analyzed the kinematics of this torus material and found that the radial velocities of the most extended material are low relative to the escape velocity of the original donor star (roughly 100 km s*-1* for our fiducial model). Thus, the majority of these (earlier) ejecta around bound to the merging binary. Some of the material at smaller radii (the later ejecta) is moving more rapidly, at velocities similar to the escape velocity. It therefore collides with the earlier, slow moving ejecta (MacLeod et al., 2018b). Qualitatively, these velocities are similar to other sources of stellar mass loss like winds or non-terminal outbursts, in that they are similar to the giant star’s escape velocity and are much less than the later explosion’s blast wave velocity.

Though the axisymmetric torus structure discussed above is clearly structured in polar angle, for the sake of computational efficiency, we model the interaction of the explosive blast wave with a one-dimensional (spherically-symmetric) density distribution derived from these models. In future work, it may be extremely interesting to relax this simplification. To derive one-dimensional profiles, we spherically-average our model results about the donor star’s core.

These 1D density profiles are shown in Figure 3, in which we compare the unperturbed envelope profile to the cases disturbed by binaries of $q=0.1$ and $q=0.3$. Where the hydrostatic profile has a distinct limb at the donor’s radius, the post-merger profiles show a roughly power-law slope in radius, with approximate scaling of $r^{-3}$, as shown in the lower panel. Comparing the $q=0.1$ and $q=0.3$ results shows that in the higher mass-ratio coalescence, more of the envelope material has been expelled beyond the donor’s original radius, yielding a shallower density fall-off with radius and a profile with more mass at large radii. In both cases, we see that the distribution of ejecta extends to roughly $10^{15}$ cm, with of order a solar mass ($q=0.3$) or a tenth of a solar mass ($q=0.1$) on these scales.

4 Merger-Driven Light Curves

In this section we use analytic and numerical models to understand the properties of merger-driven light curves. We find that the CSM distribution plays a crucial role in shaping these light curves. In Section 4.1 we provide some analytic context for the potential role of CSM. In Section 4.2, we describe our numerical method for 1D radiative hydrodynamics calculations and light curve generation. Finally, in Sections 4.3 and 4.4 we describe the key features and variations across parameter space of these numerical model light curves.

4.1 Analytic Context for Contribution to Radiative Efficiency from CSM Interaction

In order to provide context for the interpretation of our numerical light curve models, in this section we analyze simplified analytic models of emission powered by CSM interaction. In general, CSM interaction enhances the intrinsic luminosity of “cooling emission” from heated, expanding material like supernova ejecta. As ejecta shock-heat by collisions with CSM, their kinetic energy is dissipated and converted into radiation. Depending on the location of the interaction (within or outside the photosphere), this radiation may either escape immediately or adiabatically decay with expansion in the outflow prior to being free to stream out.

In particular, we focus on different radial density profiles of CSM material and the role this plays in shaping transient light curves. In doing so, we summarize and build on a considerable literature that describes how CSM interaction can form a significant contribution or even dominate the radiative power of a transient under certain conditions (e.g. Chugai & Danziger, 1994; Smith & McCray, 2007; Chevalier & Irwin, 2011; Chevalier, 2012; Chevalier & Irwin, 2012; Ginzburg & Balberg, 2012; Moriya & Maeda, 2012; Moriya et al., 2013a, b; Pan et al., 2013; Ofek et al., 2014; Ginzburg & Balberg, 2014; Morozova et al., 2015b, 2017, 2018; Morozova & Stone, 2018; Kleiser et al., 2018; Chandra, 2018).

4.1.1 Thin Shell of CSM

In the simplest version of a CSM interaction, an additional internal energy $\Delta E$ is added to ejecta by sweeping up a thin shell of CSM,

[TABLE]

where $M_{\rm tot}$ is the sum of the explosive ejecta and the swept-up CSM mass internal to the shell, $dM_{\rm CSM}$ is the CSM shell mass, and $E$ is the kinetic energy of the explosive ejecta. If this deposition of internal energy occurs in optically thin regions, all of this energy is radiated, and $\Delta E_{\rm rad}\approx\Delta E$. If the CSM shell lies interior to the photosphere radius, the heated ejecta must continue to expand, with gas internal energy decaying along an adiabat, before they are free to radiate. If we assume that gas specific internal energy decays adiabatically as $r^{-1}$ prior to reaching the photosphere (which is the case when radiation pressure dominates and $P\propto\rho^{4/3}$ along an adiabat), then $\Delta E_{\rm rad}\approx(r/R_{\rm ph})\Delta E$.

4.1.2 Continuous Distributions of CSM

The differential radiated energy due to sweeping up a CSM mass $dM_{\rm CSM}$ at radius $r$ interior to the photosphere radius is

[TABLE]

where $E_{\rm rad}$ is the contribution to the radiated energy due to CSM interaction alone. As before, we are also assuming that gas specific internal energy decays adiabatically as $r^{-1}$. Given a continuous distribution of CSM material distributed between $R_{\ast}$ and $R_{\rm ph}$, we can integrate this expression over radius. In what follows, we will assume that $E$ is constant, which is justified only if $M_{\rm CSM}\ll M_{\rm tot}$ and $E_{\rm rad}\ll E$. Otherwise, the losses in kinetic energy to thermal energy or radiation must be taken into account. We replace $dM_{\rm CSM}=4\pi r^{2}\rho dr$, where $\rho$ is the CSM density, to write

[TABLE]

This integral shows that the dependence of $\rho(r)$ will be critical in determining the CSM contribution to the radiated luminosity.

Let us write a general, power law density form for the CSM that applies from the stellar radius, $R_{\ast}$, to the eventual photosphere radius, $R_{\rm ph}$,

[TABLE]

where $\rho_{\rm ph}$ is the density at the photosphere radius, and we have chosen $R_{\rm ph}$ as a characteristic radius to normalize the power law. We will further adopt the approximation that $R_{\rm ph}\gg R_{\ast}$, under which the total CSM mass can be written,

[TABLE]

for several representative values of $n$.

The most frequently considered form of $\rho$ is that of a steady, spherical wind, $n=2$. Then, equation (6) becomes

[TABLE]

thus retrieving the often quoted result of the increase in radiated energy scaling with the CSM mass as a fraction of the total mass (Chevalier & Irwin, 2011; Ginzburg & Balberg, 2012; Pan et al., 2013).

In our $q=0.3$ model, the one-dimensional profile approximates $n=3$, which yeilds constant mass per logarithmic increase in radius. For $n=3$, equation (6) evaluates to

[TABLE]

where, in the last line, we have used $M_{\rm CSM}$ from equation (8). First, we emphasize that the ratio of $R_{\rm ph}$ to $R_{\ast}$ now affects the radiated luminosity arising from CSM interaction, which is not the case for $n=2$. Therefore, when $R_{\rm ph}\gg R_{\ast}$, the radiated energy from this CSM profile is considerably less than that of the $n=2$ profile. A final way to interpret this result is in terms of the CSM mass at radii similar to the photosphere radius, $r\sim R_{\rm ph}$. From equation (8), this is approximately $M_{\rm CSM}/\ln(R_{\rm ph}/R_{\ast})$. Equation (10) shows that this is the fraction of the CSM mass that contributes significantly to the radiated energy.

In our $q=0.1$ merger model, the scaling of the CSM is steeper, approximately $\rho\propto r^{-4}$. We reevaluate equation (6) for $n=4$ to find,

[TABLE]

Thus, for $R_{\rm ph}\gg R_{\ast}$, the CSM contribution to radiated luminosity is less for $n=4$ than either $n=3$ or $n=2$ given a CSM mass. This result can be interperted in light of equation 8, which shows that the fraction of CSM mass with $r\sim R_{\rm ph}$ is $R_{\ast}/R_{\rm ph}$ for $n=4$.

4.1.3 Interpretation

In the preceding subsection, we have shown that for CSM density profiles that are sufficiently steep, $n\geq 3$, the CSM contribution to the radiated luminosity $E_{\rm rad}/E$ depends on the ratio of the stellar radius over the photosphere radius – the radial extent of the CSM. This important ratio varies in explosions with different size stars of similar mass, or over the time evolution of a given transient as the photosphere radius increases. If $E_{\rm rad}/E$ becomes too small, then CSM interaction does not contribute significantly to the light curve of the transient at a given phase and the bulk of the radiated luminosity comes instead from the adiabatically-expanding blast wave. In this case, the transient assumes more typical supernovae type IIP properties.

These scalings indicate that we expect the CSM to be an important contribution to the $q=0.3$ merger case light curve, because the mass in the CSM is a significant fraction of the total envelope mass, and with $n=3$ there is only logarithmic dependence on $R_{\rm ph}/R_{\ast}$, equation (10). In the $q=0.1$ merger case, in which $n=4$, we expect preferential contribution from CSM interaction at early times in the transient light curves or for particularly extended donors, because either situation maximizes the ratio $R_{\ast}/R_{\rm ph}$, see equation (11).

Finally, we have so far discussed the case in which the CSM distribution extends out to, and perhaps beyond the photosphere radius. This is not necessarily realized. In the case where the CSM terminates at a radius $R_{0}$, for which $R_{0}<R_{\rm ph}$, the radiated luminosity from CSM interaction is computed much as before, but only integrating the mass distribution out to $R_{0}$. This reduces the radiated luminosity due to CSM interaction by a factor similar to $R_{0}/R_{\rm ph}<1$. In the sections that follow, we use this framework to interpret 1D radiative transfer models of explosions interacting with our model CSM distributions.

4.2 1D Radiation Hydrodynamic Models

While the analytic approach highlighted above is useful, it is necessarily simplified. To extend these calculations of the CSM imprint on transient light curves to slightly more realistic scenarios, we need to perform the associated integrations numerically. We utilize the publicly available spherically symmetric (1D) Lagrangian hydrodynamics code SuperNova Explosion Code (SNEC) to calculate bolometric and filtered light curves (Morozova et al., 2015a, b, 2017). The code uses equilibrium-diffusion radiation transport to follow the time dependent radiation hydrodynamics of the expanding blast wave. We choose to set the equation of state using the built-in version of the Paczynski (1983) equation of state, which includes contributions to the total pressure from radiation, ions and electrons based on the composition.

We map our one-dimensional, spherically-averaged profiles, shown in Figure 3, into SNEC. Mass grid cells are customized by the user’s choice of binning in mass, and we have found that the density profiles’ steep decline is best simulated with increasingly fine mass resolution at larger radii. This means that the shock break out is not well resolved (Ensman & Burrows, 1992; Morozova et al., 2015a), but light curve calculations after the first day are robust as shown by Morozova et al. (2015b). We run the simulations for this paper with 456 grid cells.

Even so, the code does not function well for the lowest densities. We therefore are required to restrict the CSM density profile to $\rho>10^{-12}$ g cm*-3* (Morozova & Stone, 2018). The outer radius is therefore $2.6\times 10^{14}$ cm for the fiducial simulation with $q=0.1$ and $1.0\times 10^{15}$ cm for the simulation with $q=0.3$. We note that this restriction is not ideal because it limits the potential interaction-driven luminosity of our models, see equations (10) and (11). In practice, this implies that the CSM maximum radius is often less than the photosphere radius, $R_{0}<R_{\rm ph}$, and the CSM contribution to the eventual radiated luminosity is reduced accordingly, see the discussion of Section 4.1.3. Nonetheless, the truncated profiles do retain more than 95% of the CSM mass in all parameter variations. We assume a roughly solar isotopic composition that matches the hydrogen envelope of a presupernova stellar model of an initially $15M_{\odot}$ star evolved with the MESA code that is included in the SNEC distribution.

To drive the explosion of our models, we adopt a thermal bomb at the inner edge of the envelope domain. This broadly mimics the energy deposition of the quenched jet into the hydrogen envelope, as described in Section 2. Based on Morozova et al. (2015b), we deposit the energy over the innermost $\Delta\mathrm{M_{bomb}}=0.1\mathrm{M_{\odot}}$ over a duration $\Delta\mathrm{t_{bomb}}=0.1$ s (Morozova et al., 2015b, has shown that the model light curves are not particularly sensitive to these parameters, see their Figure 5). The code also offers the option of adding Nickel to the composition. Though Nickel and lanthanide production is possible in merger-driven explosions (e.g. Siegel et al., 2018; Grichener & Soker, 2018), the quantities are uncertain. Because decay of radioactive 56Ni mainly powers the late-time emission, here we choose to focus on Ni-free models of the early light curve dominated by the CSM and the hydrogen envelope (Morozova et al., 2015b). Lacking radioactive material in the ejecta, our models decline rapidly after the ejecta become fully transparent. Emission from the photosphere in SNEC is assumed to follow a thermal blackbody, and thus neglects some line-blanketing effects that may be important for iron-rich ejecta in the U and B bands.

Finally, in Appendix A, we test the sensitivity of our model results to these choices by varying the inner mass (or equivalently, radius) at which energy is deposited, as well as the mass-resolution of the SNEC calculation.

4.3 Imprint of CSM

We begin to explore the imprint of the CSM mass distribution on the explosion light curve in Figure 4, in which we plot luminosities, along with photosphere radii, effective temperature, and gas bulk velocity at the photosphere radius for our fiducial case of a $30M_{\odot}$ and $1000R_{\odot}$ donor star in its initial, hydrostatic equilibrium state (labeled HSE), and following merger with a $3M_{\odot}$ ($q=0.1$) or $9M_{\odot}$ ($q=0.3$) black hole. In each case, the explosion energy is taken to be $10^{51}$ erg, and is injected at 0.1 times the radius of the original donor star, the innermost radius resolved in our hydrodynamical models.

As predicted by the analytic scalings in Section 4.1, the models with CSM are significantly more luminous at peak (by a factor of roughly 100) than the hydrostatic model. The CSM interaction models also show delayed time of peak brightness, modified colors, and light curve shapes, as we discuss in what follows. Morozova et al. (2017) have recently discussed how rapid mass loss immediately pre-supernova can transform light curves from a IIP shape (at low pre-supernova mass loss rates) to the more luminous type IIL (at higher mass loss rates). This occurs when additional internal energy is added to the ejecta by shock-heating due to the CSM mass as it is swept up. Because the CSM lies outside the stellar radius, this new internal energy does not adiabatically decay as much prior to being radiated from the transient’s photosphere. As a result, the radiative efficiency of the models, $E_{\rm rad}/E$, ranges from 1.4% for the HSE model, to 8.5% for the $q=0.1$ model, to 23% for the $q=0.3$ model. We note that these radiative efficiencies are with a factor of two of those prediced by the scaling models of Section 4.1, for $n=4$ and $n=3$, respectively, equations (10) and (11).

Many features of our model light curves with merger-ejecta are similar to Morozova et al. (2017)’s model suites including dense CSM distributions of varying mass and power-law slope. In particular, elevated early “plateau” luminosities that decay down to the unperturbed plateau are representative of significant CSM at radii less than the transient’s eventual maximum photosphere radius of approximately $10^{15}$ cm. Comparing to Figure 3, we note that large masses of relatively close-in CSM are the distinguishing features of our models. The resultant light curves, therefore, have typical duration of hundreds of days like normal IIP, not the thousands of days observed in some IIn supernovae with extended CSM distributions like that observed for SN 1988Z and SN 2005ip (Smith, 2017).

Comparing the two mergers, the $q=0.3$ scenario with larger $M_{\rm CSM}$ has a later and more luminous peak, along with a higher luminosity during the plateau. Because the CSM mass is related to the merger, we find $M_{\rm CSM}/M_{\rm tot}\sim 1.5q$, see Section 3. The models of Figure 4 for a hydrostatic explosion, $q=0.1$, and $q=0.3$, thus provide a context for interpreting the apparent variations in CSM contribution. While the CSM mass plays a primary role in determining the light curve brightness, the distribution of CSM is crucial in shaping the light curves. In the case of the hydrostatic explosion (labeled HSE in Figure 4), there is no contribution from CSM interaction. The $q=0.3$ model shows a light curve that is always elevated by approximately three magnitudes above the HSE model due to CSM interaction (at $t>30$ d). By contrast, the $q=0.1$ model is significantly elevated above the HSE model only earlier in the lightcurve, and converges to the HSE plateau luminosity around $100$ d. The distinction between these cases lies in the slope of the CSM density and in the outer CSM radius.

For $q=0.3$, the CSM has $\rho\propto r^{-3}$ ($n=3$). The total radiative efficiency due to the CSM, $E_{\rm rad}/E$, equation (10), only decreases with the logarithm of the ratio of the expanding ejetca photosphere radius, seen in the right hand panels of Figure 4, to initial stellar radius (which is the base of the CSM distribution). Further, the outermost radius of the CSM at the moment of energy injection in our numerical model is $R_{0}=10^{15}$ cm. This is larger than the ejecta photosphere radius early in the light curve, and similar to the ejecta photosphere radius later in the light curve, implying that there is not a significant adiabatic degradation of the CSM contribution before light can escape from the expanding ejecta.

We can compare these trends to the $q=0.1$ model, in which case the overall CSM mass is lower and $\rho\propto r^{-4}$ ($n=4$). As the photosphere grows with time across the transient duration, the contribution of CSM interaction decreases approximately as $R_{\ast}/R_{\rm ph}$, see equation (11). Of similar importance, the outermost CSM radius at the moment of energy deposition, $R_{0}=2.6\times 10^{14}$ cm is a factor of a few less than the ejecta photosphere radius late in the plateau phase. This effect also decreases the contribution of the CSM to the late-stage light curve compared to an $n=4$ CSM of infinite extension, e.g. equation (11). As a result, the light curve converges to a similar magnitude as the plateau of the hydrostatic explosion late in the light curve.

In addition to brightening the explosion, CSM interaction modifies the object’s colors at timescales of days to weeks on which transients are typically discovered. CSM interaction yields bluer colors at time of peak (effective temperatures of several $10^{4}$ K on timescales of tens of days). The photosphere cools to more typical IIP temperatures of thousands of Kelvin only after 50 to 100 days (for the $q=0.1$ and $q=0.3$ models, respectively). These higher temperatures are directly representative of the extra internal energy injection due to shock heating of the ejecta by the CSM density distribution.

4.4 Implications of Varying Energetics and Donor Star Properties

We expect mergers between compact objects and giant stars to occur at a wide range of donor star and compact object properties because the binaries from which they form have broad distributions of mass, semi-major axis, and mass ratio. Further, the energy of the central engine is unknown, and may, in fact, vary from merger to merger. Here we explore the implications of the parameter space of merger properties on the resultant light curves.

Figure 5 shows V-band light curves for a range of models, all $q=0.1$, in which we vary energy (top panel), mass (center panel) and radius (bottom panel) around our fiducial, $30M_{\odot}$, $1000R_{\odot}$ case with $10^{51}$ erg explosion energy. In all of these cases, because $q=0.1$, the CSM density profile is roughly $\rho\propto r^{-4}$, and equation (11) predicts the approximate contribution of CSM interaction to the radiated energy.

Varying explosion energy with other properties kept fixed yields the qualitatively expected variation in light curve luminosity and duration – higher energy explosions give rise to faster ejecta, with more luminous but shorter duration transients. We note that the relative contribution of the CSM interaction, shown by the early bump in the lightcurve, decreases in the more energetic supernovae. When the explosion energy changes, one consequence is that the photosphere radius during the plateau phase changes, roughly as $R_{\rm ph}\propto E^{5/12}$ (Popov, 1993; Kasen & Woosley, 2009). For higher energies, the larger photosphere radii imply smaller contributions from CSM interaction, because the photosphere is further outside the outermost CSM radius, $R_{0}$. In Figure 5, we observe that the light curve shape transforms as the energy increases. This is reflective of the decreasing contribution of CSM interaction to the light curve as $E$ increases. Consequently, the radiative efficiency decreases from 10% for the $3\times 10^{50}$ erg explosion to 4.7% for the the $10^{52}$ erg explosion. With this smaller CSM contribution, the $10^{52}$ erg explosion light curve shows a relatively typical IIP shape, with a small, early bump due to the CSM.

Varying donor star properties, in the form of mass and radius, similarly changes light curve duration, peak brightness, and shape. More massive donor stars yield higher ejecta masses, but constant $M_{\rm CSM}/M_{\rm tot}\propto q$. At fixed energy, the ejecta velocities are lower and light curve durations are correspondingly longer. As mass varies in Figure 5, the plateau photosphere radius varies only mildly because the higher ejecta masses are balanced by lower ejecta velocities. Therefore, though these models show different characteristic timescales, they all have very similar peak magnitudes and degrees of CSM contribution to their overall radiated luminosity (total radiative efficiencies range from 11.7% to 9.0%).

Varying donor star radius changes not only the extent of the donor itself but the extent of the CSM, which extends to tens of stellar radii. This, in turn, varies the crucial ratio of maximum CSM radius to transient photosphere radius (because varying donor radius has little effect on $R_{\rm ph}$). When the donor is more compact, for example $150R_{\odot}$, the CSM extends to approximately $10^{14}$ cm, and largely affects only the early lightcurve. Progressively larger donors of $500R_{\odot}$ and $1000R_{\odot}$ scale the radial size of the CSM distribution. This scaling yields more CSM material at radii closer to the transient’s photosphere radius at later times (for example near peak), in turn implying higher radiative efficiencies, and brighter transients. For example, the $150R_{\odot}$ model has a radiative efficiency of only 1.9%, while the $500R_{\odot}$ model radiates 5.3% of the explosion energy and the $1000R_{\odot}$ model radiates 8.5% of the explosion energy.

Figure 6 summarizes the parameter space of merger-driven explosions in luminosity, timescale, effective temperature, and radiative efficiency. The majority of merger-driven explosions have radiative efficiency on the order of 10%, much higher than the hydrostatic model with no CSM. The $q=0.3$ model has even higher radiative efficiency of 25%. However, the more compact $150R_{\odot}$ donor model and the highest explosion energy model, $10^{52}$ erg, both show relatively minimal CSM-interaction features in Figure 5 and have somewhat lower radiative efficiency (because $R_{\ast}/R_{\rm ph}$ and $R_{0}/R_{\rm ph}$ are reduced, see section 4.1).

Together, Figure 6 shows that merger-driven explosions occupy a somewhat restricted phase space of luminosity and timescale. Typical models are more luminous than standard type IIP, but less luminous than super-luminous supernovae. In all of the models bearing significant CSM-interaction features, effective temperature varies systematically with time of peak brightness with longer-duration transients appearing redder and shorter-duration transients appearing bluer.

5 Population Synthesis of Merger-Driven Explosions

We use population synthesis models of stellar binary evolution to explore the statistical properties of binary systems at the time of a common envelope phase leading to merger between a compact object and a giant star’s core. We then use these models to estimate the population of observable merger-driven explosions.

5.1 Population Model

We analyze rapid population synthesis models from the Compact Object Mergers: Population, Astrophysics and Statistics (COMPAS) suite (Stevenson et al., 2017; Barrett et al., 2018; Vigna-Gómez et al., 2018). These models employ approximate stellar evolution tracks and parameterized physics in order to facilitate exploration of the statistical properties of binary stellar evolution – including rare outcomes like the formation of double compact object binaries (for a full description of the approach, see Stevenson et al., 2017; Barrett et al., 2018; Vigna-Gómez et al., 2018). In particular, we adopt the model parameters of Vigna-Gómez et al. (2018)’s “Fiducial” case and we capitalize on a recent development by Vigna-Gomez et. al. (in preparation) to record the characteristics and outcomes of all common envelope phases experienced by modeled binaries.

Initial distributions of binary properties are sampled at the zero-age main sequence (ZAMS) in COMPAS. In the models we study, the mass of the primary star is drawn from an initial mass function in the form $dN/dm\propto m^{-2.3}$ (Salpeter, 1955) with masses between $5\leq m/\rm{M_{\odot}}\leq 100$. The mass of the secondary star is then chosen from a flat distribution in mass ratio with $0.1<q_{\rm{ZAMS}}\leq 1$ (Sana et al., 2012). The initial separation is drawn from a flat-in-the-log distribution, $dN/da\propto a^{-1}$, with separations between $0.01<a_{\rm{ZAMS}}/\rm{AU}<1000$ (Öpik, 1924; Sana et al., 2012). All stars in our model population adopt solar metalicity ($Z=0.0142$). A total of $10^{6}$ binary systems are simulated.

Common envelope phases are identified by conditions for dynamically unstable mass transfer in COMPAS. When a common envelope episode occurs, an energy criterion is used to evaluate the outcome. In particular, the final change in orbital energy is related to the energy needed to unbind the giant star’s hydrogen envelope from its core, $\Delta E_{\rm orb}=-\alpha E_{\rm bind}$, where $\alpha=1$ is an efficiency parameter (Webbink, 1984). If the maximal change in orbital energy (defined on the basis of the minimal separation at which the core fills its Roche lobe) is insufficient to unbind the envelope, $|\Delta E_{\rm orb}|<\alpha|E_{\rm bind}|$, then a merger between the companion and the core is assumed to result. This scaling implies that more compact stars have higher binding energies and, for a given companion mass, are more likely to result in merger. More extend stars (nearer to the tip of their giant-branch evolution) have lower binding energies and their common envelope phases are more likely to result in envelope ejection (de Kool, 1990; Kruckow et al., 2016).

5.2 Compact Object-Core Mergers

The most common evolutionary channel leading to a compact - giant star merger and a merger-driven explosion is as follows. A binary pair in an initially relatively wide orbit evolves, likely going through a dynamically stable mass transfer from the initially more-massive star onto its companion. That initially more massive star undergoes core collapse, leaving behind either a neutron star or black hole remnant. Because neutron star kicks tend to be large in magnitude, a relatively small fraction of systems containing newly-formed neutron stars – less than 4% (Vigna-Gómez et al., 2018) – remain bound following the supernova. Of those that remain binaries, a large fraction will undergo a common envelope phase during a reverse episode of mass transfer onto the compact object, initiated by the expansion of the intially less massive companion after it completes core hydrogen fusion. This may result in either a merger or a common envelope ejection. Those that eject their envelopes may go on to form a double compact object binary, as discussed by Vigna-Gómez et al. (2018).

The population of common envelope phases involving compact objects in these models is depicted in Figure 7. In what follows, we report on and show only events that involve post-main-sequence donor stars more massive than $10M_{\odot}$. We show the distribution of these sources in the Hertzsprung-Russell Diagram (HRD) as well as in mass, radius, and mass ratio. We highlight the distinction between all common envelope phases involving compact objects (labeled CO CE), and events resulting in mergers between a neutron star and the donor core (labeled NS-Core) and a black hole and the donor core (labeled BH-Core).

A number of interesting trends emerge from these distributions. While common envelope phases occur throughout the donor star’s post-main sequence evolution, and therefore also the HRD, particular criteria are most likely to result in a merger. Merging sources tend to have the more compact radii compared to the overall distribution of common envelope phases ($T_{\rm eff}\gtrsim 10^{4}$ K). Neutron stars interact with a broad range of stellar companion masses, while black holes common envelope phases tend to involve massive $M\gtrsim 30M_{\odot}$ and thus luminous donors. Typical mass ratios of compact-object common envelope phases range from $0.02\lesssim q\lesssim 0.6$; those resulting in mergers tend to have $q\lesssim 0.2$. The upper limits of these ranges reflect the conditions of dynamical mass transfer stability and envelope ejection, respectively. Of the mergers, the black holes form the higher mass-ratio population, $0.1\lesssim q\lesssim 0.2$, while neutron stars typically have $q\lesssim 0.1$.

5.3 Event Rate

We can estimate the event rate of compact object - giant star mergers using the results of these population synthesis models. We simulate $10^{6}$ binary systems, or approximately $1.96\times 10^{7}M_{\odot}$ of binary mass. Each solar mass of modeled stars represents $3.8M_{\odot}$ of stars formed (Vigna-Gómez et al., 2018). From our models, common envelope phases involving compact objects and donors more massive than $10M_{\odot}$ occur with a frequency of $1.5\times 10^{-4}M_{\odot}^{-1}$, where the unit denotes mergers per solar mass of stars formed. Of these, approximately 22% result in mergers. Among the mergers, 49% involve black holes, while 51% involve neutron stars. The rate of black hole-core mergers and neutron-star mergers are thus each approximately $1.6\times 10^{-5}M_{\odot}^{-1}$. These events occur in a spread of ages between 3 and 40 Myr, and are thus strongly correlated with recent star formation. By comparison, core collapse supernovae occur with a frequency of $5.8\times 10^{-3}M_{\odot}^{-1}$ in the model systems. Merger-driven explosions therefore represent on the order of 0.6% of all core collapse events.

5.4 Outburst Population

Having assessed the population of donor stars and compact object companions that undergo mergers, we now extend the results of our light curve models to estimate the properties of the population of observable transients. Guided by the results of Sections 4.4 and 4.1, we note that CSM interaction is most important when the binary mass ratio is larger (yielding more merger ejecta and higher CSM mass) and when the radius is extended (yielding less adiabatic degradation of CSM-interaction energy, proportional to $R_{\ast}/R_{\rm ph}$). Comparison to the population properties in Figure 7 shows that the systems that tend to have high mass ratios and large radii are predominantly the BH-Core merger group, in which a black hole merges with its giant star companion. By contrast, the typical radii, $R_{\ast}\sim 100R_{\odot}$, and mass ratios, $q<0.1$, for the neutron star-core mergers are such that we expect less dramatic signatures of CSM interaction (see Figure 5).

To map our parameter variations onto the modeled population, we estimate the following scalings of $M_{V}$ with changing model parameters from the results of Figure 6,

[TABLE]

An important caveat is that, given our limited model parameter coverage, these numerical scalings represent the individual dependencies on binary properties about our fiducial model rather than the full parameter covariance. We will compare these luminosities to those of standard type IIP supernovae (Popov, 1993),

[TABLE]

We note that the Popov (1993) model accurately predicts the peak V-band luminosity of our hydrostatic model (Figures 4 and 6).

In Figure 8, we apply these scalings to the population of compact-object core mergers. Again we divide the population on the basis of whether a neutron star or a black hole is merging with the core. We assume that all events have $10^{51}$ erg explosion energy for the sake of this illustration. We find that CSM interaction (as predicted by equation (5.4)) brightens all merger-driven explosions relative to their hydrostatic equivalents (as estimated from equation (5.4)). The neutron star-core mergers are brighter by approximately 1 magnitude than their Popov-model equivalents. However, the black hole mergers are brightened significantly more, by approximately 3 magnitudes. In this diagram, we observe that the black hole merger-driven explosions form a distinct population more luminous than the non-CSM-interacting IIP population.

6 Discussion

6.1 Production of Supernovae-like Transients With Massive, Close CSM

It has recently become apparent that a large fraction of type II supernovae show signs of interaction with CSM of densities much larger than that implied by nominal stellar-wind mass loss. Type IIn supernovae have long been acknowledged to have CSM due to the persistent narrow lines in their spectra. This otherwise diverse class of supernovae occupies approximately 10% of the overall core collapse rate (e.g. Kiewe et al., 2012). More recently, evidence has been emerging that a majority (up to 70%) of type II supernovae show evidence of having at least $0.1M_{\odot}$ of CSM imprinted on their light curves (Morozova et al., 2018; Morozova & Stone, 2018). For example, Forster et al. (2018) has argued for systematic evidence that most type II shock breakouts are delayed by interaction with dense CSM. A shared feature of the CSM in many type IIP and IIL supernovae is that it is very close to the donor star, indicating its loss in the years immediately prior to the explosion (e.g. Ofek et al., 2013c; Smith & Arnett, 2014).

One proposed explanation for the presence of pre-supernova CSM ejection lies in the phenomenological comparison to luminous blue variable (LBV) outbursts, which are non-terminal outbursts of massive O-type stars. Though the precise cause of these outbursts remains uncertain (e.g. Justham et al., 2014), as does their potential correlation with the evolutionary trend of the core toward collapse, in at least one dramatic example, SN2009ip, both LBV outbursts and a terminal supernova were observed in the same object over the course of a decade (Smith et al., 2010; Prieto et al., 2013; Ofek et al., 2013a; Smith et al., 2014; Margutti et al., 2014; Mauerhan et al., 2014).

Another possible explanation links the CSM to the vigorous convection due to accelerating nuclear burning in the pre-supernova core. In this case, convection launches gravity waves at the interface between the convective core and an overlying radiative layer. These waves propagate through the radiative zone and dissipate near the base of the convective hydrogen envelope (Quataert & Shiode, 2012; Shiode & Quataert, 2014; Fuller, 2017). The luminosity of these dissipating waves can be highly super Eddington in the year prior to core-collapse, driving extensive mass loss (Quataert et al., 2016) or outbursts (Fuller, 2017).

An explosion driven by the merger itself also naturally links merger ejecta and CSM with the explosive fate of the star, as we have described in the preceding sections. However, a merger-driven model cannot explain the full diversity of type II supernovae or their CSM properties. In practice, some combination of these processes must be at play in order to explain the abundance and diversity of CSM observed in type II supernovae.

6.2 Comparison to Observed Supernovae

We compare our model light curves to two representative, well-studied supernovae. Photometric similarity is insufficient to demonstrate the origin of a given transient, as we discuss further in Section 6.3. In this section, we contextualize our model merger-driven explosions by showing that they bear similarities to transients already in the supernovae archives.

6.2.1 1979c

SN1979c, classified as a type IIL, was discovered in April 1979, several weeks after explosion (Mattei et al., 1979). It has been observed extensively at radio wavelengths (Weiler et al., 1986; Montes et al., 2000; Bartel & Bietenholz, 2008; Marcaide et al., 2009), and early modeling suggested a very large progenitor radius of $R\sim 6000R_{\odot}$ and CSM extending out to $R\sim 10^{5}R_{\odot}$ (Bartunov & Blinnikov, 1992).

Observations from the following 20 years gave rise to more theoretical discussion. Bartel & Bietenholz (2003) suggested that the remnant is expanding into low density CSM with $\rho\sim r^{-n}$, with $n=1.94^{+0.10}_{-0.05}$ decreasing to $n<1.5$ at larger radii. Later Kasen & Bildsten (2010) suggested that the light curve was brightened by the spindown of a magnetar at the center of the SN remnant. Patnaude et al. (2011) note that, rather than a magnetar, a $5-10M_{\odot}$ BH accreting from fallback material can also explain X-ray data seen from 1995 and 2007.

In the top panel of figure 9 we have plotted optical data from the first 100 days of observation (from the open supernova catalog; Guillochon et al., 2017). On top we plot absolute magnitudes from our simulations with $M_{\ast}=30M_{\odot}$ and $R_{\ast}=1000R_{\odot}$ and $E=3\times 10^{51}$erg. The plot is not a fit, but shows that the outcome of our simulations can closely replicate observed transients. This, in addition to the potential for a remnant black hole (Patnaude et al., 2011) make SN1979c an interesting candidate for further investigation under the merger-driven hypothesis.

6.2.2 1998s

SN1998s is one of the most studied type IIn supernovae (Shivvers et al., 2015). From spectral lines, two shells of CSM were identified. Fassia et al. (2001) found that the inner CSM was within 90AU from the center and the outer CSM extended from 185 AU to over 1800 AU. The classification of SN1998s as a IIn is a direct result of the very early spectral observations; the narrow line features disappeared and morphed into the broad lines of a type IIL or IIb within weeks Smith (2017).

SN1998S has later been interpreted as having a red supergiant progenitor with possibly asymmetric CSM consisting of the two separate shells, caused by separate mass-loss events. (Kangas et al., 2016) claims that this type of SN is very common and that many IIL and IIP share spectral features with IIn in early spectra.

In the bottom panel of figure 9 we plot the $V,R$ and $I$ band data from SN1998s (data from the open supernova catalog; Guillochon et al., 2017) plotted on top of absolute magnitude from simulations with $M_{\ast}=30M_{\odot}$ and $E=1\times 10^{51}$ erg, solid lines for radius scaled to $R_{\ast}=1000R_{\odot}$, and dashed lines for $R_{\ast}=500R_{\odot}$. The light curve we see from from our simulations with $R_{\ast}=1000R_{\odot}$ is similar to SN1998s, though the rate of decline perhaps fits better with our $R_{\ast}=500R_{\odot}$ simulation. Just as with SN1979c, the overall light curve shape, duration, and brightness are well-approximated by our models. The presence of nearby CSM is also consistent with a merger-driven explosion. However, the explanation for two distinct shells of CSM is not immediately apparent given our model predictions, and may be in tension with the merger-driven hypothesis for this transient (though see the discussion of Clayton et al., 2017).

6.3 Identification in Optical Surveys

Having shown that merger-driven explosion models can reproduce the basic, photometric properties of several observed supernovae, we now focus on the prospects for their more secure identification.

The prevalence of merger-driven explosions (of order 0.5% of the core-collapse rate) begs questions about their prior detection in existing datasets and their imprints on future surveys. Current surveys, such as the Zwicky Transient Factory (Bellm & Kulkarni, 2017) and All-Sky Automated Survey for Supernovae (e.g. Holoien et al., 2019) are presently discovering hundreds of new core-collapse supernovae per year. This discovery rate suggests that one or more merger-driven explosions is currently being discovered per year. Efforts at early discovery and spectroscopy of these transients aim to reveal CSM properties through “flash spectroscopy” in which the CSM is ionized prior to being swept up by the blast wave. The Large Synoptic Survey Telescope (LSST) will discover on the order of $10^{5}$ core collapse events per year (LSST Science Collaboration et al., 2009, chapter 8), implying hundreds to thousands of merger-driven explosions detected in a given observing year.

Among this flood of optical transients, the challenge will be unambiguous identification of merger-driven explosions rather than detection. A full consideration is beyond the scope of our initial study, but we speculate on several potential signatures here. As discussed in Section 3, the ejecta from the pre-merger common envelope phase are densest in the equatorial plane of the binary. When the supernova explodes into this aspherical density distribution, the blast wave will be shaped by these asymmetrical surroundings (Blondin et al., 1996). Emission from the photosphere will, as a result, be polarized by one to several percent, as has been described in the case of SN2009ip (Mauerhan et al., 2014).

The interacting binary progenitor of the explosion may also offer clues in the identification of merger-driven supernovae, as in ongoing progenitor-monitoring efforts described by Kochanek et al. (2008); Adams et al. (2017). Drawing parallels to low-mass, Galactic stellar merger events like V1309 Sco (Mason et al., 2010; Tylenda et al., 2011), increasing rates of non-conservative mass transfer (seen in the panels of Figure 1) may enshroud the merging binary in dust and cause an optical fading of the progenitor star prior to merger. In V1309 Sco, such a phase of optical dimming was observed in the phase of 100 to 1000 orbital periods prior to coalescence. In the last orbits leading into the merger (the portion captured by Figure 1), V1309 Sco brightened in optical bands as more-and-more emission arose from the outflow from the binary (Pejcha, 2014; Pejcha et al., 2016a, b, 2017). Future work is needed to extend these scenarios to detailed predictions for pre-explosive behavior in massive star coalescence.

Multiwavelength, particularly X-ray, signatures, while less frequently available than optical photometry, provide a powerful tool for probing early CSM interaction (e.g. Chevalier & Irwin, 2012; Margutti et al., 2017; Morozova & Stone, 2018). These data can probe the CSM distribuiton in great detail, including the density distribution through the hard to soft emission ratio (Morozova & Stone, 2018). If the CSM is as steep as predicted in the merger-driven models (steeper than $\rho\propto r^{-3}$), it will accelerate the leading edge of the ejecta to high velocities and produce hard x-ray emission (Morozova & Stone, 2018).

Finally, merger-driven explosions will leave a black hole as the remnant of the rapid accretion phase following merger of the compact object with the stellar core. Though black hole formation is common in core-collapse events, it is typically believed to accompany implosion rather than explosions and luminous supernovae (e.g. Sukhbold et al., 2018). If detected, the coexistence of a supernova-like transient and a remnant black hole would thus be consistent with the merger-driven explosion scenario. In theory we might distinguish neutron star and black hole central x-ray sources on the basis of their x-ray specta. In practice, this identification can be ambiguous when the surrounding, absorbing medium is substantial. One such example of a transient harboring an embedded x-ray source is AT 2018cow (Margutti et al., 2019).

7 Summary and Conclusion

In this paper, we have presented models for merger-driven explosions that arise from the plunge of a compact object within the helium core of its giant star companion following a common envelope phase (Chevalier, 2012). When a compact object merges with the helium core of a massive, post main-sequence star, the conditions for rapid, neutrino-cooled accretion are met (e.g. Zhang & Fryer, 2001). The accompanying release of energy may deposit approximately $10^{51}$ erg into the surrounding hydrogen envelope, leading to a merger-driven explosion (Chevalier, 2012). Some key findings of our investigation are:

The binary coalescence leading to the merger of the compact object with the core expels slow-moving material into the surrounding environment, forming a dense, toroidal CSM (Figures 1 and 2). The spherically-averaged density profile has a steep radial slope of $\rho\propto r^{-3}$ or $\rho\propto r^{-4}$ (Figure 3). 2. 2.

Using 1D radiation hydrodynamic models of the explosions, we find that the CSM distribution is crucial in shaping the transient light curves. Merger-driven explosions are brightened by up to three magnitudes relative to their counterparts in hydrostatic stars (Figure 4), with timescale and light curve shape that vary with donor-star mass and radius and explosion energy (Figures 5 and 6). 3. 3.

From population models, we find that black hole and neutron star mergers with giant star companions occur with similar frequency, each with a rate per mass of stars formed of $1.6\times 10^{-5}M_{\odot}^{-1}$. The combined rate is 0.6% of the core-collapse rate in our models. Merger-driven explosions occur across a roughly flat distribution of donor-star masses from $10M_{\odot}$ to $100M_{\odot}$ (Figure 7). CSM interaction brightens neutron star mergers by approximately one magnitude, but brightens the population of black hole mergers by approximately three magnitudes relative to type IIP models with the same energy injection and pre-supernova stellar mass and radius (Figure 8). 4. 4.

The most luminous transients, those involving black hole mergers, are at least as common as their less luminous neutron star counterparts. Black hole mergers have $M_{V,{\rm peak}}\sim-18$ to $-19$ with $t_{\rm peak}\sim 20$ to $30$ d. The implication for optical surveys is that the brightest, easiest-to-detect events comprise a significant fraction of the entire population.

The calculations presented in this paper have demonstrated that merger-driven explosions provide a natural mechanism for the production of supernovae-like transients with close-in, slow-moving CSM. Future work could improve on the treatment of the stellar model (a polytropic envelope in our approximation) and the details of energy injection into this envelope. At present, we inject energy spherically into the envelope at one tenth the star’s overall radius. In practice, the unknown location and asymmetry of energy injection might play a key role in shaping transient light curves, colors, and peak luminosities with respect to the estimates of our current models.

We compare our models to two representative supernovae, SN1979c and SN1998s in Figure 9. However, we note that more work is needed to provide unambiguous confirmations of merger-driven explosions. In Section 6.3, we discuss additional strategies for the identification of merger-driven explosions including their asymmetry and polarization due to the toroidal CSM, the properties of their progenitor binaries, and their early spectra and X-ray emission. In future work, these signatures can be investigated through multi-dimensional calculations of the explosive evolution and emergent light curve, as well as more detailed modeling of the progenitor system’s plunge toward merger.

We thank R. Margutti, E. Ramirez-Ruiz, and M. Rees for advice and helpful discussions in the development of this work. We gratefully acknowledge the support of E. Ostriker and J. Stone in the development and analysis of the hydrodynamic models, and V. Morozova for support with SNEC. S. S. acknowledges support by the Danish National Research Foundation (DNRF132). Part of the simulations used in this paper were performed on the University of Copenhagen high-performance computing cluster funded by a grant from VILLUM FONDEN (project number 16599). M.M. is grateful for support for this work provided by NASA through Einstein Postdoctoral Fellowship grant number PF6-170169 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. This work was supported in part by the Black Hole Initiative at Harvard University, which is funded by a JTF grant. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center.

Appendix A Validation of Light Curve Calculations

In this appendix, we discuss the validation of several numerical choices in the 1D radiation hydrodynamics calculations with SNEC that we use to produce model light curves.

In mapping the 3D hydrodynamics calculation of the merger (Section 3) to the 1D explosive calculation, we need to make an assumption about the location (described by radius or mass coordinate) where the explosion energy is deposited. We have argued in Section 2 that this deposition location is somewhere within the hydrogen envelope. Here we explore the sensitivity to that choice as follows. Beginning with our fiducial case of a $30M_{\odot}$ donor star and $3M_{\odot}$ black hole in a $q=0.1$ merger, we deposit $10^{51}$ erg of thermal energy spread over $0.1M_{\odot}$ at different mass coordinate locations. Our default assumption is $M_{\rm in}=10.75M_{\odot}$, which corresponds to the enclosed mass of at $0.1R_{\ast}$ of $7.75M_{\odot}$ plus a $3M_{\odot}$ black hole. This model in Figure 10 corresponds to the fiducial simulation presented in Figure 4, which is labeled $q=0.1$. We then vary the mass coordinate at which thermal energy is deposited, moving outward in the star’s Hydrogen envelope. Material inside $M_{\rm in}$ acts as a gravitational point mass for the remainder of the calculation. We find that for $M_{\rm in}=10.75M_{\odot}$, $M_{\rm in}=14M_{\odot}$, and $M_{\rm in}=18M_{\odot}$ (corresponding to radius coordinates of 0.1, 0.16, and 0.24 times the donor star’s original radius) the model light curves are very similar indicating only weak dependence on how the energy is spatially deposited within the hydrogen envelope. All of these radii are significantly outside the star’s more compact Helium core. We note that for $M_{\rm in}>20M_{\odot}$, the case in which $>2/3$ of the donor star forms a black hole while only $<1/3$ is expelled, we do observe departures in the model light curves, with the bulk of the thermal energy radiated early, and not coupling efficiently to driving envelope expansion.

We also test the dependence of our model results on spatial resolution within the 1D SNEC calculations. Our fiducial case divides the mass into 456 elements. Figure 11 compares this case to models with twice and four times as many zones (912, 1824, respectively). These tests confirm that our results are converged to within 1% across the light curve duration with any of these resolution choices.