Constraining the Rate of Protostellar Accretion Outbursts in the Orion Molecular Clouds

TL;DR

This study searches for protostellar outbursts in the Orion molecular clouds using infrared data over 6.5 years, finding only two bursts and estimating an average interval of about 1000 years between such events.

Contribution

It provides the first statistical estimate of outburst frequency in protostars using infrared variability over several years.

Findings

Detected two protostellar outbursts, including a new one.

Estimated an average interval of ~1000 years between outbursts.

Suggests outbursts are more frequent in protostars than in older stars.

Abstract

Outbursts due to dramatic increases in the mass accretion rate are the most extreme type of variability in young stellar objects. We searched for outbursts among 319 protostars in the Orion molecular clouds by comparing 3.6, 4.5, and 24 micron photometry from the Spitzer Space Telescope to 3.4, 4.6, and 22 micron photometry from the Wide-field Infrared Survey Explorer (WISE) obtained ~ 6.5 yr apart. Sources that brightened by more than two standard deviations above the mean variability at all three wavelengths were marked as burst candidates, and they were inspected visually to check for false positives due primarily to the reduced angular resolution of WISE compared to Spitzer. We recovered the known burst V2775 Ori (HOPS 223) as well as a previously unknown burst, HOPS 383, which we announced in an earlier paper. No other outbursts were found. With observations over 6.5 yr, we estimate an interval of about 1000 yr between bursts with a 90% confidence interval of 690 to 40,300 yr. The most likely burst interval is shorter than those found in studies of optically revealed young stellar objects, suggesting that outbursts are more frequent in protostars than in pre-main-sequence stars that lack substantial envelopes.

Figures

Tables

Table 1. Table 1: Two-Epoch MIPS 24 µm Photometry

	2004–2005aaUsed in this paper.		2008bbUsed in Furlan et al. (2016).
HOPS ID	Mag.	Unc.ccAll uncertainties include a 5% floor.	Mag.	Unc.ccAll uncertainties include a 5% floor.
86	1.788	0.055	2.105	0.057
87	2.490	0.055	2.396	0.061
88	3.268	0.055	2.840	0.066
89	3.779	0.055	3.645	0.064
90	1.842	0.055	1.911	0.072
91	4.095	0.056	3.821	0.069
92	0.451	0.055	0.571	0.057
93	4.311	0.055	4.476	0.066
94	1.332	0.055	1.391	0.077
95	5.167	0.056	5.174	0.062
96	4.252	0.055	4.222	0.056
99	5.256	0.057	5.459	0.077
100	5.977	0.057	5.779	0.062
102	3.700	0.056	3.706	0.058
105	5.349	0.057	5.442	0.066
107	1.517	0.056	1.526	0.058
204	3.388	0.061	3.507	0.055
206	4.087	0.055	3.746	0.055
207	3.810	0.055	3.582	0.055
208	6.787	0.059	7.082	0.061
209	4.473	0.055	4.481	0.055
210	2.854	0.055	2.901	0.055
211	5.996	0.057	6.071	0.057
213	2.622	0.056	2.729	0.057
214	5.030	0.055	5.180	0.056
215	3.545	0.055	3.563	0.055
216	3.147	0.055	3.231	0.056
219	2.546	0.056	2.505	0.056
220	4.755	0.057	4.463	0.056
221	0.981	0.055	1.022	0.055
223	2.548	0.055	0.389	0.056
224	4.580	0.056	4.423	0.056
225	3.597	0.059	3.444	0.055
226	3.752	0.059	3.514	0.055
227	4.378	0.055	4.445	0.055
228	0.410	0.055	0.714	0.055
232	2.536	0.056	2.290	0.061
235	1.429	0.056	1.580	0.067
383	7.875	0.057	4.015	0.057

Table 2. Table 2: HOPS Protostars with WISE  Counterparts

	3.6 µm	4.5 µm	24 µm	All
Spitzer and WISE	248	260	249	233
Spitzer only	3	0	20	$\cdots$
WISE only	16	10	0	$\cdots$
No detection	4	1	2	$\cdots$
Mean difference (mag)aaMean and standard deviation of Spitzer magnitude minus WISE magnitude for sources detected in both catalogs.	$-0.283$	$0.356$	$0.013$
Std. dev. of difference (mag)aaMean and standard deviation of Spitzer magnitude minus WISE magnitude for sources detected in both catalogs.	$0.646$	$0.590$	$0.620$

Table 3. Table 3: Outburst Candidates

HOPS	RA	Dec	Magnitude change
ID	(deg)	(deg)	3.6 µm	4.5 µm	24 µm
223	85.7019	$-$8.2762	1.69	2.23	1.77
383	83.8742	$-$4.9975	$>$2.9	2.40	3.96

Equations

$$r_b=\frac{n_b}{N_p\Delta t},$$

$$\begin{gathered} P\left(n_b|I_b,N_p,\Delta t\right)= \\ \frac{N_p!}{n_b!(N_p-n_b)!}(1-e^{-\Delta t/I_b}{)}^{n_b}(e^{-\Delta t/I_b}{)}^{(N_p-n_b)}. \end{gathered}$$

Full text

circumstellar matter — infrared: stars — stars: formation — stars: protostars

1 INTRODUCTION

Young stellar objects (YSOs) display a wide range of photometric variability. Common sources of variability include the rotation of the star, which can carry hot or cool spots across the line of sight (Herbst et al., 1994), the rotation of the protoplanetary disk, which can carry orbiting disk structures across the line of sight (Bouvier et al., 2007; Morales-Calderón et al., 2011), changes in the structure of the disk (Rice et al., 2012; Wolk et al., 2013; Flaherty et al., 2012, 2013), the influence of a companion star (Muzerolle et al., 2013), or changes in the foreground extinction (Chiang et al., 2015). The most extreme examples of variability, however, are due to changes in the luminosity of the system. While the luminosity of a low-mass pre-main-sequence star is expected to change gradually over $\sim 10^{6}$ yr as it contracts toward the main sequence, the luminosity generated by disk material accreting onto the star may vary on much shorter timescales due to changes in the accretion rate.

Episodic accretion is a phenomenon in which the accretion rate of a forming star rapidly increases by as much as several orders of magnitude (Hartmann & Kenyon, 1996). In the case of protostars, where accretion can dominate the total luminosity, it is one potential solution to the classic luminosity problem, where luminosities predicted for protostars that accrete at a constant rate throughout the formation period are greater than those typically observed (Kenyon et al., 1990; Evans et al., 2009). With episodic accretion, the protostar is usually less luminous than such predictions indicate, in agreement with observations, but it spends enough time at an elevated luminosity for its luminosity averaged over the entire period to agree with predictions. Many examples of episodic accretion have been observed since the first discovery of an outburst, that of FU Ori, in 1936 (Wachman, 1939). While outbursts clearly occur, the fraction of a typical star’s main-sequence mass that is accumulated during bursts is still poorly constrained.

Investigating the overall importance of bursts to star formation, Dunham et al. (2010) showed that models featuring episodic accretion were successful in reproducing the bolometric luminosity and temperature (BLT) distribution, a common evolutionary diagram, of 1024 nearby YSOs. Offner & McKee (2011), on the other hand, found that models with a limited role for episodic accretion but with constant accretion times reproduced observations better than those with constant accretion rates. To explain the BLT distribution of 315 Orion protostars, Fischer et al. (2017) invoked exponentially declining accretion rates over the star formation period with, again, a limited role for episodic accretion.

The reviews of Reipurth & Aspin (2010) and Audard et al. (2014) tabulate the known outbursts, discuss their observational characteristics, and explore potential triggering mechanisms. The two main types of bursts are EX Lupi and FU Orionis events. EX Lupi bursts are relatively small in amplitude and short-lived, but they are also known to recur in the same star (Aspin, 2010).

Conversely, the amplitudes of FU Orionis bursts are greater. Their decay timescales are far longer, more than several decades, but with significant variation in their post-burst light curves. Bursts of this type may recur in the same star, but the timescales are too long for this to have yet been observed. Several outbursts do not fit cleanly into either class, such as V2492 Cyg (Covey et al., 2011) or V1647 Ori (Muzerolle et al., 2005; Fedele et al., 2007; Contreras Peña et al., 2017a), indicating that the usual division of bursts into two classes is a simplification.

The predominant explanation for the jump in accretion rate invokes disk instabilities, either thermal or gravitational, triggered by growth in the disk mass due to accumulation of infalling envelope material. In thermal instability models, the gas accumulates in the disk as it is blocked from accreting onto the star by an orbiting planet or a dead zone, eventually leading to a rapid rise in temperature and accretion rate as the material falls onto the star (Lodato & Clarke, 2004; Armitage et al., 2001; Zhu et al., 2010). In the gravitational instability models, the growing disk becomes unstable and fragments into clumps, which migrate into the inner disk. The system begins a luminosity outburst when the clumps accrete onto the star (Vorobyov & Basu, 2005, 2015).

In this work we use the common classification system for YSOs. Classes I, II, and III were defined by Lada (1987), and the system was extended with the Class 0 (André et al., 1993) and flat-spectrum (Greene et al., 1994) categories. For this study, we consider Class 0, Class I, and flat-spectrum YSOs to be protostars, because their circumstellar envelopes make an important contribution to their observed properties and subsequent evolution. Furlan et al. (2016) provide detailed criteria for distinguishing among these classes and compare the properties of the different classes of protostars in Orion.

Though more widely studied within the last two decades, many characteristics of accretion bursts are still not well constrained. This is due in part to the serendipitous nature of the discovery of most outbursts, which limits the effectiveness of statistical techniques in estimating their frequency. Further, detections of protostellar outbursts are rare compared to detections of other types of YSO variability. Observational studies of large samples of YSOs are essential in constraining the characteristics of the outburst process and in guiding theoretical work.

The advent of wide-area, IR surveys of large numbers of diverse star-forming regions makes systematic searches for accretion bursts feasible. Scholz et al. (2013), for example, used two epochs of photometry for about 8000 YSOs covering a wide range of masses and ages, searching for luminosity increases that satisfy specific outburst criteria, with the ultimate goal of constraining the frequency of outbursts in their sample. They found a burst interval of 5000 to 50,000 yr. For all sample sizes and time baselines studied to date, however, the upper bound of the burst interval is poorly constrained (Hillenbrand & Findeisen, 2015).

Here we present the results from an analysis of 319 well characterized protostars in the Orion molecular cloud complex, using two epochs of mid-IR photometry from the Spitzer Space Telescope (Werner et al., 2004) and the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010). This survey is the largest to date for protostellar outbursts, given the large population of young YSOs in the Orion region. Like Scholz et al. (2013), we will use the results of this survey to estimate the frequency of outbursts between the two epochs.

2 DATA AND SAMPLE SELECTION

As part of the Spitzer Orion Survey (Megeath et al., 2012, 2016), fields in the Orion A and B molecular clouds were mapped in 2004 and 2005 with the Infrared Array Camera (IRAC; Fazio et al. 2004) and the Multiband Imaging Photometer for Spitzer (MIPS; Rieke et al. 2004). Images and photometry were obtained at five wavelengths, including bands centered at 3.6 (IRAC 1), 4.5 (IRAC 2), and 24 µm (MIPS 1). The FWHM of the point-spread function at each wavelength is approximately 2″ for the IRAC channels and 6″ for MIPS. Megeath et al. (2012) show that the $5\sigma$ detection limits in each band are typically 17 mag in IRAC 1 and IRAC 2 and 9 mag in MIPS 1; these limits vary across the Orion survey due to confusion with spatially varying nebulosity. They also show the coverage of the maps and the dates on which they were obtained.

The Spitzer data were used to identify targets for the Herschel Orion Protostar Survey (HOPS; Manoj et al. 2013; Stutz et al. 2013; Furlan et al. 2016; Fischer et al. 2017), an open-time key program of the Herschel Space Observatory to obtain 70 and 160 µm imaging and 50–200 µm spectra of Orion protostars. Furlan et al. (2016) presented 1 to 870 µm SEDs and model fits for 330 HOPS targets in Orion. These are protostellar candidates that were imaged by HOPS and detected in the Herschel 70 µm images. Of the 330 YSOs, Furlan et al. (2016) concluded from SED analysis that 319 are protostars and the remaining 11 are Class II YSOs. These 319 protostars, those that appear in Table 1 of Furlan et al. (2016) and have Class designations of 0, I, or flat, constitute the sample for this work. Their coordinates and SED properties are tabulated there as well.

Our Spitzer photometry is from the 2004–2005 campaigns discussed above. The data are published in Table 2 of Furlan et al. (2016), except that those authors used more recent MIPS photometry for 39 of the protostars and more recent IRAC photometry for HOPS 383. In Table 1, we list the original and previously unpublished 2004–2005 MIPS photometry used here for these sources as well as their previously published 2008 MIPS photometry. Both epochs of IRAC photometry for HOPS 383 appear in Safron et al. (2015).

The WISE mission surveyed the entire sky in four bands, including 3.4 (W1), 4.6 (W2), and 22 µm (W4). WISE mapped Orion at all wavelengths in 2010 March, and after the exhaustion of cryogen, the region was mapped again at W1 and W2 in 2010 September. This work uses the AllWISE version of the WISE point-source catalog, which incorporates both epochs. The FWHM of the point-spread function for W1 and W2 is approximately 6″, while that for W4 is approximately 12″. The quoted $5\sigma$ sensitivities for AllWISE are 16.9 mag in W1, 15.9 mag in W2, and 8.0 mag in W4, although these sensitivities can be considerably worse in a high-background region such as Orion.444See the discussion at http://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec2_3a.html. The three WISE bands of interest have response functions similar to Spitzer’s IRAC 1, IRAC 2, and MIPS 1 bands, respectively. We find at most a weak correlation between the intrinsic Spitzer color of a source and the difference between its Spitzer and WISE magnitudes, so we conclude that differences in the Spitzer and WISE response functions are a minor effect compared to true outbursts, and we directly compare the two sets of magnitudes. For simplicity we use 3.6, 4.5, and 24 µm to refer to the central wavelengths for both telescopes.

Comparison of the two data sets allows the discovery of variability on a timescale of 5 to 6 yr. Due to the larger point-spread functions of WISE, blending with nearby sources is a concern, and visual inspection is needed to conclusively identify variable sources.

Our process for matching WISE point sources to the HOPS protostars is that described in Fischer et al. (2016a). We used the NASA/IPAC Infrared Science Archive555See http://irsa.ipac.caltech.edu. to search AllWISE for the closest match to each Spitzer position. For 227 of 319 protostars, there was a match within 1″; these were automatically paired with their Spitzer counterparts. For 32 protostars, there was no WISE source within 10″; these were automatically marked as having no counterpart. For the 60 sources at intermediate separations, we manually inspected the images and found that the larger offset was due either to scattered light from the Spitzer point source or to a blend in WISE of two or more distinct Spitzer sources. We retained the 44 matches where the WISE flux was judged to be mainly due to the Spitzer source in question (all with separations less than 7″) and rejected the other 16. Therefore, 271 of the 319 protostars have WISE counterparts.

A comparison of the detections by Spitzer and WISE is shown in Table 2. Of the 271 protostars with WISE counterparts, 233 are detected by both Spitzer and WISE in all three bands of interest. Considering the bands independently, in each band 91% to 96% of the 271 protostars are detected in both catalogs. Of those that were detected in both, the mean difference in magnitude is slightly negative at 3.6 µm, slightly positive at 4.5 µm, and nearly zero at 24 µm. The standard deviations in these differences are about 0.6 mag in all bands. Throughout this work, we subtract WISE magnitudes from Spitzer magnitudes so that positive numbers correspond to an increase in brightness from one epoch to the next.

To identify outburst candidates, we use these typical differences as a guide and flag protostars that became substantially brighter between the Spitzer and WISE maps or were detected in the later WISE map after not being detected by Spitzer. We then manually inspect these candidates to assess whether the different angular resolutions of the two telescopes play a role in the reported changes in brightness. In the following section we present and justify outburst criteria.

3 CANDIDATE CRITERIA

Depending on the underlying physics, YSO variability is characterized by a range of amplitudes and color changes. Here we are interested in major increases in the accretion rate from the disk onto the star. In this section we discuss the expected effect of such increases on the SED at 3.6, 4.5, and 24 µm and develop criteria for identifying outburst candidates. We argue that a large fractional change in the accretion rate will induce fractional changes of a similar magnitude in the flux densities at these wavelengths.

The accretion rate from the disk onto the star $\dot{M}$ is proportional to the accretion luminosity $L_{\rm acc}$, where $L_{\rm acc}=\eta\,GM_{*}\dot{M}/R_{*}$. Here $G$ is the gravitational constant, $M_{*}$ and $R_{*}$ are the mass and radius of the central object, and $\eta$ is a constant of order unity that depends on the details of the accretion process. A large increase in $\dot{M}$ will drive a proportional change in the accretion luminosity, since the stellar parameters do not change on timescales short enough to be relevant and changes in $\eta$ are expected to be small ($\lesssim 20\%$) compared to the change in $\dot{M}$, which is typically one order of magnitude or more (Baraffe et al., 2012).

The total luminosity of the system $L$ is the sum of the accretion luminosity and the luminosity of the central star $L_{*}$. Like other parameters intrinsic to the star, $L_{*}$ is not expected to change appreciably due to the outburst. The ratio of the total luminosity in outburst to that in the quiescent state is then $(L_{{\rm acc},b}+L_{*})/(L_{{\rm acc},q}+L_{*})$. This is a lower limit to the ratio of accretion luminosities, $L_{{\rm acc},b}/L_{{\rm acc},q}$, and it approaches the latter ratio when both accretion luminosities are much larger than that of the star, the expectation early in protostellar evolution (Adams et al., 1987; André et al., 1993; Fischer et al., 2017). Therefore, a large fractional increase in the total luminosity of the system implies even larger fractional increases in the accretion luminosity and accretion rate.

In the SED, a change in the total luminosity may affect the flux density differently at different wavelengths. Figure 11 of Furlan et al. (2016) shows how, as the luminosity of a protostar increases from 0.1 to 303 $L_{\sun}$ in steps of about half an order of magnitude, the fractional change in flux density predicted by a radiative transfer model at 3.6, 4.5, and 24 µm is similar to the fractional change in luminosity. This was also indicated by Kenyon et al. (1993), who calculated that the peak wavelength of a protostellar SED varies as $L^{-1/12}$; i.e., a change in luminosity primarily affects the overall flux in the SED, not its shape. Finally, Scholz et al. (2013) analyzed model SEDs for prototypical Class I and II sources at various accretion rates. They found that, between 2 and 5 µm, the flux increases by a factor of 10 or more when the accretion rate increases from zero to $10^{-6}~{}M_{\sun}~{}{\rm yr}^{-1}$ or more. These lines of reasoning suggest that dramatic increases in luminosity can be detected by searching for simultaneous large increases in flux density at all three wavelengths of interest.

The typical short-term near- or mid-IR variations in large samples of YSOs, which may be due to rotational modulation by spots on the star or inner disk inhomogeneities, are in the range of 0.1–0.6 mag (Morales-Calderón et al., 2011; Flaherty et al., 2012; Megeath et al., 2012) and can be strongly wavelength dependent. The slightly different Spitzer and WISE bandpasses can mimic variability between the two surveys of up to a few tenths of a magnitude. Low-level variability and differences in the bandpasses are likely responsible for the distributions in magnitude differences presented in Table 2, where the mean and standard deviation in each band are of order a few tenths of a magnitude.

To filter out low-level variability and differences due to the different bandpasses, we search for protostars with brightenings greater than two standard deviations above the mean at all three wavelengths of comparison, or 1.01 mag at 3.6 µm, 1.54 mag at 4.5 µm, and 1.25 mag at 24 µm. We also consider protostars that exceed the threshold at one or two wavelengths while changing from a non-detection to a detection at the other(s). These changes correspond to an increase in luminosity of a factor of approximately 2.5–4. Since the variations in the photometry are dominated by low-level variability as well as by potential source confusion in the WISE data, the two-sigma threshold should not be considered the statistical significance of the bursts. Instead, it gives criteria for distinguishing large bursts from smaller-scale variability and systematic effects in the data. The thresholds are expected to pick up moderate outbursts as well as those in which neither epoch of photometry exactly catches the minimum or maximum brightness.

4 ANALYSIS OF BRIGHTNESS CHANGES

In this section we plot the difference in magnitude of each protostar at the three wavelengths of comparison, identify those that brightened by more than the threshold at one or more wavelengths, and examine false positives.

Figure 1 plots the IRAC 1 $-$ W1 magnitude against the IRAC 1 magnitude for the 248 protostars detected at both epochs. Nine protostars lie above the threshold for outbursts; these points are labeled with their HOPS numbers. Five of them satisfy the requirement only at this wavelength; HOPS 28, 335, and 370 satisfy the requirement at this and one other wavelength; and HOPS 223 satisfies the requirement at all three wavelengths.

Figure 2 plots the IRAC 2 $-$ W2 magnitude against the IRAC 2 magnitude for the 260 protostars detected in both bands. Eight protostars lie above the threshold for outbursts. Five of them satisfy the requirement only at this wavelength, and two of them satisfy the requirement at this and one other wavelength. Of these two, HOPS 370 brightened at 3.6 and 4.5 µm but not 24 µm. The other, HOPS 383, was not detected at 3.6 µm by Spitzer but brightened at the other two wavelengths. Again, HOPS 223 satisfies the requirement at all three wavelengths.

Figure 3 plots the MIPS 1 $-$ W4 magnitude against the MIPS 1 magnitude for the 249 protostars detected in both bands. Eight protostars lie above the threshold for outbursts. Four of them brightened only at this wavelength, and three of them brightened at one other wavelength. These are HOPS 28, 335, and 383, which are all mentioned above. Again, HOPS 223 satisfies the requirement at all three wavelengths.

Differences at one wavelength are plotted against those at a second wavelength in the next three figures. Figure 4 plots the difference in magnitude at 3.6 µm against the difference at 4.5 µm for the 248 protostars detected by Spitzer and WISE at both wavelengths. These differences are moderately well correlated, with a correlation coefficient of 0.60. The wavelengths are close enough that sources are unlikely to brighten or fade substantially at one but not the other. The most flagrant exception, with a putative brightening of more than 2 mag at 3.6 µm but a fading at 4.5 µm, is HOPS 40, above and slightly to the left of the main cloud of points. It is faint in IRAC 1 (14.7 mag) and contaminated by background emission in W1.

Figure 5 plots the difference in magnitude at 3.6 µm against the difference at 24 µm for the 233 protostars detected by Spitzer and WISE at both wavelengths, and Figure 6 plots the difference in magnitude at 4.5 µm against the difference at 24 µm for the 241 protostars detected by both telescopes at both wavelengths. These differences are not well correlated, with correlation coefficients of 0.21 and 0.37. The Spitzer 24 µm and WISE 22 µm bands are far enough in wavelength space from the 3–5 µm bands, and the WISE 22 µm band is sufficiently affected by confusion due to its lower angular resolution, that changes in the former do not necessarily track those in the latter.

4.1 False Positives

Here we look at the three protostars that were detected at all three wavelengths by both telescopes and satisfy the outburst criteria at just two of the three wavelengths. These are the protostars that lie above and to the right of the blue lines, but are not circled, in the figures that plot one magnitude difference against another. Figure 4 shows that HOPS 370 is part of this collection, and Figure 5 highlights HOPS 28 and 335. With visual inspection, their putative outbursts can all be rejected as false positives due to contamination of the relatively large WISE point-spread function by nebulosity or nearby point sources. The remaining sources that satisfy only one threshold are also all false positives. They will not be examined in detail here, but we found that they suffer from the same issues described below.

Figure 7 shows thumbnails of HOPS 370, which satisfies the criteria only at 3.6 and 4.5 µm. At those wavelengths, it became brighter by 3.9 and 2.7 mag, respectively, while at 24 µm, it apparently faded by 0.1 mag. HOPS 370 is a binary protostar (Nielbock et al., 2003) that is barely resolved by IRAC. The saturation evident in the MIPS image does not affect the analysis, since the photometry was obtained by fitting a point-spread function with the saturated pixels masked (Kryukova et al., 2012). In W1 the components are blended, while in W2, the sources are visibly separated. In both, extended emission and the bright source just to the northwest (HOPS 66) contribute to the photometry, while in W4, the source is so bright that contamination by these factors is minimal. The erroneously elevated photometry in W1 and W2 make this candidate a false positive.

Figure 8 and Figure 9 show thumbnails of HOPS 28 and HOPS 335, the two protostars that satisfy the criteria only at 3.6 and 24 µm. HOPS 28 apparently brightened by 1.1 mag at 3.6 µm, 1.2 mag at 4.5 µm, and 1.3 mag at 24 µm. These represent increases of 2.1, 1.4, and 2.1 standard deviations; i.e., it slightly exceeded the criteria at two bands and fell short at the other. HOPS 335 apparently brightened by 1.3 mag at 3.6 µm, 1.4 mag at 4.5 µm, and 1.9 mag at 24 µm, for increases of 2.4, 1.7, and 3.0 standard deviations, coming close to meeting all criteria. In both of these cases, extended emission appears to be contaminating the photometry in W1 and W2, while a nearby, much brighter source contaminates the photometry in W4.

In summary, the larger WISE point spread functions tend to attribute environmental nebulosity to point sources. Blending with nearby point sources is also a concern, particularly at 22 µm. Visual inspection is a crucial step to identify false outbursts in the comparison of Spitzer and WISE photometry. Alternatively, blending of sources in a particularly unlucky way may result in the non-detection of a real burst.

5 OUTBURSTS

The only protostars that fulfill all criteria, i.e., HOPS 223 and HOPS 383, have already been discussed in the literature in references such as those given below. Visual inspection confirms that their brightness changes from Spitzer to WISE are reflective of their outbursts, not source confusion in the large WISE beam. These changes are summarized in Table 3. HOPS 223 brightened by more than the threshold in all three bands and is marked with a purple circle in Figures 1 through 6. HOPS 383 changed from a non-detection to 14.1 mag at 3.6 µm, which is consistent with a brightening of at least 2.9 mag given the IRAC sensitivity limit of 17 mag. The minimum change at the shortest wavelength was larger than the 1.0 mag threshold there, and the source also easily satisfied the criteria at the other two wavelengths. It is marked with an orange circle in the relevant figures. Although HOPS 223 and HOPS 383 are both in regions that were mapped by MIPS in 2004 and 2008, the full comparison of Spitzer to WISE data indicates that these are the only two sources that erupted during the 6.5 yr interval. Their fortuitous location at the overlap of the two MIPS fields serves only to confirm the discovery of their bursts.

5.1 HOPS 223

Our comparison of Spitzer and WISE data was motivated by the large population of protostars in Orion and the announcement by Caratti o Garatti et al. (2011) that HOPS 223, a flat-spectrum protostar in the Lynds 1641 region, had undergone an outburst. The protostar was identified as [CTF93]216-2 in the discovery paper due to its proximity to the protostar [CTF93]216 (Chen et al. 1993; HOPS 221 in our survey). It was subsequently given the variable-star identifier V2775 Ori (Kazarovets et al., 2011). Here we refer to it by its HOPS number for consistency with the rest of our sample.

Our technique recovers the HOPS 223 outburst. The magnitude changes at 3.6, 4.5, and 24 µm were 1.7, 2.2, and 1.8 mag, respectively, corresponding to increases in flux density by factors of 4.7, 7.9, and 5.3. These correspond to magnitude changes of 2.8 to 3.2 standard deviations above the mean. The source offers solid detections in every band at both epochs, and the outburst is obvious in the corresponding images, shown in Figure 10. Before the outburst, its 3.6 and 4.5 µm flux densities were similar to those of HOPS 221 to the southwest, and at 24 µm it was fainter than HOPS 221. After the outburst, it was the brightest source in the field at all three wavelengths.

Any blending of HOPS 223 with the star 11″ to its northwest is insignificant in the post-outburst WISE images. This star was found by Megeath et al. (2012) not to have an IR excess, and it is undetected at 24 µm. In the pre-outburst Spitzer images, it is increasingly fainter than HOPS 223 with increasing wavelength, going from 2.6 mag fainter at 3.6 µm to $>6.5$ mag fainter at 24 µm if the detection limit there is 9 mag. The contrast would be even greater post-outburst.

Based on photometry and imaging spanning a range of regimes from submillimeter to near-IR, HOPS 223 was determined to have risen in luminosity from 4.5 to $\sim 51$ $L_{\sun}$, a factor of about 12, sometime between 2005 April 2 and 2007 March 12. Fischer et al. (2012) used criteria developed by Connelley & Greene (2010) to determine that the near-IR spectrum of HOPS 223 after its outburst was consistent with that of an FU Orionis object. Though not as extreme as classical examples of FU Orionis outbursts such as FU Ori itself or V1057 Cyg, the accretion rate of HOPS 223 increased from about $10^{-6}$ to $10^{-5}~{}M_{\sun}~{}{\rm yr}^{-1}$ during the outburst.

Fischer et al. (2016b) reported additional observations of HOPS 223. In 2015 November they obtained mid-IR photometry with the Stratospheric Observatory for Infrared Astronomy (SOFIA), and in 2016 January they obtained near-IR spectroscopy at Palomar Observatory. These observations, to be discussed in detail in future work, revealed that HOPS 223 had faded slightly but retained the spectral characteristics of an FU Orionis star roughly a decade after its initial outburst.

5.2 HOPS 383

HOPS 383 is a non-detection in the IRAC 1 catalog, but it is firmly detected in all three WISE bands analyzed here. Its magnitude changes at 4.5 and 24 µm were 3.5 and 6.4 standard deviations above the mean, and its magnitude change at 3.6 µm was at least 4.9 standard deviations above the mean. Visual inspection of this source revealed a clear brightening in all three wavelength bands. Images appear in Figure 11.

We previously announced the brightening of this protostar in Safron et al. (2015). We estimated that its outburst began between 2004 October 12 and 2006 October 20. Analysis of the HOPS 383 SED, including its detection at submillimeter wavelengths, unambiguously indicates that it is a Class 0 YSO, with 1.5% of its total luminosity emitted at 350 µm and longer, a bolometric temperature of 43 K, and a peak wavelength of 106 µm. Its large submillimeter fluxes imply that it is embedded in a massive envelope. Although the coverage of the pre-outburst SED is sparse, the bolometric luminosity was estimated to increase from 0.2 $L_{\sun}$ to 7.5 $L_{\sun}$, consistent with a legitimate accretion burst. HOPS 383 appears to be the least evolved source known to undergo an outburst.

Subsequent study of HOPS 383 has revealed a range of behaviors at different wavelengths. Galván-Madrid et al. (2015) analyzed archival observations of HOPS 383 at 3 cm. They found two components that are separated by about 0.45″ and connected by a ridge of faint emission, likely due to a bipolar jet or a binary. They found only mild variations of the radio flux between 1998 and 2014, interpreting the lack of a radio outburst as evidence that an ejection enhancement did not follow the accretion outburst. Mairs et al. (2017) compared 850 µm data obtained between 2012 August and 2013 August to data obtained at the same wavelength between 2015 December and 2017 Feburary. They found that HOPS 383 was a “possible” variable candidate, undergoing a brightness decrease of $2.66\pm 0.64\%$ yr*-1*. Imaging of the field at $J$, $H$, and $K_{s}$ in 2015 December resulted in non-detections (Fischer & Hillenbrand, 2017), suggesting a decline in the near-IR flux density after the acquisition of post-outburst data between 2008 and 2011. Further study of this deeply embedded outburst is essential for understanding episodic accretion in the youngest protostars.

5.3 Other Outbursts

Additional protostellar bursts were known in the Orion molecular clouds before the Spitzer observations that begin our baseline. They are not considered in our statistical treatment, but we review them here for completeness. Just before the beginning of Spitzer’s scientific operations in late 2003, V1647 Ori (HOPS 388) began an outburst that illuminated McNeil’s Nebula in Lynds 1630. It peaked in both 2003 and 2008 (Reipurth & Aspin, 2004; Aspin, 2011). This outburst is not cleanly classified as either an FU Orionis or an EX Lupi outburst (Muzerolle et al., 2005; Fedele et al., 2007; Contreras Peña et al., 2017a), but its detailed behavior is less relevant to the luminosity problem than the finding of an increase in its accretion rate.

Two other Orion protostars brightened well before the epoch of our comparison: Reipurth 50 and V883 Ori (Strom & Strom, 1993). The outburst of Reipurth 50 began between 1955 and 1979, and it has shown recent variability that has been interpreted as variable dust obscuration, not a change in its accretion rate (Chiang et al., 2015). The protostar was saturated in the IRAC and MIPS images used to define the HOPS sample, so it does not have a HOPS number. No precise estimate exists for the beginning of the V883 Ori outburst, but it has been bright since at least 1888 (Strom & Strom, 1993). It has the number HOPS 376.

6 THE BURST INTERVAL

Our search reveals two accretion bursts between the 2004–2005 Spitzer photometry and the 2010 WISE photometry: HOPS 223 and HOPS 383. Based on this finding, here we estimate the rate of accretion outbursts $r_{b}$ in a single protostar. This is

[TABLE]

where $n_{b}=2$ protostars is the number of bursts found within our time interval $\Delta t$, and $N_{p}$ is the number of protostars in our sample.

The time interval $\Delta t$ is from the first Spitzer observations to the last WISE observations incorporated in the AllWISE photometry. This varies depending on the location of each protostar within Orion. For 311 of the 319 protostars, the first Spitzer images were obtained in 2004 February or March, and the last relevant WISE images were obtained in 2010 September, giving baselines that range from 6.53 to 6.59 yr. For the remaining eight protostars, all in the Lynds 1622 region of Orion B, the first Spitzer images were acquired in 2005 October, reducing the interval to 4.90 yr. To an appropriate level of precision, $\Delta t=6.5$ yr for the vast majority of the sample, and we use this value for subsequent calculations.

We carried out the comparison for $N_{p}=319$ protostars. While 48 of them lack WISE counterparts, we interpret that as evidence against a burst. Even the HOPS 383 outburst, which reached a relatively low post-outburst luminosity of 7.5 $L_{\sun}$, was detected unambiguously in the WISE images. If we were to use only the 271 protostars with WISE counterparts in the denominator of Equation 1, we would derive a burst interval that is 15% shorter, which is a small change compared to the uncertainties discussed below. Inserting the values presented here into Equation 1 gives $r_{b}=9.65\times 10^{-4}$ yr*-1* and an interval between bursts of $I_{b}=1/r_{b}=1036.75$ yr.

To determine the uncertainty in the estimated interval, we model the bursts as a Poisson process over the time interval $\Delta t$ with a rate equal to $1/I_{b}$. If the burst interval is $I_{b}$ yr and we watch a protostar for $\Delta t$ yr, the chance we will detect one or more bursts from a given protostar is $1-e^{-\Delta t/I_{b}}$. The likelihood $P$ of observing exactly $n_{b}$ bursts in $N_{p}$ protostars is then a binomial distribution, where

[TABLE]

Adopting uniform priors, this gives the posterior probability $P\left(I_{b}|n_{b},N_{p},\Delta t\right)$ for the burst interval.

The top panel of Figure 12 plots $P$ against $I_{b}$ for our observational result of $n_{b}=2$ bursts in $N_{p}=319$ protostars over $\Delta t=6.5$ yr. The probability peaks at 27.2% for an interval of 1033.50 yr. This evaluation of the probability distribution function for $I_{b}$ gives a most likely interval that is very slightly less than the interval estimated above; henceforth, we will use 1040 yr. The 90% confidence interval, further illustrated with the cumulative probability distribution in the bottom panel, is from 690 to 40,300 yr; i.e, there is a 5% chance that the interval is below this range and a 5% chance that the interval is above it. For the $N_{p}=271$ case discussed above, where we exclude the Spitzer protostars without WISE counterparts, the most likely interval is 880 yr with a 90% confidence interval extending from 586 to 34200 yr, all about 15% shorter than the adopted values.

The broad 90% confidence interval is consistent with the findings of Hillenbrand & Findeisen (2015), who demonstrated that a sample of far more than 319 protostars is needed to precisely estimate the burst interval with high confidence given the available time intervals. If $I_{b}$ is of order one thousand years and the protostellar lifetime is $5\times 10^{5}$ yr (Evans et al., 2009), then a given protostar will undergo several hundred outbursts.

7 DISCUSSION

7.1 Dependence of Outburst Interval on Evolution

The typical burst interval of 1040 yr in Section 6 is much less than the interval calculated from results of similar analyses of Class II YSOs. Comparing WISE and Spitzer data for samples of several thousand primarily Class II YSOs in various star-forming regions, Scholz et al. (2013) used criteria similar to our own to find that the most likely burst interval is between 5000 and 50,000 yr. Carpenter et al. (2001) analyzed near-IR photometry spread across 2 yr of a region centered on the Orion Nebula Cluster with 2700 YSOs and detected no bursts. This puts a lower limit on the burst interval of 5400 yr. Though the upper end of our 90% confidence interval does overlap with the lower limits of these studies, the difference is intriguing.

These findings are consistent with outburst modeling where the disk instabilities are ultimately due to mass infall from the protostellar envelope. In these models most bursts would be expected to occur in Class 0 and Class I protostars, where the remaining envelope mass is significant, and bursts would be possible but rare for Class II YSOs, where the envelope is at most tenuous. Consistent with this outline, Vorobyov & Basu (2015) count bursts that occur in their gravitational instability simulations and find that most bursts occur in Class I protostars, while a smaller fraction take place in Class 0 protostars, and only a few happen in Class II YSOs.

7.2 Protostellar Outbursts and the Luminosity Problem

To consider the implications of our findings for the luminosity problem, we use a toy model to explore the fraction of a star’s ultimate main-sequence mass that is accumulated in outbursts. If the protostellar phase can be divided unambiguously into periods of quiescence and outburst, and the fraction of the time spent in outburst is $f_{b}$, then the total mass accumulated in quiescence is $M_{q}=\left(1-f_{b}\right)\int\dot{M_{q}}\,dt$. If the accretion rate during bursts is a constant factor $A$ greater than the slowly varying accretion rate during quiescence, then the mass accumulated during outbursts is $M_{b}=Af_{b}\int\dot{M_{q}}\,dt$. In both cases the integral is over the duration of the protostellar phase.

The fraction of the mass accumulated in bursts is then $M_{b}/\left(M_{b}+M_{q}\right)$, which is simply $Af_{b}/\left[\left(A-1\right)f_{b}+1\right]$. This quantity is plotted against $f_{b}$ for three values of $A$ in Figure 13. The fraction $f_{b}$ is $d_{b}/I_{b}$, the duration of each outburst divided by the interval between their onsets. We have argued that $I_{b}$ is of order 1000 years, but the $d_{b}$ of the outbursts we have examined are as yet unknown, since they are still ongoing. SOFIA observations of the HOPS 223 and 383 outbursts indicate that they were still at elevated luminosities about a decade after their onsets (W. J. Fischer et al., in preparation).

The outbursts included in our calculation of the burst frequency experienced luminosity increases of factors 12 (HOPS 223) and 35 (HOPS 383). These are lower limits to the factors by which the mass accretion rates increased, although we argued in Section 3 that the two factors are similar if the accretion luminosities are much greater than the stellar luminosities, as expected for protostars. If outbursts like that of HOPS 223 are typical, representing fifteen-fold increases in the mass accretion rate and persisting for 20 to 40 years ($0.02<f_{b}<0.04$), then between 23% and 38% of a star’s mass would be accumulated in outbursts. These values are larger than but not dramatically different from the fraction of 25% estimated by Offner & McKee (2011). If the HOPS 383 outburst is typical, then forty-fold increases in the mass accretion rate are more likely, and more than half of a star’s mass would be accumulated in outbursts for any $f_{b}>0.024$.

To investigate the luminosity problem, Fischer et al. (2017) plotted bolometric luminosities and temperatures for 315 Orion protostars. They found that the spread in luminosities at each bolometric temperature is three orders of magnitude, and that this spread is consistent with a slowly, exponentially decreasing accretion rate that depends on the final mass. Offner & McKee (2011) also found that models with an approximately constant accretion time and a limited role for outbursts could produce a broad spread in luminosities. Even if it is responsible for the accumulation of a nontrivial percentage of a star’s main-sequence mass, episodic accretion is not required to explain the spread in BLT diagrams or to solve the luminosity problem, which assumes a constant, mass-independent accretion rate.

Finally, if outbursts that result in luminosity increases of factors of 10 to 100 do occur with the frequency we calculate here, then these excursions are small compared to the factor of 1000 spread in protostellar luminosities. They do, however, imply that the path of each protostar across the BLT diagram does not vary smoothly with age and mass, and that stars with the same main-sequence mass may not have followed the same trajectory across the BLT diagram.

7.3 Other Constraints on the Protostellar Outburst Frequency

The separation of knots in protostellar outflows is an alternative way to estimate the frequency of enhanced accretion events. Arce et al. (2013) observed clumps in the HH 46/47 molecular outflow, which is driven by a Class I protostar. They inferred from the spacing and velocities of three clumps that they were ejected approximately 300 yr apart. If these are due to episodic ejection events that are driven by accretion, this points to an outburst interval in this protostar that is at the low end of the range we derive for the Orion protostars.

Another way to estimate the outburst frequency is to look for evidence of past outbursts in freeze-out chemistry (Visser et al., 2015). Enhanced luminosity alters the chemistry of the protostellar envelope, evaporating CO throughout the envelope and H2O to a radius 10$\times$ larger than usual. The CO evaporation enhances the abundance of HCO+ and reduces the abundance of N2H+; this altered chemistry persists for 103 to 104 yr. Jørgensen et al. (2013) used such arguments to conclude that the protostar IRAS 15398$-$3359 underwent an accretion outburst 100 to 1000 yr ago. The composition of circumstellar dust may also be altered by outbursts. Poteet et al. (2011) detected crystalline silicates in the spectrum of HOPS 68, the first such detection in a cold protostellar envelope. They argued that the heat necessary to crystallize the silicates must have been provided by an episodic accretion event.

A significant reduction in the uncertainty of the interval requires a timescale or a sample that is at least an order of magnitude larger (Hillenbrand & Findeisen, 2015), including protostars at distances of up to $\sim 2$ kpc instead of only 500 pc. Progress on this is being enabled by near-IR all-sky monitoring from the ground, e.g., with the Visible and Infrared Survey Telescope for Astronomy (VISTA) surveys (Contreras Peña et al., 2017b) and a United Kingdom Infrared Telescope Infrared Deep Sky Survey (UKIDSS) search for eruptive YSOs (Lucas et al., 2017). It will continue with the planned Wide-Field Infrared Survey Telescope (WFIRST).

While surveys at optical and near-IR wavelengths are valuable, mid-IR to submillimeter surveys provide more direct evidence of changes in the luminosities of the most deeply embedded protostars. Submillimeter monitoring with the James Clerk Maxwell Telescope is proving fruitful in understanding the occurrence rate and duration of outbursts (Herczeg et al., 2017; Johnstone et al., 2018). A sensitive far-IR space telescope with fast mapping capabilities such as the proposed Origins Space Telescope (Meixner et al., 2017) would probe deeply embedded protostars at the peaks of their SEDs. Billot et al. (2012) showed the value of far-IR monitoring in the vicinity of the Orion Nebula, where, for an initial look at the far-IR variability of protostars, they monitored 43 protostars with Herschel over six weeks to find that variations at the 10–20% level were typical at 70 and 160 µm. SOFIA is also essential for mid-to-far-IR follow-up of protostellar outbursts to determine their durations and to extend their SEDs across the mid- and far-IR regimes for accurate luminosities.

8 CONCLUSIONS

We compared three-band Spitzer and WISE photometry for 319 Orion protostars to search for luminosity outbursts in the $\sim 6.5$ yr between the acquisition of the two data sets. Although catalog photometry often indicates minor brightening in one or two bands, visual inspection reveals that these are usually false positives due to confusion in WISE. Those protostars that brightened by at least two standard deviations more than the mean difference at all three wavelengths of comparison were identified as outbursts and visually confirmed as such.

This approach turned up one previously known outburst, V2775 Ori (HOPS 223), and one new burst, HOPS 383, which was announced in our earlier paper (Safron et al., 2015). No other new sources were discovered. The latter is the first Class 0 accretion burst to be discovered. We estimate the most likely interval between bursts in a given protostar to be 1040 yr, with a 90% chance that the interval lies between 690 and 40,300 yr.

This and related studies of more evolved, Class II YSOs paint a picture where outbursts are frequent in the protostellar stage, occuring once every $\sim 10^{3}$ yr, and they are much less frequent in the optically revealed, disk dominated stage, occuring every $\sim 10^{4}$ yr. A similar dependence on evolutionary class was also reported by Barsony et al. (2005) and Contreras Peña et al. (2017b). There is emerging evidence that protostellar outbursts are less dramatic than Class II outbursts; those reported here represent brightenings of factors of 10–40 rather than the best known FU Orionis events observed in Class II YSOs, which feature brightenings by factors of 100 or more. These details are consistent with the predictions of models where the disk becomes unstable due to infall from a surrounding envelope, it fragments, and the accretion rate spikes when these fragments are accreted onto the star.

The duration of the bursts is not well constrained, but it appears to be typically longer than a decade. For a range of assumptions, episodic accretion events account for 25% or more of a star’s ultimate main-sequence mass, but they are not required to explain the luminosity distributions of large samples of protostars. Further surveys using a variety of observational techniques are needed to reduce the large uncertainty in the outburst interval.

The National Aeronautics and Space Administration (NASA) provided support for this work through awards issued by the Jet Propulsion Laboratory/California Institute of Technology (JPL/Caltech). Support for the contributions of WJF was provided in part by NASA through awards 03_0143 and 04_0134 issued by the Universities Space Research Association. This work is based in part on observations made with the Spitzer Space Telescope, which is operated by JPL/Caltech under a contract with NASA. It also makes use of data products from the Wide-field Infrared Survey Explorer, which is a joint project of the University of California, Los Angeles, and JPL/Caltech, funded by NASA. Finally, it has made extensive use of the NASA/IPAC Infrared Science Archive, which is operated by JPL/Caltech under a contract with NASA.