FINDING η CAR ANALOGS IN NEARBY GALAXIES USING Spitzer. II. IDENTIFICATION OF AN EMERGING CLASS OF EXTRAGALACTIC SELF-OBSCURED STARS

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 through mass loss, and they are an important source of heavier elements contributing to the chemical enrichment of galaxies (e.g., Maeder 1981). The deaths of these massive stars are associated with some of the highest energy phenomena in the universe such as core-collapse supernovae (ccSNe; e.g., Smartt 2009), 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., Ott 2009). 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 Smith 2014). In addition, there is also evidence that some supernova (SN) progenitors undergo major mass ejection events shortly before exploding (e.g., Gal-Yam et al. 2007; Smith et al. 2008; Ofek et al. 2013), 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., Meynet et al. 1994; Heger et al. 2003) and the SNe requiring a very dense circumstellar medium (e.g., Schlegel 1990; Filippenko 1997) 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. A full understanding of the pre-SN evolution of these stars requires exploring galaxies beyond the Milky Way (see, e.g., Humphreys et al. 2013, 2014 for recent examples of such studies).

Understanding the evolution of massive (M ≳ 30 M_☉) 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 & Davidson 1984; Smith & Owocki 2006) as they undergo periods of photospheric instabilities leading to stellar transients (M_V ≲ −13) followed by rapid ($\dot{M}\gtrsim 10^{-3}\,M_\odot$ yr) mass loss in the last stages of their evolution (e.g., Kochanek et al. 2012; Smith 2014). Dense winds and massive outflows tend to form dust, although for hot stars the wind must be dense enough to form a pseudo-photosphere in the wind (Davidson 1987) that shields the dust formation region from the UV emission of the star (Kochanek 2011). The star will then be heavily obscured by dust for an extended period after the eruption (see, e.g., Humphreys & Davidson 1994). 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. The production and lifetime of the mid-IR component may also be driven by the interaction between explosive stellar ejecta and pre-existing circumstellar material (Smith 2013).

A possible Galactic analog of these extragalactic events is the Great Eruption of η Carinae (η Car) in the mid-1800 s. η Car is one of the most massive (100–150 M_☉) and most luminous (∼3 × 10⁶ L_☉) stars in our Galaxy. Following its Great Eruption between 1840 and 1850, η Car faded dramatically, and only began to become steadily brighter in ∼1950 (see, e.g., Humphreys et al. 2012). The Great Eruption was a period of enormously greater mass loss whose origin could be radiative (e.g., Smith & Owocki 2006) or explosive (e.g., Smith 2013). Whatever the cause, η Car began forming dust shortly after the eruption, and is now surrounded by roughly 10 M_☉ of dusty ejecta that absorbs ∼90% of the stellar radiation and re-emits it in the mid-IR with an SED peaking near 24 μm. For these first few 100 yr, the dust is warm enough to make η Car a very bright source in the warm Spitzer bands (3.6 and 4.5 μm). As the ejecta expands, the dust cools, and the emission in the warm Spitzer bands rapidly drops. η Car will remain a very luminous source at longer wavelengths (e.g., 24 μm) for millenia.

In fact, Galactic stars with resolved shells of dust formation are easily found at 24 μm (Wachter et al. 2010; Gvaramadze et al. 2010). The advantage of the 24 μm band is that it can be used to identify dusty ejecta up to 10³–10⁴ yr after formation. A minority of these objects are very luminous stars (L  ≳ 10^5.5 L_☉) with massive (∼0.1–10 M_☉) shells (see summaries by Humphreys & Davidson 1994; Humphreys et al. 1999; Smith & Owocki 2006; Smith 2009; Vink et al. 2012). These include AG Car (Voors et al. 2000), the Pistol Star (Figer et al. 1999); G79.29+0.46 (Higgs et al. 1994); Wray 17−96 (Egan et al. 2002); IRAS 18576+0341 (Ueta et al. 2001); and WR 122 (Crowther & Smith 1999). These systems are significantly older (10³–10⁴ yr) than η 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 η Car-like eruptions is on the order of a modest fraction of the ccSN rate (Kochanek 2011).

While ongoing studies are helping us further analyze the Great Eruption (see, e.g., Rest et al. 2012; Prieto et al. 2014), deciphering the rate of such events and their consequences is challenging because no true analog of η Car in mass, luminosity, energetics, mass lost or age has been found (see, e.g., Smith et al. 2011; Kochanek et al. 2012). Moreover, the associated transients are significantly fainter than SN explosions (see, for example, Mauerhan et al. 2014) and can be easily missed. It is difficult to quantify searches for such objects in our Galaxy because it is hard to determine the distances and the survey volume as 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, 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; hereafter Paper I) we used archival Spitzer IRAC images of seven ≲ 4 Mpc galaxies (closest to farthest: NGC 6822, M 33, NGC 300, NGC 2403, M 81, NGC 0247, NGC 7793) in a pilot study to search for extragalactic analogs of η Car. We found 34 candidates with flat or rising mid-IR SEDs and total mid-IR luminosity L_mIR ≳ 5 × 10⁵L_☉. Here in Paper II, 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 Hubble Space Telescope (HST), Two Micron Sky Survey (2MASS), and Herschel data as well as ground-based data (Section 2). 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 (Section 3). 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 (Section 4). Finally, we consider the implications of our findings for theories and observations of massive star evolution and their fates (Section 5).

In this section, we describe the details of how we obtained the photometric measurements at various wavelengths to determine the properties of the candidates from Paper I. The optical through far-IR photometry are reported in Table 1, and the extended SEDs are shown in Figures 1 and 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We utilized VizieR³ (Ochsenbein et al. 2000) to search for other observations of the candidates, in particular for WISE (Wright et al. 2010, 12 μm), 2MASS (Cutri et al. 2003, JHK_s), Sloan Digital Sky Survey (SDSS; Abazajian et al. 2009, ugriz) and X-ray detections. For M 33, we used the UBVRI images from the Massey et al. (2006) optical survey, and archival HST images of NGC 300, NGC 2403, M 81, NGC 247, and NGC 7793. Finally, we used Herschel PACS data to supplement the Spitzer measurements.

For the Spitzer IRAC 3.6, 4.5, 5.8, and 8 μm as well as MIPS (Rieke et al. 2004) 24, 70, and 160 μm data, we use the measurements reported in Paper I. For M 33, our measurements were based on IRAC data from McQuinn et al. (2007) and MIPS data from the Spitzer Heritage Archive.⁴ Data from the LVL survey (Dale et al. 2009) were used for NGC 300 and NGC 247, and data from the SINGS survey (Kennicutt et al. 2003) for NGC 6822, NGC 2403, and M 81.

We used the Herschel PACS (Poglitsch et al. 2010) 70, 100, and 160 μm images available from the public Herschel Science Archive.⁵ Although both MIPS and PACS cover the same far-IR wavelength range (70–160 μm), Herschel has significantly higher resolution (see Figure 3). All three PACS band data were available for M 33 and NGC 7793, 70 and 160 μm data were available for NGC 2403 and M 81, and 100 and 160 μm data were available for NGC 300. There are no publicly available PACS images of the candidates in NGC 247. We used aperture photometry (IRAF⁶ ApPhot/Phot) with the extraction apertures and aperture corrections from Balog et al. (2014) and given in Table 2. As with our treatment of the MIPS 70 and 160 μm measurements in Paper I, we treat the measurements obtained in the PACS bands as upper limits because the spatial resolution of these bands requires increasingly large apertures at longer wavelengths. For similar reasons, we also treat the WISE 12 μm fluxes, where available, as upper limits.

Zoom In Zoom Out Reset image size

Table 2. PACS Aperture Definitions

Band	Pixel Scale	Rap	Rin	Rout	Ap. Corr.
(μm)
70 μm	32	64	608	704	0.72−1
100 μm	32	64	608	704	0.69−1
160 μm	64	128	1216	1408	0.78−1

Download table as:  ASCIITypeset image

For the optical photometry of the candidates in M 33, we used the Local Group Galaxies Survey UBVRI images (Massey et al. 2006). First we verified that the coordinates match with the IRAC images to within few× 01 and then used 10 radius extraction apertures centered on the IRAC source locations. We transformed the aperture fluxes to Vega-calibrated magnitudes using zero point offsets determined from the difference between our aperture magnitudes and calibrated magnitudes for bright stars in the Massey et al. (2006) catalog of M 33.

For the candidates in NGC 300, M 81, NGC 2403, and NGC 247, we searched the ACS Nearby Galaxy Survey (Dalcanton et al. 2009) B, V, and (where available) I-band point source catalogs derived using DOLPHOT (Dolphin 2000). We verified that the IRAC and HST astrometry of the NGC 300, NGC 2403 and NGC 247 images agree within (mostly) ≲ 01 to (in a few cases) 03. We corrected the astrometry of the M 81 HST images using the LBT images described later in this section to achieve similar astrometric accuracy. We also used the HST I-band photometry of M 81 from HST program GO-10250 (P.I. J. Huchra). We retrieved all publicly available archival HST images of NGC 7793 overlapping the IRAC source locations along with the associated photometry tables from the Hubble Legacy Archive.⁷ The HST and Spitzer images have a significant (few× 10) astrometric mismatch, and there are too few reference stars in the HST images to adequately improve the astrometry. Therefore, we utilized the IRAF GEOXYMAP and GEOXYTRAN tasks to locally match the overlapping HST and Spitzer images of NGC 7793 within uncertainties of 01 ∼ 03.

We have variability data for the galaxies M 81 and NGC 2403 from a Large Binocular Telescope survey in the UBVR bands that is searching for failed SNe (Kochanek et al. 2008; Gerke et al. 2014), studying SN progenitors and impostors (Szczygieł et al. 2012), and Cepheid variables (Gerke et al. 2011; M. Fausnaugh et al. (2014, in preparation). We analyzed 27 epochs of data for M 81 and 28 epochs of data for NGC 2403, spanning a 5 yr period. The images were analyzed with the ISIS image subtraction package (Alard & Lupton 1998; Alard 2000) to produce light curves (see Figure 4).

Zoom In Zoom Out Reset image size

In this section, we first discuss how we classify the candidates based on their SEDs. Next, we describe the non-stellar and stellar sources. Finally, we model the SEDs of the stellar sources to determine their physical properties. Figure 1 shows the SEDs of the stellar sources with the best-fit SED models over plotted and Figure 2 shows the SEDs of the non-stellar sources.

3.1. Source Classification

We classify the candidates either as stellar or non-stellar based on their photometric properties. We focus on identifying two tell-tale signatures of the SED of a luminous star obscured by warm circumstellar dust—low optical fluxes or flux limits compared to the mid-IR luminosities and signs of the SEDs turning over between 8 μm and 24 μm. Toward longer wavelengths, emission from warm circumstellar dust should peak between the IRAC 8 μm and MIPS 24 μm bands. It is almost impossible for mass lost from a single star to both have a significant optical depth and a dust temperature cold enough to peak at wavelengths longer than ∼24 μm. Such systems are almost certainly star clusters with significant amounts of cold dust. Therefore, any SED that appears to have a steep slope between 8 and 24 μm is considered to be a likely cluster, rather than a single dust obscured star. Frequently, these sources are also too luminous to be a single star. At the shorter wavelengths, we expect a dusty star to have relatively lower luminosity compared to its mid-IR luminosity and redder optical colors.

We examine the HST B − V/V and the V − I/V color–magnitude diagrams (CMDs) for each source for which HST data is available. The presence of a very red optical counterpart or the absence of a luminous star supports the existence of significant dust obscuration. On the other hand, the presence of a blue or bright optical counterpart makes it likely that the source is a star cluster, a background galaxy/active galactic nucleus (AGN), or a foreground star. We first search for bright and/or red optical sources within the 03 matching radius that can be the obvious counterpart of the bright and red IRAC source. Next, if multiple bright and/or red optical matches are found, we identify the best astrometric match to the IRAC location. Finally, if no reasonable match is found, we adopt the flux of the brightest of the nearby sources as a conservative upper limit on the optical luminosity of the candidate.

To demonstrate these, we discuss the case of M 81-12 in detail. M 81-12 has a steeply rising optical and mid-IR SED (Figure 5) with two distinct peaks—one in the near-IR, between the R band and 3.6 μm, and another in the mid-IR between 8 and 24 μm. Figure 6 shows the HST optical CMD for sources near the location of M 81-12. Besides the sources within the 03 matching radius, it also shows all sources within 03–20 of the candidate using a different symbol to emphasize the absence of any other unusual nearby sources. We detect a very red (B = 23.95, V = 21.98, I = 19.07, B − V ≃ 2, V − I ≃ 2.9) HST counterpart with an excellent astrometric match (<01, Figure 7) to the IRAC position. This source is the brightest, red HST point source within 20 of the IRAC location (Figure 6) and so we define it to be the counterpart used in the SED. The LBT V- and R-band light curves show a variable source with the correlated irregular variability (∼0.4 mag, Figure 4) typical of many evolved massive stars (e.g., Kourniotis et al. 2014). Based on the SED shape and the unambiguous detection of a red, variable optical counterpart, we conclude that M 81-12 is a massive, dust-obscured, single star.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In addition to Object X (M 33-1), we identified 17 additional dust obscured stars and classified 16 others as non-stellar. We left one source (N 7793-12) unclassified due to a lack of sufficient optical data (it falls on an HST/Advanced Camera for Surveys (ACS) chip gap). It could well be a dusty star, but we do not discuss it further.

3.2. The 18 Stars and 16 Non-stellar Sources

3.3. SED Modeling

We fit the SEDs of the 18 self-obscured stars using DUSTY (Ivezic & Elitzur 1997; Ivezic et al. 1999; Elitzur & Ivezić 2001) to model radiation transfer through a spherical dusty medium surrounding a star and Figure 1 shows the best-fit models. We estimate the properties of a blackbody source obscured by a surrounding dusty shell that would produce the best-fit to the observed SED (see Figure 8 for an example). We considered models with either graphitic or silicate (Draine & Lee 1984) dust. We distributed the dust in a shell with a ρ∝1/r² density distribution. The models are defined by the stellar luminosity (L_*), stellar temperature (T_*), the total (absorption plus scattering) V-band optical depth (τ_V), the dust temperature at the inner edge of the dust distribution (T_d), and the shell thickness ζ = R_out/R_in. The exact value of ζ has little effect on the results, and after a series of experiments with 1 < ζ < 10, we fixed ζ = 4 for the final results. We embedded DUSTY inside a Markov Chain Monte Carlo driver to fit each SED by varying T_*, τ_V, and T_d. We limit T_* to a maximum value of 30,000 K to exclude unrealistic temperature regimes.

Zoom In Zoom Out Reset image size

The parameters of the best-fit model determine the radius of the inner edge of the dust distribution (R_in). The mass of the shell is

$$\begin{equation} M_e = \frac{4 \pi R_{\rm in}^2 \tau _V}{\kappa _V}, \end{equation} \tag{ 1 }$$

where we simply scale the mass for a V-band dust opacity of κ_V = 100 κ₁₀₀ cm² g⁻¹ and the result can be rescaled for other choices as $M_e \propto \kappa _V^{-1}$. Despite using a finite width shell, we focus on R_in because it is well-constrained while R_out (or ζ) is not. We can also estimate an age for the shell as

$$\begin{equation} t_e = \frac{R_{\rm in}}{v_e}, \end{equation} \tag{ 2 }$$

where we scale the results to v_e = 100 v_e100 km s⁻¹.

For a comparison sample, we followed the same procedures for the SEDs of three well-studied dust obscured stars: η Car (Humphreys & Davidson 1994); the Galactic OH/IR star IRC+10420 (Jones et al. 1993; Humphreys et al. 1997; Tiffany et al. 2010); and M 33's Variable A, which had a brief period of high mass loss leading to dust obscuration over the last ∼50 yr (Hubble & Sandage 1953; Humphreys et al. 1987, 2006). We use the same SEDs for these stars as in Khan et al. (2013). In Table 3, we report χ², τ_V, T_d, T_*, R_in, L_*, M_e (Equation (1)), and t_e (Equation (2)) for the best-fit models for these three sources as well as the newly identified stars. The stellar luminosities required for both dust types are mutually consistent because the optically thick dust shell acts as a calorimeter. However, because the stars are heavily obscured and we have limited optical/near-IR SEDs, the stellar temperatures generally are not well constrained. In some cases, for different dust types, equally good models can be obtained for either a hot (>25, 000 K, such as a LBV in quiescence) or a relatively cooler (<10, 000 K, such as a LBV in outburst) star. Indeed, for many of our 18 sources, the best fit is near the fixed upper limit of T_* = 30, 000 K. To address this issue, we also tabulated the models on a grid of three fixed stellar temperatures, T_* = 5000 K, 7500 K, 20, 000 K, for each dust type. The resulting best-fit parameters are reported in Tables 4 and 5.

Table 3. Best-fit Models

			Graphitic										Silicate
ID			χ2	τV	Td	T*	log (Rin)	log L*	Me	te			χ2	τV	Td	T*	log (Rin)	log L*	Me	te
					(K)	(K)	(cm)	(L☉)	(M☉)	(yr)					(K)	(K)	(cm)	(L☉)	(M☉)	(yr)
M 33 − 1			46	8.46	708	6749	16.04	5.60	0.026	34.79			18	10.70	1218	24602	15.80	5.68	0.064	19.8
M 33 − 3			34	1.16	416	29927	16.94	5.81	0.096	278.5			55	2.63	649	29955	16.38	5.88	0.565	76.21
M 33 − 4			47	2.00	390	29977	16.99	5.76	0.156	307.3			75	3.88	599	29991	16.40	5.80	1.182	80.14
M 33 − 7			29	2.67	381	29931	17.06	5.85	0.222	360.5			57	4.79	599	29973	16.43	5.86	2.170	85.98
N 300 − 1			4	5.97	506	17129	16.65	5.74	0.083	141.1			8	8.57	959	28956	16.09	5.83	0.744	39.24
N 2403 − 2			30	0.90	440	29988	16.85	5.76	0.076	227			26	2.56	676	29942	16.34	5.84	0.292	68.79
N 2403 − 3			8	76.36	455	2533	16.38	5.88	0.263	76.89			2	25.27	499	4981	16.11	5.88	2.821	40.8
N 2403 − 4			100	5.39	398	3790	16.62	5.86	0.045	131.7			146	10.00	568	4545	15.93	5.89	0.585	26.78
N 2403 − 5			27	1.97	429	16146	16.85	5.84	0.109	224.6			5	3.56	662	20750	16.34	5.98	0.621	70.06
M 81 − 5			1	7.60	428	3050	16.40	5.71	0.117	78.89			2	35.00	1117	29778	15.86	5.58	0.296	23.14
M 81 − 6			0.9	8.61	504	4897	16.31	5.55	0.086	64.76			0.5	46.97	985	13885	15.73	5.49	0.226	17.15
M 81 − 11			17	2.56	463	10825	16.64	5.68	0.042	138.4			3	4.54	746	15015	16.09	5.82	0.307	38.59
M 81 − 12			12	4.31	365	5548	16.84	5.89	0.105	217.5			10	7.84	529	7910	16.16	5.92	1.275	46.34
M 81 − 14			14	20.08	341	4839	16.91	5.97	0.591	260.4			8	29.47	416	4528	16.25	5.93	8.513	56.65
N 7793 − 3			⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅			⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
N 7793 − 9			46	4.61	432	14632	16.73	5.61	0.101	169.2			27	7.04	713	22072	16.18	5.70	0.825	47.86
N 7793 − 10			32	3.80	396	7406	16.85	5.92	0.121	222.3			25	6.36	609	12175	16.24	5.99	1.174	55.18
N 7793 − 13			32	1.90	369	29943	17.15	5.99	0.377	452.2			42	4.01	546	29989	16.59	6.05	2.431	122.6
IRC+10420			222	8.06	399	8157	16.79	5.76	0.030	194.2			240	11.86	835	11780	15.80	5.55	1.900	20.17
M 33 Var A			43	2.29	536	11741	16.27	5.23	0.003	58.86			8	4.11	1046	14549	15.56	5.29	0.050	11.54
η Car			490	4.55	361	18134	17.41	6.57	5.611	809			853	7.99	468	26164	17.02	6.74	18.615	335.1

Download table as:  ASCIITypeset image

Table 4. Best-fit Models for Graphitic Dust and Fixed Temperature

		T* = 5000 K						T* = 7500 K						T* = 20, 000 K
ID		τV	Td	log (Rin)	log L*	Me		τV	Td	log (Rin)	log L*	Me		τV	Td	log (Rin)	log L*	Me
			(K)	(cm)	(L☉)	(M☉)			(K)	(cm)	(L☉)	(M☉)			(K)	(cm)	(L☉)	(M☉)
M 33 − 1		8.17	595	16.15	5.58	0.102		8.05	701	16.09	5.62	0.077		5.88	900	15.96	5.66	0.031
M 33 − 3		2.50	400	16.61	5.69	0.261		2.67	460	16.54	5.59	0.202		1.89	400	16.92	5.76	0.821
M 33 − 4		3.14	400	16.56	5.58	0.260		3.52	400	16.68	5.58	0.507		4.53	614	16.26	5.72	0.094
M 33 − 7		3.69	400	16.59	5.64	0.351		4.11	400	16.71	5.66	0.679		3.18	400	16.91	5.74	1.321
N 300 − 1		6.52	501	16.41	5.72	0.271		6.89	504	16.52	5.74	0.475		5.47	500	16.69	5.75	0.824
N 2403 − 2		3.13	401	16.63	5.72	0.358		2.99	493	16.45	5.56	0.149		1.89	408	16.88	5.73	0.683
N 2403 − 3		18.58	404	16.71	5.86	3.070		9.41	405	16.81	5.86	2.466		6.55	417	16.92	5.83	2.849
N 2403 − 4		5.91	400	16.74	5.93	1.121		7.09	300	17.32	6.33	19.435		4.92	1182	15.75	5.91	0.010
N 2403 − 5		2.78	400	16.71	5.89	0.460		2.96	400	16.82	5.87	0.812		1.55	400	16.99	5.91	0.928
M 81 − 5		14.32	514	16.32	5.58	0.393		16.83	548	16.35	5.55	0.530		19.68	587	16.41	5.54	0.817
M 81 − 6		14.55	535	16.25	5.52	0.289		16.18	559	16.30	5.50	0.405		13.09	558	16.45	5.49	0.653
M 81 − 11		2.81	406	16.66	5.80	0.369		2.67	498	16.47	5.64	0.146		1.47	432	16.85	5.78	0.463
M 81 − 12		3.77	399	16.67	5.79	0.518		3.99	400	16.78	5.79	0.910		2.59	400	16.95	5.81	1.295
M 81 − 14		21.78	350	16.91	5.97	9.040		21.43	400	16.81	5.84	5.613		22.28	430	16.90	5.83	8.834
N 7793 − 3		⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅		⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅		⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
N 7793 − 9		7.50	400	16.61	5.66	0.506		8.00	400	16.77	5.76	0.629		7.00	402	17.02	5.96	1.516
N 7793 − 10		2.73	500	16.19	5.30	0.584		2.74	500	16.27	5.25	1.030		1.12	500	16.45	5.30	1.371
N 7793 − 13		5.74	300	17.34	6.60	0.528		5.89	400	17.16	6.53	0.923		4.29	400	17.32	6.55	1.224
IRC+10420		5.08	401	16.60	5.66	0.782		5.25	424	16.64	5.63	1.743		4.00	400	16.89	5.70	4.823
M 33 Var A		3.37	400	16.72	5.90	0.041		3.59	400	16.83	5.89	0.060		2.18	400	17.00	5.92	0.056
η Car		3.67	400	16.68	5.82	17.258		3.86	400	16.79	5.81	7.732		2.45	400	16.95	5.83	11.767

Download table as:  ASCIITypeset image

Table 5. Best-fit Models for Silicate Dust and Fixed Temperature

		T* = 5000 K						T* = 7500 K						T* = 20, 000 K
ID		τV	Td	log (Rin)	log L*	Me		τV	Td	log (Rin)	log L*	Me		τV	Td	log (Rin)	log L*	Me
			(K)	(cm)	(L☉)	(M☉)			(K)	(cm)	(L☉)	(M☉)			(K)	(cm)	(L☉)	(M☉)
M 33 − 1		14.42	816	15.58	5.56	0.013		14.26	973	15.59	5.59	0.014		10.86	1204	15.75	5.67	0.022
M 33 − 3		4.30	592	15.83	5.75	0.012		4.66	654	15.88	5.68	0.017		3.36	696	16.20	5.77	0.053
M 33 − 4		5.26	513	15.90	5.65	0.021		5.64	601	15.89	5.58	0.021		2.58	400	16.87	5.66	0.892
M 33 − 7		5.91	523	15.91	5.68	0.025		6.49	600	15.93	5.64	0.030		5.41	604	16.31	5.79	0.142
N 300 − 1		10.72	724	15.76	5.77	0.022		10.99	800	15.80	5.74	0.027		9.17	922	16.03	5.80	0.066
N 2403 − 2		5.56	583	15.89	5.80	0.021		5.66	641	15.92	5.73	0.025		3.65	702	16.20	5.78	0.058
N 2403 − 3		21.93	469	16.18	5.90	0.316		19.50	501	16.24	5.91	0.370		16.98	651	16.31	5.85	0.445
N 2403 − 4		9.81	505	16.09	5.93	0.093		9.72	1114	15.58	5.82	0.009		7.95	1499	15.64	5.85	0.010
N 2403 − 5		5.30	502	16.10	6.03	0.053		5.35	600	16.06	5.93	0.044		3.41	628	16.39	6.00	0.129
M 81 − 5		28.38	616	15.86	5.63	0.094		28.71	705	15.85	5.59	0.090		24.38	913	15.94	5.59	0.116
M 81 − 6		37.84	655	15.79	5.55	0.090		35.40	722	15.81	5.53	0.093		29.17	957	15.87	5.52	0.101
M 81 − 11		5.26	586	15.93	5.90	0.024		5.38	650	15.96	5.83	0.028		3.36	702	16.25	5.88	0.067
M 81 − 12		7.16	410	16.26	6.00	0.149		7.33	500	16.21	5.93	0.121		5.62	517	16.55	6.05	0.444
M 81 − 14		26.85	400	16.34	5.95	0.808		25.47	496	16.23	5.85	0.462		21.68	601	16.36	5.85	0.715
N 7793 − 3		⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅		⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅		⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
N 7793 − 9		8.42	585	15.85	5.68	0.027		8.67	617	15.93	5.65	0.039		6.80	701	16.16	5.70	0.089
N 7793 − 10		6.24	472	16.17	6.05	0.086		6.56	512	16.23	6.03	0.119		4.56	600	16.44	6.05	0.217
N 7793 − 13		6.53	484	16.10	5.94	0.065		6.88	500	16.21	5.93	0.114		4.81	589	16.41	5.95	0.200
IRC+10420		10.99	600	15.73	5.44	0.020		11.51	700	15.77	5.48	0.025		10.00	1000	15.86	5.57	0.033
M 33 Var A		5.04	793	15.44	5.36	0.002		5.03	867	15.48	5.28	0.003		2.67	962	15.75	5.33	0.005
η Car		9.94	250	17.08	6.76	9.026		10.42	300	17.05	6.78	8.245		8.49	400	17.12	6.80	9.272

Download table as:  ASCIITypeset image

Figure 9 shows the integrated luminosities of the newly identified self-obscured stars described in Section 3.2 as a function of M_e for the best-fit graphitic models of each source. Object X, IRC+10420, M 33 Var A, and η Car are shown for comparison. Figure 10 shows the same quantities, but for various dust models and temperature assumptions. It is apparent from Figure 10 and Tables 3–5 that the integrated luminosity and ejecta mass estimates are robust to these uncertainties. The exceptions are N 2403-4 and N 7793-3. Without any optical or near-IR data, many of the models of N 7793-3 are unstable so we simply drop it. The only models having a luminosity in significant excess of 10⁶ L_☉ are some of the fixed temperature models of N 2403-4. These models have a poor goodness of fit and can be ignored.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One check on our selection methods is to examine the distribution of shell radii. Crudely, we can see a shell until it either becomes optically thin or too cold, so the probability distribution of a shell's radius assuming a constant expansion velocity is

$$\begin{equation} \frac{dN}{dR}_{\rm in} = \frac{1}{R_{\rm max}} = {\rm constant} \end{equation} \tag{ 3 }$$

for R_in < R_max. An ensemble of shells with similar R_max should then show this distribution. Figure 11 shows the cumulative histogram (excluding N 7793-3) of the inner shell radii (R_in). The curves show the expected distribution where we simply normalized to the point where F(< R_in) ≃ 0.5. The agreement shows that our sample should be relatively complete up to R_max ≃ 10^16.5–10¹⁷ cm which corresponds to a maximum age of

$$\begin{equation} t_{\rm max} \simeq 300\,{v_{e100}^{\,\,-1}}\; {\rm yr}. \end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

Figure 12 shows the age ($t_e=R_{\rm in}{{v_{e100}^{-1}}}$) of the shells as a function of M_e. We also show lines corresponding to optical depths of τ_V = 1, 10, 100. As expected, we see no sources with very low or high optical depths, as we should have trouble finding sources with τ_V < 1 due to a lack of mid-IR emission and τ_V ≳ 100 due to the dust photosphere being too cold (peak emission in the far-IR). Indeed, most of the dusty stars have 1 < τ_V < 10 and none has τ_V > 100. The large t_e estimate for η Car when scaled by v_e100 is due to its unusually large ejecta velocities (∼600 km s⁻¹ along the long axis; Cox et al. 1995; Smith 2006) compared to typical LBV shells (∼50–100 km s⁻¹; Tiffany et al. 2010).

Zoom In Zoom Out Reset image size

The advantage of surveying external galaxies with a significant SN rate is that we can translate our results into estimates of abundances and rates. We scale our rates using the observed SN rate of R_SN = 0.15 yr⁻¹ (0.05 < R_SN < 0.35 at 90% confidence). As we discussed in Paper I, this is significantly higher than standard star formation rate estimates for these galaxies, but the SN rate is directly proportional to the massive star formation rate rather than an indirect indicator, and similar discrepancies, although not as dramatic, have been noted in other contexts (e.g., Horiuchi et al. 2011). In this section we first outline how we will estimate rates, and then we discuss the constraints on analogs of η Car and the implications of our sample of luminous dusty stars.

We are comparing a sample of N_SN = 3 SNe observed over t_SN = 20 yr to a sample of N_c candidate stars that are detectable by our selection procedures for a time t_d. In Paper I we used DUSTY to model the detection of expanding dusty shells and found that a good estimate for the detection time period was

$$\begin{equation} t_d = t_w + 66 \left({ 100\, {\rm km\; s}^{-1} \over v_e } \right) \left({ L_* \over 10^6 L_\odot } \right)^{0.82} \left({ M_e \over M_\odot } \right)^{0.043} \;{\rm yr}\; \end{equation} \tag{ 5 }$$

for shells with masses in the range −1 ⩽ log M_e/M_☉ ⩽ 1 around stars of luminosity 5.5  ⩽  log L_*/L_☉  ⩽  6.5 where t_w is the duration of the "wind" phase and the second term is an estimate of how long the shell will be detected after the heavy mass loss phase ends. The principle uncertainty lies in the choice of the velocity, v_e. If the rate of events in the sample is R_e, then we expect to find N_e = R_et_d candidates.

The transient rate in a sample of galaxies is less interesting than comparing the rate to the SN rate. Let f_e be the fraction of massive (M_ZAMS  >  8 M_☉) stars that create the transients, where f_e  =  (M_C/8 M_☉)^−1.35 if we assume a Salpeter initial mass function (Kennicutt 1998), that all stars more massive than 8 M_☉ become SNe and that all stars more massive than M_C cause the transients. If each star undergoes an average of N_e eruptions, then the rate of transients is related to the rate of SNe by R_e = N_ef_eR_SN = F_eR_SN. The interesting quantity to constrain is F_e = N_ef_e rather than R_e. Poisson statistics provide constraints on the rates, where P(D|R)∝(Rt)^Nexp (− Rt) for N events observed over a time period t. This means that the probability of the rates given the data is

$$\begin{eqnarray} P(R_{\rm SN},R_e|D) &\propto& P(R_{\rm SN}) P(R_e) (R_{\rm SN} t_{\rm SN})^3 (R_e t_d)^{N_c}\nonumber\\ &&\times \exp (-R_{\rm SN}t_{\rm SN}-R_e t_d), \end{eqnarray} \tag{ 6 }$$

where P(R_SN) and P(R_e) are priors on the rates that we will assume to be uniform and we have set N_SN  =  3. If we now change variables to compute F_e and marginalize over the unknown SN rate, we find that the probability distribution for the ratio of the rates is

$$\begin{equation} P(F_e |D) \propto F_e^{N_c} \left(F_e t_d + t_{\rm SN}\right)^{-5-N_c}, \end{equation} \tag{ 7 }$$

with the standard normalization that ∫P(F_e|D)dF_e ≡ 1. For our estimates of F_e we present either 90% confidence upper limits or the value corresponding to the median probability and symmetric 90% probability confidence regions. Note that the probability distribution really just depends on the product F_et_d, so the results for any given estimate of t_d are easily rescaled.

4.1. No η Car Analog Is Found

4.2. An Emerging Class of Dust Obscured Stars

In our survey, we have found no true analogs of η Car. This implies that the rate of Great Eruption-like events is of order $F_e \,{=}\, 0.046\, t_{d200}^{-1}$ (0.0083  <  F_et_d200  <  0.19) of the ccSN rate, which is roughly consistent with each M ≳ 70 M_☉ star undergoing 1 or 2 such outbursts in its lifetime. This is scaled by an estimated detection period of order t_d = 200 t_d200 yr. We do identify a significant population of lower luminosity dusty stars that are likely similar to IRC+10420. Stars enter this phase at the rate $F_e = 0.20\, t_{d500}^{-1}$ (0.086 < F_et_d500 < 0.55) compared to the ccSN rate and for a detection period of t_d = 500 t_d500 yr. Here the detection period is assumed longer because the expansion velocities are likely slower. This rate is comparable to having all stars with 25 < M < 60 M_☉ undergoing such a phase once or twice.

If the estimated detection periods and mass ranges are roughly correct, and our completeness is relatively high, there are two interesting implications for both populations. First, these high optical depth phases represent a negligible fraction of the post-main-sequence lifetimes of these stars, at most lasting a few thousand years. This implies that these events have to be associated with special periods in the evolution of the stars. The number of such events a star experiences is also small, 1 or 2, not 10 or 20. Second, while a significant amount of mass is lost in the eruptions, they cannot be a dominant contribution to mass loss. For these high mass stars, standard models (e.g., Maeder 1981; Maeder & Meynet 1987, 1988; Stothers & Chin 1996; Meynet et al. 1994) typically strip the stars of their hydrogen envelopes and beyond, implying total mass losses of all but the last 5–10 M_☉. The median mass loss in Figure 9 is M_e  ∼  0.5 M_☉ and if every star underwent two eruptions, the typical total would be N_eM_e ∼ M_☉. Clearly there are some examples that require significantly larger M_e, but we simply do not find enough heavily obscured stars for this phase to represent more than a modest fraction of the total mass loss (∼10% not ∼50%).

For the stars similar to IRC+10420, this is consistent with the picture that the photospheres of blue-ward evolving RSGs with log (L_*/L_☉)  =  5.6  ∼  6.0 and T_star  ≃  7000–12, 500 K, become moderately unstable, leading to periods of lower effective temperature and enhanced mass loss as the stars try to evolve into a "prohibited" region of the HR diagram that de Jager & Nieuwenhuijzen (1997), de Jager (1998), and Nieuwenhuijzen & de Jager (2000) termed the "yellow void." In this phase, the stars lose enough mass to evolve into a hotter, less massive star on the blue side of the HR diagram. This is also the luminosity regime of the "bistability jump" in wind speeds driven by opacity changes, which Smith et al. (2004) hypothesizes can explain the absence of LBVs and the existence of YHGs with high mass loss rates and dust formation (Vink et al. 2012) in this luminosity range. In fact, Humphreys et al. (2002) propose that IRC+10420, which is identified by our selection criterion, is such a star. While these arguments supply a unique, short-lived evolutionary phase, there may be problems with the absolute scale of the mass loss, since estimates are that IRC+10420 started with a mass of ∼40 M_☉ and has lost all but 6 ∼ 15 M_☉ (Nieuwenhuijzen & de Jager 2000).

The only other similarly unique phase in the lives of these stars is the final post-carbon ignition phase. There are now many examples of stars that have had outbursts shortly before exploding as SNe (e.g., Pastorello et al. 2007, 2013; Mauerhan et al. 2013; Prieto et al. 2013; Ofek et al. 2013) and superluminous SNe that are most easily explained by surrounding the star with a large amount of previously ejected mass (Smith & McCray 2007; Gal-Yam et al. 2007; Smith et al. 2008; Kozlowski et al. 2010; Ofek et al. 2013). Powering these SNe requires mass ejected in the last years to decades of the stellar life (e.g., Chevalier & Fransson 1994; Chugai & Danziger 2003; Smith 2009; Moriya et al. 2014) and it seems natural to associate these events with the mass ejections of LBVs like η Car (e.g., Smith & McCray 2007; Gal-Yam & Leonard 2009). The statistical properties and masses of either of the classes of dusty stars we discuss are well-matched to the statistical requirements for explaining these interaction powered SNe if the instability is associated with the onset of carbon burning (see Kochanek 2011). If there is only one eruption mechanism, it must be associated with a relatively long period like carbon burning (thousands of years) rather than the shorter, later nuclear burning phases, because we observe many systems like η Car that have survived far longer than these final phases last. If the mechanism for producing the ejecta around the superluminous SNe is associated with nuclear burning phases beyond carbon, then we must have second eruption mechanism to explain η Car or other still older LBVs surrounded by massive dusty shells. If there indeed are two mass loss mechanisms—one commencing ≳ 10³ yr from core-collapse and the other occurring in the ∼1 yr prior to core-collapse—then the self-obscured stars identified in this work may very well be experiencing the earlier of these two mechanisms. Otherwise, in a larger sample of ∼100 such stars, one should be exploding as a ccSN every ∼10 yr.

The dusty stars can be further characterized by their variability, which will help to follow the evolution of the dust. For the optically brighter examples, it may be possible to spectroscopically determine the stellar temperatures, although detailed study may only become possible with the James Webb Space Telescope (JWST). It is relatively easy to expand our survey to additional galaxies. For very luminous sources like η Car analogs, this is probably feasible to distance of 10 Mpc, while for the lower luminosity IRC+10420 analogs this is likely only feasible at the distances of the most distant galaxies in our sample (∼4 Mpc). Larger galaxy samples are needed not only to increase the sample of dusty luminous stars (and hopefully find a true η Car analog!), but also to have a sample with a larger number of SNe, or equivalently a higher star formation rate. Our estimate of the abundance of IRC+10420 analogs is limited by the small number of ccSN (3) in our sample more than by the number of dusty stars (18) identified. Finally, while we have shown that surveys for the stars are feasible using archival Spitzer data, the JWST will be a far more powerful probe of these stars. The HST-like resolution (Gardner et al. 2006) of the JWST will be enormously useful to either greatly reduce the problem of confusion or to expand the survey volume. Far more important will be the ability to carry out the survey at 24 μm, which will increase the time over which dusty shells can be identified from hundreds of years to thousands of years, greatly improving the statistics and our ability to survey the long term evolution of these systems and the relationship between stellar eruptions and SNe.

We thank the referee for insightful comments. Hendrik Linz for helping us analyze the Herschel PACS data and John Beacom for numerous productive discussions. We extend our gratitude to the SINGS Legacy Survey and LVL Survey for making their data publicly available. This research has made use of observations made with the NASA/ESA Hubble Space Telescope, and obtained from the Hubble Legacy Archive, which is a collaboration between the Space Telescope Science Institute (STScI/NASA), the Space Telescope European Coordinating Facility (ST-ECF/ESA) and the Canadian Astronomy Data Center (CADC/NRC/CSA). R.K. and K.Z.S. are supported in part by NSF grant AST-1108687.