Introduction

Despite being very rare, massive stars such as luminous blue variable (LBVs), red super giants (RSGs), and Wolf-Rayet stars (WRs) play a pivotal role in enriching the interstellar medium (ISM) through mass loss, and they are an important source of heavier elements contributing to the chemical enrichment of galaxies (e.g., Maeder1981). The deaths of these massive stars are associated with some of the highest energy phenomena in the universe such as core-collapse supernovae (ccSNe, Smartt2009), long-duration gamma-ray bursts (e.g., Stanek et al.2003), neutrino bursts (e.g., Bionta et al.1987) and gravitational wave bursts (e.g., Ott2009). The physical mechanism, energetics and observed properties of these events depend on the structure and terminal mass of the evolved stars at core-collapse, which in turn are determined by stellar mass loss (see, e.g., review by Smith2014). In addition, there is also evidence that some supernova (SN) progenitors undergo major mass ejection events shortly before exploding (e.g., Smith et al.2008; Ofek et al.2013; Gal-Yam et al.2007), further altering the properties of the explosion and implying a connection between some eruptive mass-loss events and death. It is generally agreed that the effects of winds are metallicity dependent (e.g., Heger et al.2003; Meynet et al.1994) and the SNe requiring a very dense circumstellar medium (e.g., Schlegel1990; Filippenko1997) predominantly occur in lower metallicity galaxies (e.g., Stoll et al.2011). This strongly suggests that the nature and distribution of stars undergoing impulsive mass loss will also be metallicity dependent and a full understanding requires exploring galaxies beyond the Milky Way.

Understanding the evolution of massive (M$\gtrsim30$M$_\odot$) stars is challenging even when mass loss is restricted to continuous winds (e.g., Fullerton et al.2006). However, shorter, episodic eruptions, rather than steady winds, may be the dominant mass loss mechanism in the tumultuous evolutionary stages toward the end of the lives of the most massive stars (e.g., Humphreys & Davidson1984; Smith & Owocki2006) as they undergo periods of photospheric instabilities leading to stellar transients ( $M_V\lesssim-13$) followed by rapid ( $\dot{M}\gtrsim 10^{-4}\,M_\odot$year) mass-loss in the last stages of their evolution (see Smith2014; Kochanek et al.2012). Deciphering the rate of these eruptions and their consequences is challenging because no true analog of $\eta$Car in mass, luminosity, energetics, mass lost and age has been found (see Smith et al.2011; Kochanek et al.2012), and the associated transients are significantly fainter than supernova explosions and are easily missed. These phases are as difficult to model theoretically as they are to simulate computationally.

Dense winds tend to form dust, although for hot stars the wind must be dense enough to form a pseudo-photosphere in the wind (Davidson1987) that shields the dust formation region from the UV emission of the star (Kochanek2011). The star will then be heavily obscured by dust for an extended period after the eruption (see, e.g., Humphreys & Davidson1994). The Great Eruption of $\eta$Car between 1840 and 1860 is the most studied case of a stellar outburst (see, e.g., Humphreys et al.2012). The $\sim10 M_{\odot}$ ejecta are now seen as a dusty nebula around the star absorbing and then reradiating $\sim90\%$ of the light in the mid-IR. This means that dusty ejecta are a powerful and long-lived signature of eruption. The emission from these dusty envelopes peaks in the mid-IR with a characteristic red color and a rising or flat spectral energy distribution (SED) in the Spitzer IRAC (Fazio et al.2004) bands.

In the Galaxy, stars with resolved shells of dust emission are easily found at 24$\micron$ (Gvaramadze et al.2010; Wachter et al.2010). The advantage of the 24$\micron$ band is that it can be used to identify dusty ejecta up to $10^3 - 10^4$years after formation. A minority of these objects are very luminous stars (L $\gtrsim10^{5.5}\,$L$_\odot$) with massive ($\sim0.1-10\,$M$_\odot$) shells (see summaries by Humphreys & Davidson1994; Smith2009; Smith & Owocki2006; Humphreys et al.1999; Vink et al.2009). These include AGCar (Voors et al.2000), the Pistol Star (Figer et al.1999), G79.29$+$0.46 (Higgs et al.1994), Wray17$-$96 (Egan et al.2002), and IRAS18576$+$0341 (Ueta et al.2001). These systems are significantly older ($10^3 - 10^4$years) than $\eta$Car, which makes it difficult to use the ejecta to probe the rate or mechanism of mass-loss. Still, the abundance of Galactic shells implies that the rate of $\eta$Car-like eruptions is on the order of a modest fraction of the ccSN rate (Kochanek2011). Their emission peaks in the shorter IRAC bands when they are relatively young ($\sim10-100$years) because the dust becomes cooler and the emission shifts to longer wavelengths as the ejected material expands (Kochanek et al.2012). It is difficult to quantify searches for such objects in our Galaxy because it is hard to determine the distances and the survey volume because we have to look through the crowded and dusty disk of the Galaxy. Surveys of nearby galaxies are both better defined and can be used to build larger samples of younger systems whose evolution can be studied to better understand the mechanism. We previously demonstrated in Khan et al. (2010); Khan et al. (2011) that it is possible to identify post-eruptive massive stars in galaxies beyond the Local Group using the mid-IR excess created by warm circumstellar dust despite the crowding problems created by the limited spatial resolution of Spitzer at greater distances.

In Khan et al. (2013) (``PaperI'' hereafter) we used archival Spitzer IRAC images of seven $\lesssim4$Mpc galaxies (closest to farthest: NGC6822, M33, NGC300, NGC2403, M81, NGC0247, NGC7793) in a pilot study to search for extragalactic analogs of $\eta$Car. We found 34 candidates with flat or rising mid-IR spectral energy distributions (SEDs) and total mid-IR luminosity $L_{mIR}\gtrsim 5\times10^5 L_\odot$. Here, in PaperII, we characterize these sources and quantify the rate of episodic mass loss from massive stars in the last stages of evolution. First, we construct extended optical through far-IR SEDs using archival HST, 2MASS, and Herschel data as well as ground based data (Section2). Then, we classify the sources as either stellar or non-stellar based on properties of the extended SEDs and model the SEDs to infer the properties of the underlying star and the obscuring circumstellar medium (Section3). Next, we relate these properties to the observed ccSN rate of the targeted galaxies to quantify the rate of episodic mass loss in the last stages of massive star evolution (Section4). Finally, we consider the implications of our findings for theories and observations of massive star evolution and their fates (Section5).