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

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 (Flaherty et al. 2012, 2013; Rice et al. 2012; Wolk et al. 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. Although the luminosity of a low-mass pre-main-sequence star is expected to change gradually over ∼10⁶ 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 (Wachmann 1939). Although 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 et al. 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 the 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 (Armitage et al. 2001; Lodato & Clarke 2004; Zhu et al. 2010). In 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, infrared (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–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.

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) showed that the 5σ 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 showed 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 spectral energy distribution (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 to 2005 campaigns discussed above. The data are published in Table 2 of Furlan et al. (2016), except 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σ 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.⁴ 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–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 Archive⁵ 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%–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.

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_acc, where ${L}_{\mathrm{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 η 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 η are expected to be small (≲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_acc,b + L_*)/(L_acc,q + L_*). This is a lower limit to the ratio of accretion luminosities, L_acc,b/L_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 L_☉ to 303 L_☉ 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⁻⁶ M_☉ yr⁻¹ 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 nondetection 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.

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.

Zoom In Zoom Out Reset image size

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 at 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 to 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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, because 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 makes this candidate a false positive.

Zoom In Zoom Out Reset image size

Figures 8 and 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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 nondetection of a real burst.

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 nondetection 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 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–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.

Zoom In Zoom Out Reset image size

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 L_☉ to ∼51 L_☉, 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⁻⁶ to 10⁻⁵ M_☉ yr⁻¹ 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 nondetection 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.

Zoom In Zoom Out Reset image size

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_☉ to 7.5 L_☉, 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 045 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 February. They found that HOPS 383 was a "possible" variable candidate, undergoing a brightness decrease of 2.66% ± 0.64% yr⁻¹. Imaging of the field at J, H, and K_s in 2015 December resulted in nondetections (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

Our search reveals two accretion bursts between the 2004 and 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

$\begin{eqnarray}&&{r}_{b}=\displaystyle \frac{{n}_{b}}{{N}_{p}{\rm{\Delta }}t},\end{eqnarray} \tag{ 1 }$

where n_b = 2 protostars is the number of bursts found within our time interval Δt, and N_p is the number of protostars in our sample.

The time interval Δ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, Δ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. Although 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_☉, 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 × 10⁻⁴ yr⁻¹ 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 Δt with a rate equal to 1/I_b. If the burst interval is I_b yr and we watch a protostar for Δt yr, the chance we will detect one or more bursts from a given protostar is $1-{e}^{-{\rm{\Delta }}t/{I}_{b}}$. The likelihood P of observing exactly n_b bursts in N_p protostars is then a binomial distribution, where

$\begin{eqnarray}&&\begin{array}{l}P\left({n}_{b}| {I}_{b},{N}_{p},{\rm{\Delta }}t\right)\\ \quad =\,\displaystyle \frac{{N}_{p}!}{{n}_{b}!({N}_{p}-{n}_{b})!}{\left(1-{e}^{-{\rm{\Delta }}t/{I}_{b}}\right)}^{{n}_{b}}{\left({e}^{-{\rm{\Delta }}t/{I}_{b}}\right)}^{({N}_{p}-{n}_{b})}.\end{array}\end{eqnarray} \tag{ 2 }$

Adopting a uniform prior for I_b, this gives the posterior probability $P\left({I}_{b}| {n}_{b},{N}_{p},{\rm{\Delta }}t\right)$ for the burst interval.

The top panel of Figure 12 plots $P\left({I}_{b}| {n}_{b},{N}_{p},{\rm{\Delta }}t\right)$ against I_b for our observational result of n_b = 2 bursts in N_p = 319 protostars over Δt = 6.5 yr. The probability peaks at 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 34,200 yr, all about 15% shorter than the adopted values.

Zoom In Zoom Out Reset image size

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 1000 yr and the protostellar lifetime is 5 × 10⁵ yr (Evans et al. 2009), then a given protostar will undergo several hundred outbursts.

7.1. Dependence of Outburst Interval on Evolution

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 yr, but the d_b of the outbursts we have examined are as yet unknown, because 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. 2019, in preparation).

Zoom In Zoom Out Reset image size

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–40 yr (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–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

We compared three-band Spitzer and WISE photometry for 319 Orion protostars to search for luminosity outbursts in the ∼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, occurring once every ∼10³ yr, and they are much less frequent in the optically revealed, disk-dominated stage, occurring every ∼10⁴ 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 W.J.F. was provided in part by NASA through awards SOF 03-0143 and SOF 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.