# The Rate, Amplitude and Duration of Outbursts from Class 0 Protostars in
Orion

Wafa Zakri  Ritter Astrophysical Research Center, Dept. of Physics and
Astronomy, University of Toledo, Toledo, OH 43606, USA  Department of Physics,
Jazan University, Jazan, Saudi Arabia  S. T. Megeath Ritter Astrophysical
Research Center, Dept. of Physics and Astronomy, University of Toledo, Toledo,
OH 43606, USA  William J. Fischer Space Telescope Science Institute, 3700 San
Martin Drive, Baltimore, MD 21218, USA  Robert Gutermuth Department of
Astronomy, University of Massachusetts, Amherst, MA 01003 USA  Elise Furlan
NASA Exoplanet Science Institute, Caltech/IPAC, Mail Code 100-22, 1200 E.
California Blvd., Pasadena, CA 91125, USA  Lee Hartmann Department of
Astronomy, University of Michigan, 1085 S. University Street, Ann Arbor, MI
48109, USA  Nicole Karnath SOFIA-USRA, NASA Ames Research Center, MS 232-12,
Moffett Field, CA 94035, USA  Mayra Osorio Instituto de Astrofisica de
Andalucia, CSIC, Camino Bajo de Huétor 50, E-18008 Granada, Spain  Emily
Safron Department of Physics and Astronomy, Lousiana State University, Baton
Rouge, LA 70803 USA  Thomas Stanke European Southern Observatory, Karl-
Schwarzschild-Str. 2, D-85748 Garching b. München, Germany  Amelia M. Stutz
Departamento de Astronomía, Universidad de Concepción,Casilla 160-C,
Concepción, Chile  Max-Planck-Institute for Astronomy, Königstuhl 17, 69117
Heidelberg, Germany  John J. Tobin National Radio Astronomy Observatory, 520
Edgemont Rd., Charlottesville, VA 22903, USA  Thomas S. Allen  Lowell
Observatory, 1400 W. Mars Hill Road, Flagstaff AZ 86001, USA  Sam Federman
Ritter Astrophysical Research Center, Dept. of Physics and Astronomy,
University of Toledo, Toledo, OH 43606, USA  Nolan Habel Ritter Astrophysical
Research Center, Dept. of Physics and Astronomy, University of Toledo, Toledo,
OH 43606, USA  P. Manoj Department of Astronomy and Astrophysics, Tata
Institute of Fundamental Research, Homi Bhabha Rd, Colaba, Mumbai 400005,
India  Mayank Narang Department of Astronomy and Astrophysics, Tata Institute
of Fundamental Research, Homi Bhabha Rd, Colaba, Mumbai 400005, India  Riwaj
Pokhrel Ritter Astrophysical Research Center, Dept. of Physics and Astronomy,
University of Toledo, Toledo, OH 43606, USA  Luisa Rebull Infrared Science
Archive (IRSA), IPAC, 1200 E. California Blvd., California Institute of
Technology, Pasadena, CA 91125, USA  Patrick D. Sheehan Center for
Interdisciplinary Exploration and Research in Astronomy, 1800 Sherman Rd.,
Evanston, IL 60202, USA  Dan M. Watson Department of Physics and Astronomy,
University of Rochester, Rochester, NY 14627, USA

###### Abstract

At least half of a protostar’s mass is accreted in the Class 0 phase, when the
central protostar is deeply embedded in a dense, infalling envelope. We
present the first systematic search for outbursts from Class 0 protostars in
the Orion clouds. Using photometry from Spitzer/IRAC spanning 2004 to 2017, we
detect three outbursts from Class 0 protostars with $\geq 2$ mag changes at
3.6 or 4.5 $\mu$m. This is comparable to the magnitude change of a known
protostellar FU Ori outburst. Two are newly detected bursts from the
protostars HOPS 12 and 124. The number of detections implies that Class 0
protostars burst every 438 yr, with a 95% confidence interval of 161 to 1884
yr. Combining Spitzer and WISE/NEOWISE data spanning 2004-2019, we show that
the bursts persist for more than nine years with significant variability
during each burst. Finally, we use $19-100$ $\mu$m photometry from SOFIA,
Spitzer and Herschel to measure the amplitudes of the bursts. Based on the
burst interval, a duration of 15 yr, and the range of observed amplitudes,
3–100% of the mass accretion during the Class 0 phase occurs during bursts. In
total, we show that bursts from Class 0 protostars are as frequent, or even
more frequent, than those from more evolved protostars. This is consistent
with bursts being driven by instabilities in disks triggered by rapid mass
infall. Furthermore, we find that bursts may be a significant, if not
dominant, mode of mass accretion during the Class 0 phase.

infrared: stars — stars: variables: general — stars: formation — stars:
protostars

##  1 INTRODUCTION

The outburst of HOPS 383, identified in multi-epoch Spitzer and WISE
observations, was the first detection of a luminous outburst toward a Class 0
protostar (Safron et al., 2015). This outburst appears similar to the
outbursts driven by rapid increases in mass accretion detected toward more
evolved protostars. Class 0 protostars sample the earliest phase of
protostellar evolution (Andre et al., 1993; Furlan et al., 2016), and more
than half of mass accretion is thought to occur in the Class 0 phase (Fischer
et al., 2017). If outbursts are common for Class 0 protostars, they may
contribute a substantial fraction of a star’s mass. Indeed, simulations of
protostars and modeling of chemical abundances in protostellar envelopes
predict that outbursts are frequent in the Class 0 phase (Bae et al., 2014;
Visser et al., 2015; Hsieh et al., 2019).

The luminosity functions of protostars also suggest bursts are a common, if
not dominant, mode of accretion for protostars. In particular, models of
protostellar evolution invoking isothermal sphere collapse require
luminosities higher than those observed. This discrepancy, known as the
protostellar luminosity problem, can be resolved if most of the stellar mass
is accreted during bursts (Kenyon et al., 1994; Dunham & Vorobyov, 2012).
There are alternative models of infall and accretion, however, that do not
require episodic bursts to reproduce the luminosity function (Offner & McKee,
2011).

Establishing the importance of bursts for mass accretion requires multi-epoch,
photometric surveys that provide direct measurements of their rate,
amplitudes, and durations. While ground-based 2 µm observations can detect
bursts from Class I and flat-spectrum protostars (e.g. Contreras Peña et al.,
2017), reliable detections of Class 0 protostars and their bursts requires
observations at $>2.5$ µm. The Spitzer Space Telescope (Werner et al., 2004)
detected all but the youngest, most deeply embedded protostars (Kryukova et
al., 2012; Stutz et al., 2013). The Spitzer survey of the Orion molecular
clouds in 2004-2005 provided baseline measurements of protostellar flux
densities from 3.6 to 70 $\mu$m (Megeath et al., 2012). An extension of the
Spitzer survey in 2008, still during the cryogenic mission, led to detections
of outbursts in the overlapping regions between the 2008 and 2004/2005 MIPS 24
$\mu$m maps. This revealed outbursts of the flat-spectrum protostar V2775 Ori
(HOPS 223; Caratti o Garatti et al., 2011; Fischer et al., 2012) and the Class
0 protostar HOPS 383 (Safron et al., 2015). Examination of the WISE survey
3.4–22 $\mu$m photometry in 2010 showed no additional bursts, leading to an
estimated interval between bursts of 1000 yr for each protostar (Fischer et
al., 2019).

In this paper, we search for outbursts from the 92 identified Class 0
protostars in the Orion clouds by combining the 2004 cryogenic data of Orion
with warm mission data taken during the YSOVAR exploration program and the
OrionTFE (The Final Epoch) survey, which duplicated the original survey fields
in 2016–2017. These data provide sparsely sampled 3.6 and 4.5 $\mu$m light
curves that span 13 yr with a higher angular resolution and sensitivity than
WISE/NEOWISE.

We announce the discovery of outbursts from the Class 0 protostars HOPS 12 and
124, and show the light curve covering the entire burst for the previously
known Class 0 protostar HOPS 383. We supplement the light curves with
additional Spitzer and WISE/NEOWISE photometry. We use 19–100 $\mu$m mid to
far-IR data to estimate the change in luminosity during each burst, with SOFIA
providing additional epochs to earlier Spitzer and Herschel observations.

##  2 Observations and Data Analysis

###  2.1 Spitzer Cryogenic and Warm Mission Data

Table 1: Spitzer/IRAC Photometry

| Flux Density (mJy)
---|---
Date | 3.6 $\micron$ | 4.5 $\micron$
HOPS 12
2004 Mar 9 | $0.45\pm 0.02$ | $2.63\pm 0.04$
2004 Oct 12 | $0.44\pm 0.02$ | $2.30\pm 0.04$
2009 Nov 111 | $1.65\pm 0.04^{2}$ | $14.2\pm 0.4^{2}$
2010 Nov 193 | $1.85\pm 0.05^{2}$ | $14.2\pm 0.7^{2}$
2017 Jan 11 | $1.71\pm 0.03$ | $15.6\pm 0.1$
2017 Jun 1 | $1.82\pm 0.03$ | $16.2\pm 0.1$
2019 Jul 21 | $0.92\pm 0.02$ | $9.8\pm 0.1$
HOPS 124
2004 Feb 17 | $3.59\pm 0.11$ | $16.9\pm 0.3$
2004 Feb 18 | $3.64\pm 0.11$ | $17.4\pm 0.2$
2004 Oct 8 | $2.65\pm 0.13$ | $13.6\pm 0.4$
2004 Oct 27 | $2.69\pm 0.11$ | $14.5\pm 0.3$
2017 Jan 1 | $18.7\pm 0.3$ | $122\pm 2$
2017 Jan 6 | $22.0\pm 0.3$ | $129\pm 2$
2017 Jun 1 | $18.3\pm 0.3$ | $112\pm 2$
2017 Jun 6 | $19.8\pm 0.4$ | $126\pm 1$
2019 Jul 15 | $13.4\pm 0.2$ | $84.0\pm 0.9$
HOPS 383
2004 Mar 9 | $\cdots$ | $0.30\pm 0.02$
2004 Oct 12 | $\cdots$ | $1.34\pm 0.03$
2009 Nov 121 | $\cdots$ | $6.17\pm 0.07^{2}$
2010 Nov 193 | $\cdots$ | $6.46\pm 0.20^{2}$
2017 Jan 11 | $\cdots$ | $1.23\pm 0.03$
2017 Jun 1 | $\cdots$ | $0.57\pm 0.03$
2019 Jul 15 | $\cdots$ | $0.32\pm 0.02$
2019 Jul 21 | $\cdots$ | $0.33\pm 0.02$

11footnotetext: Light curve of first YSOVAR campaign taken from 2009 Oct 23 to
Dec 01

22footnotetext: Median value of light curve $\pm$ median absolute deviation of
light curve

33footnotetext: Light curve of second YSOVAR campaign taken from 2010 Oct 29
to Dec 09

Spitzer initially surveyed Orion with the IRAC camera (Fazio et al., 2004)
during the cryogenic (hereafter: cryo) mission, observing most protostars in
two epochs separated by six months in 2004, with single-epoch field extensions
in 2007 and 2009 (Megeath et al., 2012; Kryukova et al., 2012). With the onset
of the warm mission, the YSOVAR exploration program observed the Orion Nebula
Cluster (ONC) during 80 epochs in 2009 and 2010 (Morales-Calderón et al.,
2011; Rebull et al., 2014). The Orion: The Final Epoch (hereafter: OrionTFE)
program then duplicated the cryo field in two epochs in 2016 and 2017.
Finally, during the Spitzer Beyond phase of the mission, Spitzer observed YSOs
in 2019 that had been found to vary by $\geq 1$ mag between the cryogenic and
OrionTFE programs.

We extracted photometry from images downloaded from the Spitzer Heritage
Archive and processed through a custom IDL program developed for cryo mission
data and then later adapted for the warm mission (Gutermuth et al., 2009;
Morales-Calderón et al., 2011). At the positions of YSOs taken from Megeath et
al. (2012, 2016), we performed aperture photometry on the individual frames
using the “aper” utility from the IDL Astronomy Users’ Library (Landsman,
1993). The aperture radius was $2.4^{\prime\prime}$ and the sky annulus
extended $2.4^{\prime\prime}$ to $7.2^{\prime\prime}$. The adopted zero points
in DN s-1 at the IRAC native pixel scale were 19.455 and 18.699 for the cryo
data and 19.306 and 18.669 for the warm data in the 3.6 and 4.5 $\mu$m bands,
respectively (Gutermuth et al., 2008; Morales-Calderón et al., 2011).
Magnitudes were converted to flux densities using the zero magnitude fluxes
from Reach et al. (2005). The Spitzer/IRAC photometry is shown in Table 1. We
present only median values of the YSOVAR photometry for each of the two Orion
campaigns; the light curves are available in Morales-Calderón et al. (2011).

For the comparison of the HOPS 12 Spitzer photometry to the WISE/NEOWISE data,
we degraded the Spitzer data from the cryo and warm mission to a resolution
similar to that of WISE. To do this, the Spitzer 3.6 and 4.5 $\mu$m mosaics
were convolved with a kernel to reproduce the larger WISE 3.4 and 4.6 $\mu$m
PSFs. We then extracted photometry from the Spitzer images with an aperture
radius of 12.4″ and a sky annulus extending from 12.4″ to 15.1″. The data were
calibrated by comparing the photometry of isolated stars in the convolved
images with the photometry of Megeath et al. (2012).

###  2.2 WISE and NEOWISE Data

We used WISE mission 3.4 and 4.6 $\mu$m data to supplement the Spitzer
photometry (Wright et al., 2010; Mainzer et al., 2011). For HOPS 12, the
photometry was obtained from the WISE All-Sky Database, the WISE 3-Band Cryo
Database, and the NEOWISE-R single-image photometry data. For the single-image
photometry, we averaged all points taken within a single week. To separate
HOPS 124 from the neighboring HOPS 125, we performed PSF fitting on unWISE
images for every individual epoch of the WISE and NEOWISE mission from 2010 to
2017 (Lang, 2014; Meisner et al., 2017). Using the DAOPHOT implementation in
the IDL Astronomy Users’ Library, we simultaneously fit the PSF of HOPS 124
and the overlapping PSFs of nearby point sources (Landsman, 1993).

###  2.3 SOFIA FORCAST Photometry

HOPS 12, HOPS 124, HOPS 383 and their surroundings were observed with FORCAST
aboard SOFIA in programs 04_0134 and 07_0200. HOPS 12 was observed at 25.3 and
37.1 µm on 2019 October 16 and at 19.7, 25.3, 31.5, and 37.1 µm on 2019
October 25. HOPS 124 was observed at the same four wavelengths on 2019 October
16. HOPS 383 was observed at 25.3 and 37.1 µm on 2016 February 11. The
observations were performed with two-position chopping and nodding 180∘ from
the chop direction (C2N/NMC mode).

We used the flux-calibrated (Level 3) images provided by the SOFIA Science
Center, with 0.768″ pixels. We obtained flux densities of the targets in 8 or
10 pixel (6.14″ or 7.68″) apertures with subtraction of the modal sky signal
in an annulus extending from 10 pixels to 12 pixels (7.68″ to 9.22″); these
were determined from inspection of the point-spread functions. The color
corrections needed to compare the FORCAST and MIPS 24 $\mu$m photometry are
described in Appendix A. The extraction of the MIPS 24 $\mu$m fluxes is
described in Megeath et al. (2012).

###  2.4 Far-IR Data from Herschel, Spitzer and SOFIA

Data from PACS on Herschel (Poglitsch et al., 2010; Pilbratt et al., 2010),
MIPS on Spitzer (Rieke et al., 2004), and HAWC+ on SOFIA (Harper et al., 2018)
measured the 70–100 $\mu$m far-IR fluxes of the three protostars before,
during, and/or after the bursts. Fischer et al. (2020) describe the
acquisition, data processing, and photometric measurement of the PACS 70 µm
images, which come from the Herschel Orion Protostar Survey (HOPS; Furlan et
al., 2016). The 100 µm PACS photometry used for HOPS 383 is from Herschel
Gould Belt Survey data (André et al., 2010) and was extracted using the
approach described in Stutz et al. (2013). For HOPS 12, we separated the 70
$\mu$m PACS fluxes of the blended eastern and western sources using the IDL
implementation of DAOPHOT (Landsman, 1993). The positions and fluxes of the
sources were then tweaked to minimize the artifacts in the subtracted images
as judged by eye. HOPS 12 was also observed pre-outburst by MIPS at 70 $\mu$m;
photometry for this source was obtained using the procedure of Fischer et al.
(2012), which uses 70 µm PACS data to calibrate 70 µm MIPS data.

HOPS 383 was imaged at 89 µm with 8″ resolution by HAWC+ aboard SOFIA in
program 07_0200. Observations were obtained in total-intensity scan-mapping
mode on 2019 September 10. We extracted photometry from the flux-calibrated
(Level 3) image provided by the SOFIA Science Center. We used an 11 pixel
(17.05″) aperture and subtracted the modal sky signal in an annulus extending
from 11 pixels to 16 pixels (17.05″ to 24.80″); these were determined from
inspection of the point-spread function of the brighter HOPS 90.

Figure 1: Spitzer 3.6 and 4.5 $\mu$m images of the HOPS 12 outburst (left,
$\alpha=83.78584^{\circ}$, $\delta=-5.93174^{\circ}$ J2000) and the HOPS 124
outburst (right, $\alpha=84.83327^{\circ}$, $\delta=-7.43644^{\circ}$ J2000)
with the locations of the protostars marked. The 2004 cyro images (top) show
the protostars pre-outburst, while the OrionTFE images (bottom) show the
protostars during the bursts.

The 19.7 to 100 µm fluxes from Spitzer, SOFIA, and Herschel are shown in Table
2.

Table 2: Photometry from 19.7 to 100 µm Wavelength | Flux Density | Instrument | Date
---|---|---|---
(µm) | (mJy) |  |
HOPS 12
19.7 | 276±22 | FORCAST | 2019 Oct 25
24 | 276±14 | MIPS | 2004 Mar 20
25.3 | 1130±60 | FORCAST | 2019 Oct 16
25.3 | 1140±90 | FORCAST | 2019 Oct 25
31.5 | 2560±30 | FORCAST | 2019 Oct 25
37.1 | 3160±140 | FORCAST | 2019 Oct 16
37.1 | 3170±200 | FORCAST | 2019 Oct 25
70 | 5780±580 | MIPS | 2004 Mar 20
70 | 15600±800 | PACS | 2010 Sep 09
HOPS 124
19.7 | 3170±60 | FORCAST | 2019 Oct 16
24 | 1710±90 | MIPS | 2005 Apr 02
25.3 | 15600±200 | FORCAST | 2019 Oct 16
31.5 | 38100±100 | FORCAST | 2019 Oct 16
37.1 | 57800±400 | FORCAST | 2019 Oct 16
HOPS 383
24 | 5.08±0.27 | MIPS | 2004 Mar 20
24 | 178±9 | MIPS | 2008 Apr 19
25.3 | 279±40 | FORCAST | 2016 Feb 11
37.1 | 569±79 | FORCAST | 2016 Feb 11
70 | 13100±200 | PACS | 2010 Sep 10
89 | 2160±170 | HAWC$+$ | 2019 Sep 10
100 | 25400±1400 | PACS | 2010 Oct 08

##  3 RESULTS

In a comparison of the cryo, OrionTFE, and YSOVAR data for 319 protostars, we
identified five protostars with bursts of $\geq 2$ mag in the 3.6 or 4.5
$\micron$ Spitzer bands. This includes the previously known outbursts from the
flat-spectrum (hereafter: FS) protostar V2775 Ori, or HOPS 223 (Caratti o
Garatti et al., 2011; Fischer et al., 2012), the Class I protostar HOPS 41
(Park et al., 2021, W. Zakri et al., in prep.), and the Class 0 protostar HOPS
383 (Safron et al., 2015). The remaining objects are newly identified
outbursts of the Class 0 protostars HOPS 12 and HOPS 124. The focus of this
paper is the Class 0 outbursts.

The 1-2.5 $\mu$m spectrum of V2775 Ori shows it to be one of the lowest
amplitude FU Ori outbursts known (Fischer et al., 2012). We thus adopt the 2
mag threshold to focus on bursts with amplitudes comparable to known FU Ori
outbursts in more evolved protostars. Furthermore, in the Orion sample, only
bursts show variations at 3.6 or 4.5 $\micron$ that exceed 2 mag. In a
companion paper, we will show that protostars with smaller amplitude
variations exhibit a more diverse range of behaviors, including bursts,
fluctuations and fades (W. Zakri et al., in prep.; see also Park et al.,
2021).

The outbursts of HOPS 12 and 124 are shown by the cryo and OrionTFE images in
Figure 1. HOPS 12 is found in the ONC and was consequently observed in two
month-long campaigns by the YSOVAR program (see Rebull et al., 2014). HOPS 124
is located in the Lynds 1641 cloud (Allen & Davis, 2008). Both protostars have
nearby YSOs that may be companions. HOPS 12 is separated by 6.46” (projected
separation of 2,510 au at 389 pc) from a YSO later resolved into a close
binary by ALMA and VLA data (Tobin et al., 2020). We refer to the outbursting
source as HOPS 12 West and the close binary as HOPS 12 East. HOPS 124 is 9.47”
(projected separation 3,770 au at 398 pc, distances for HOPS 12 and 124 from
Tobin et al., 2020) from the flat-spectrum protostar HOPS 125 (Furlan et al.,
2016). Spitzer resolves a scattered light nebula associated with HOPS 124 that
is also apparent in high angular resolution HST/NICMOS 1.60 and 2.05 $\mu$m
images (Habel et al., 2021). These two bursts were not identified by Fischer
et al. (2019) or Park et al. (2021) since HOPS 12 and 124 are blended with
their companions in the ALLWISE point source catalog. Furthermore, HOPS 124
was in the middle of its rise during the 2010 WISE observations (see Sec.
3.1).

###  3.1 The 3–5 $\mu$m Light Curves between 2004 and 2020

The light curves for the two new outbursts are shown in Figure 2. HOPS 12
brightened between 2004 and 2017 by 1.66 mag at 3.6 $\micron$ and by 2.09 mag
at 4.5 $\micron$, corresponding to factors of 4.6 and 6.8 in flux density. The
light curve shows that HOPS 12 rose rapidly between 2004 and 2009, sustained
its brightness through 2017, and faded by $\sim 0.5$ mag by 2019. Between 2004
and 2017, HOPS 124 brightened by 2.2 mag at 3.6 µm and 2.3 mag at 4.5 µm, or
factors of 7.5 and 8.2. HOPS 124 did not fall within the YSOVAR field and was
not observed between 2004 and 2016–2017 by Spitzer. The Spitzer Beyond
observations show a $\sim$ 0.5 mag fade between 2017 and 2019.

Figure 2: Left: 4.5 µm light curves of HOPS 12 (top) and HOPS 124 (bottom).
Included are the cryo observations from 2004, YSOVAR data from 2009–2010,
OrionTFE measurements from 2016–2017, and Spitzer Beyond observations from
2019. HOPS 124 did not fall within the YSOVAR field. Dashed lines show the
times of the MIPS, PACS, and FORCAST observations. Middle: Spitzer
$[3.6]-[4.5]$ colors as a function of time. Right: Spitzer (4.5 µm), WISE, and
NEOWISE (4.6 µm) light curves; these have been measured relative to the
averaged Spitzer 2016–2017 magnitude (for Spitzer) or the averaged magnitude
of the NEOWISE data bracketing the averaged date of the OrionTFE data (for
WISE/NEOWISE). The Spitzer photometry for the 2004, 2010, and the 2016–2017
epochs have been averaged. The HOPS 12 points are for the combined East and
West sources.

The middle panels in Figure 2 show the $[3.6]-[4.5]$ colors as a function of
time. They show that the colors of both protostars became $\sim 0.5$ mag
redder during the bursts. This is the opposite of what one would expect if the
brightening were due to a decrease in extinction, and we interpret this as
evidence that the bursts are due to rises in the luminosities of the
protostars.

The rightmost panels in Figure 2 show the relative changes in the Spitzer,
WISE and NEOWISE photometry. Relative changes are used to minimize offsets due
to the different filter profiles of Spitzer and WISE. The Spitzer data are
measured relative to the averaged 2016 and 2017 OrionTFE epochs. The
WISE/NEOWISE data are measured relative to the average of the two photometry
points bracketing the OrionTFE 2016–2017 averaged date. For HOPS 12, the
combined magnitudes of the East and West sources obtained from the convolved
Spitzer images (see Sec. 2.1) are displayed since the two objects are not
resolved by WISE.

The HOPS 12 NEOWISE data show that by 2015, the source started fluctuating
during its burst with peaks in brightness separated by intervals of
approximately two years. For HOPS 124, the data show that the source is rising
in 2010 during the WISE observations, and is approximately midway between the
2004 and 2015-2019 mags. It is not clear if the WISE observations caught a
rapidly rising burst, or whether the source had a more gradual increase. The
protostar peaked in 2015, and similar to HOPS 12, fluctuated during its burst.

In summary, both bursts have lasted for 9 yr, although with a potentially
slower rise for HOPS 124. The light curves show $>0.5$ mag variability during
the bursts and that both sources may be declining by 2019.

The 4.5 µm light curve of the previously known outburst for the Class 0
protostar HOPS 383 is shown in Figure 3. Since it is located in the ONC, the
source was covered by the YSOVAR data. Its light curve encompasses the entire
burst event, which had a duration of $\sim 15$ yr. The observed decline is
consistent with the report by Grosso et al. (2020), based on NEOWISE data,
that the outburst of HOPS 383 has ended. This protostar shows the largest
amplitude change of our Class 0 protostar sample; the change of 3.6 mag at 4.5
$\mu$m corresponds to a factor of 22 change in flux density.

###  3.2 Measuring the Burst Amplitude

Observations at 19-100 $\mu$m with MIPS, PACS, FORCAST, and HAWC+ measure
changes in the luminosities of the outbursting protostars. While the 4.5
$\mu$m data trace scattered light from the disks of the protostars, the 19-100
$\mu$m data cover the reprocessed luminosities of the protostars (Whitney et
al., 2003). The fluxes at these wavelengths are more tightly correlated with
the luminosities of protostars, with the correlation improving with increasing
wavelength (Dunham et al., 2008). The shape of the protostellar SED over these
wavelengths changes slowly with increasing luminosity, with the peak of the
SED trending to shorter wavelengths (Kenyon et al., 1993; Furlan et al.,
2016). Hence, to reasonable approximation, the fluxes in this wavelength range
scale linearly with changes in luminosity, particulary those in the 70-100
$\mu$m range.

Figure 3: Comparison of HOPS 383 and nearby protostars in the 70 µm
Herschel/PACS image obtained on 2010 September 10 (left) and the 89 µm
SOFIA/HAWC+ image obtained on 2019 September 10 (middle). The location of HOPS
383 is marked ($\alpha=83.87421^{\circ}$, $\delta=-4.99754^{\circ}$ J2000).
Right: Light curve of HOPS 383 at 4.5 µm. The times of the MIPS, FORCAST, PACS
and HAWC$+$ observations are shown by the vertical lines.

The 19-100 $\mu$m measurements taken before or after bursts can be compared to
those obtained during bursts to estimate the amplitude of each burst, i.e. the
multiplicative factor by which the flux changes. This amplitude is
approximately equal to the amplitude of the luminosity change. The vertical
lines in Figure 2 show when the 19-100 $\mu$m measurements were made relative
to the bursts tracked by the 4.5 $\mu$m light curves.

The pre-outburst MIPS 70 $\mu$m photometry of HOPS 12 shows a flux of 5.78 Jy,
compared to the PACS flux of 15.56 Jy during the burst. Our deblending of HOPS
12 East and West sources with the PACS data show that 14% of the 70 $\mu$m
flux during the burst, i.e. 2.2 Jy of the total flux, comes from HOPS 12 East.
If we assume that HOPS 12 East did not vary between the two epochs, then the
pre- and post-outburst fluxes of HOPS 12 West are 3.6 and 13.4 Jy. These show
an increase in the flux by a factor of 3.7 at 70 $\mu$m.

HOPS 12 was also measured by MIPS at 24 $\mu$m pre-outburst (2004) and with
FORCAST during the burst (2019; Table 1). Correcting for the differences in
the SOFIA and Spitzer bandpasses (see Appendix A), we find the flux increased
by a factor of 2.5 at 24 $\mu$m. The increase is smaller than that measured by
the 70 $\mu$m data; this is likely due to the decline in the burst by 2019
when the FORCAST data were obtained. We therefore use the PACS value of 3.7 as
the amplitude of the luminosity change for this source.

HOPS 124 was also measured by MIPS at 24 $\mu$m pre-outburst (2005) and with
FORCAST during the burst (2019; Table 2). Again correcting for the different
filter wavelengths (Appendix A), we find the flux increased by a factor of 5.4
at 24 $\mu$m. We use this factor as the amplitude of the luminosity change;
however, since HOPS 124 was declining in 2019, this is likely an underestimate
of the peak burst amplitude.

For HOPS 383, we compared Herschel/PACS data during the burst in 2010 to
SOFIA/HAWC+ data taken after the burst in 2019 (Table 2, Furlan et al., 2016).
Linearly interpolating between the 70 µm and 100 µm PACS measurements gives an
89 µm flux density of $20.9$ Jy. Comparing that to the 89 µm HAWC+ data
(Figure 3), we find that the source faded by a factor of 9.7.

This amplitude is less than the factor of 35 increase at 24 $\mu$m between
2004 and 2008 reported by (Safron et al., 2015). In 2016, FORCAST observed
HOPS 383 and found that the flux increase was a factor of 40 higher than 2004
(Appendix A). The discrepency between the 24-25 $\mu$m amplitudes and 70-89 µm
amplitude is not understood. The outburst may not have returned to the 2004
levels during the HAWC+ observations. Alternatively, the 24 $\mu$m fluxes may
have been boosted by changes in the disk/envelope during the outburst. We
adopt the 70-89 µm amplitude of 9.7 since it is more conservative and likely
more representative of the luminosity change during the outburst.
Nevertheless, the FORCAST observations shows that the outburst was still
ongoing in 2016.

##  4 Discussion

###  4.1 The Burst Interval

Among the 92 Class 0 protostars in the Orion sample (Furlan et al., 2016), we
detected three outbursts with $\geq 2$ mag increases in the Spitzer 3.6 or 4.5
$\mu$m bands. Only two bursts of this magnitude are detected among the 227
more evolved Class I and flat-spectrum protostars in Orion: HOPS 41 in the ONC
and V2775 Ori in Lynds 1641 (see Sec 3, locations in Furlan et al., 2016).
Therefore, the Orion Class 0 protostars have a higher incidence of outbursts
than their older counterparts. We infer an outburst rate of (3 bursts) / (92
protostars) / (13 yr) = $2.5\times 10^{-3}$ yr-1 for Class 0 protostars. This
corresponds to a burst interval (interval $=$ 1/rate) of one burst every 400
yr for each protostar.

Figure 4: Probability density functions (top) and cummulative probability
distributions (bottom) for the burst interval. These are calculated for the
entire clouds sample (left: 92 Class 0 protostars and 227 Class I/FS
protostars), the YSOVAR field of the ONC (middle: 27 Class 0 protostars and 45
Class I/FS protostars), and all protostars in the clouds and in the ONC
(right: 319 protostars in the clouds and 72 in the ONC). A 13 year baseline is
assumed. These are based on equation 2 from Fischer et al. (2019) (also see
Hillenbrand & Findeisen, 2015; Contreras Peña et al., 2019).

From the number of bursts, the number of protostars, and the observed time
interval, a probability distribution can be derived for the burst interval
(Scholz et al., 2013; Contreras Peña et al., 2019; Hillenbrand & Findeisen,
2015). Following the approach of Fischer et al. (2019), we present the
posterior probability density functions and cumulative probability
distributions for the burst intervals in Figure 4. These values assume a
logarithmic prior; the details of our approach and the dependence of the
results on the adopted prior are discussed in Appendix B.

We calculate burst intervals for three pairs of samples. The first pair is (i)
all the Class 0 protostars and (ii) all the more evolved Class I/FS
protostars. Since the HOPS 383 outburst subsided by the time of the OrionTFE
survey, however, it was only detected since it fell within the YSOVAR
observations of the ONC. Thus, protostars outside the ONC that burst after
2004 and declined before 2016 would be missed in our survey. For this reason,
the second pair is (i) the Class 0 protostars in the ONC and (ii) Class I/FS
prototstars in the ONC. This pair has more complete time coverage, but
contains a smaller sample of protostars than the first pair. The final pair is
(i) all protostars in the entire field and (ii) all protostars in the ONC.

Table 3: Burst Intervals Sample | Median (yr) | 95% CI (yr)aa95% confidence interval using logarithmic prior
---|---|---
Class 0 | 438 | 161 – 1884
Class I/FS | 1741 | 527 – 12011
Class 0 in ONC | 199 | 60 – 1393
Class I/FS in ONCbbThese values are strongly dependent on the prior due to the detection of only one burst. | 823 | 155 – 21694
All protostarsccClass 0, Class I and Flat Spectrum combined | 879 | 400 – 2515
All protostars in ONCccClass 0, Class I and Flat Spectrum combined | 341 | 126 – 1468

The median values and 95% confidence interval from the cumulative
distributions are given in Table 3. Our burst interval for all protostars is
consistent with previous estimates of 1000 yr. The previous estimates used a
combination of Spitzer and WISE photometry (Fischer et al., 2019) and, more
recently, NEOWISE photometry (Park et al., 2021) to infer a burst interval.
Although the burst intervals from all three studies are consistent, each used
different criteria, and the previous studies used lower magnitude thresholds
for bursts than our 2 mag cut. Furthermore, the previous studies relied on
WISE/NEOWISE data which, primarily due to the lower angular resolution of
WISE, missed the HOPS 12 and HOPS 124 bursts. The median burst interval in
Table 3 ($P(t<I_{b})=0.5$) shrinks to 341 yr for the ONC sample alone due to
the better time sampling with the inclusion of the YSOVAR data. Thus, the 1000
yr burst interval may be an overestimate, although due to the small size of
the ONC sample, it still falls within the 95% confidence interval.

These data also show that Class 0 protostars have a median burst interval of
438 yr, a factor of four shorter than that for the more evolved Class I/FS
protostars. The interval for Class 0 protostars shrinks to 199 yr for the ONC
sample, still a factor of four shorter than that of more evolved protostars in
the ONC. To test the significance of these differences, we calculate the
p-values assuming $t_{b}=879$ yrs using Equation B1. For the Class 0
protostars we calculate the probability that the number of bursts is greater
than or equal to the number of observed bursts; these are $P(n_{burst}\geq
3)=15\%$ and $P(n_{burst}\geq 2)=6\%$ for the clouds and ONC sample,
respectively. For the Class I/FS protostars we calculate the values for the
number of bursts less than or equal to the observed numbers, these are
$P(n_{burst}\leq 2)=35\%$ and $P(n_{burst}\leq 1)=86\%$. If the burst interval
is 342 yrs, as suggested for all protostars using the ONC sample, then
$P(n_{burst}\geq 2)=27\%$ and $P(n_{burst}\leq 1)=49\%$ for the ONC. Thus, we
cannot rule out that protostars of all classes have the same burst interval.

There are caveats to this analysis. Due to degeneracies between inclination,
foreground extinction, and evolution, the division into samples based on SEDs
may place a few evolved protostars into the Class 0 sample and vice versa.
Based on an analysis of models presented in Stutz et al. (2013) and Furlan et
al. (2016), Stutz & Kainulainen (2015) found that extinction introduces
contamination at less than the 10% level. It is also possible that some Class
0 protostars are more evolved protostars seen at high inclinations. Only 4 of
75 Class 0 protostars observed in HST 1.6 $\mu$m observations, however, show
edge-on inclinations (Habel et al., 2021), and none of the outburst sources
appear edge-on. Thus, contamination of the Class 0 sample should have a
smaller effect than the statistical uncertainties on the burst interval. A
more important caveat is the assumption that all protostars have similar burst
rates. This requires that bursts do not result from multiple, distinct
mechanisms and that the properties of the protostar, such as binarity, do not
affect outbursts (Reipurth & Aspin, 2004; Scholz et al., 2013; Riaz et al.,
2021). More sophisticated treatments will require larger samples and
spectroscopic data.

Our measurement of the burst interval for Class 0 protostars can be compared
to indirect estimates from the literature. Hsieh et al. (2019) argued that
N2H+ and HCO+ molecules are depleted in the inner envelopes of protostars
during bursts as the CO and H2O snowlines move outward, and then reform after
the bursts as the snowlines recede. They predict a much longer burst interval
of 2400 yr from the fraction of Perseus Class 0 protostars with depletion. Our
burst interval is similar to the $310\pm 150$ yr interval implied by the
spacing between ejection events in the outflow from the Class 0 protostar
CARMA-7, located in the Serpens South Cluster (Plunkett et al., 2015). It is
also shorter than the 1-3 kyr interval between ejection events for outflows
driven by protostars and pre-main sequence stars in nearby clouds (Ioannidis &
Froebrich, 2012; Froebrich & Makin, 2016); this is consistent with a longer
interval for more evolved sources. These data are evidence for a connection
between the outbursts and ejection events in outflows.

###  4.2 The Fraction of Mass Accreted During Bursts

A primary goal of studies of episodic accretion is to constrain the fraction
of the mass accretion that occurs during bursts ($M_{b}/M$). We can estimate
that ratio for the Class 0 phase using the interval between bursts ($I_{b}$),
the duration of bursts ($d_{b}$), and the ratio of the accretion rate during
bursts to the accretion rate during quiescent intervals between bursts ($A$).
From Fischer et al. (2019),

$M_{b}/M=\frac{Ad_{b}/I_{b}}{(A-1)d_{b}/I_{b}+1}.$ (1)

We calculate the fraction for the Class 0 protostars using the rate of bursts
determined above, $I_{b}=438$ yr and a burst duration of $d_{b}=15$ yr. For
$A$, we adopt burst amplitudes of 3.7 for HOPS 12, 5.4 for HOPS 124, and 9.7
for HOPS 383; this assumes that all the luminosity is coming from accretion
and the accretion luminosity scales with the bolometric luminosity. The
resulting $M_{b}/M$ values are 12, 16 and 26%. If we reduce $I_{b}$ to the ONC
value of 199 yr, the $M_{b}/M$ values are 23, 31 and 44%. If we increase
$I_{b}$ to 1884 yr, the upper confidence interval value for the full sample,
then $M_{b}/M$ are 3, 4 and 7%.

Since the intrinsic luminosities of the protostars also contribute to the
observed luminosities, the ratio of pre-outburst to post-outburst luminosities
is a lower bound on $A$. If the entire pre-outburst luminosity of HOPS 383,
0.5 L⊙, is due to its intrinsic luminosity, then $M_{b}/M=$ 100%. From the
$L_{bol}$ values in Furlan et al. (2016) and the burst amplitudes, we estimate
the pre-outburst luminosities of HOPS 12 and HOPS 124 are approximately 2 and
9.5 L⊙. If their intrinsic luminosities contribute 1 L⊙ to their total
luminosities, then $M_{b}/M$ may be as high as 39% and 33%, respectively,
assuming $I_{b}$ is 199 yr. Thus, $M_{b}/M$ may be as low as 3–39% for HOPS
12, or as high as 7–100% for HOPS 383.

These data show bursts can play a significant role in stellar mass assembly,
particularly in the Class 0 phase. Variability with lower amplitudes than
those considered here may also play a role. The JCMT transient survey measured
the contribution of lower amplitude accretion events over a 4 yr interval (Lee
et al., 2021). These events contributed about 12% percent of the total mass
accreted onto Class 0 protostars. (Note that Lee et al. (2021) define the
fraction of mass accreted during bursts as the excess accretion over the
quiescent level divided by the total mass accreted. Since we define the
fraction of mass accreted as all the mass accreted during bursts divided by
the total accreted mass, we have revised their estimate to be consistent with
our definition.) We will explore the contribution of lower amplitude bursts in
a subsequent paper (W. Zakri et al., in prep.).

###  4.3 Implications for Burst Mechanisms

Outbursts are thought to result from the steady accumulation of mass in disks
when the rate of infall from the protostellar envelopes onto the disks exceeds
the rate of accretion from the disks onto the protostars. As the disks grow
more massive, gravitational and/or thermal instabilities lead to episodes of
rapid accretion onto the central protostar (Bell & Lin, 1994; Bell et al.,
1995; Vorobyov & Basu, 2015; Zhu et al., 2009; Bae et al., 2014).

The high rate of bursts of Class 0 protostars supports this picture since the
rate of mass infall onto the disks is highest during the Class 0 phase (Furlan
et al., 2016; Fischer et al., 2017). Furthermore, VLA and ALMA observations by
Tobin et al. (2020) find that disks around Class 0 protostars are
systematically more massive than those of Class I/FS protostars, hinting that
they are more likely to become gravitationally unstable. Indeed, Tobin et al.
(2020) find a Toomre-Q value for HOPS 124 indicative of an unstable disk.
Sharma et al. (2020) also finds the HOPS 383 disk is close to the stability
limit; this was based on VLA observations made in 2016 just before the
outburst started subsiding.

Close binary companions have also been proposed as triggers of outbursts
(Reipurth & Aspin, 2004; Muzerolle et al., 2013; Riaz et al., 2021); although
this mechanism may still require the replenishment of the disk by infall
(Green et al., 2016). So far there is no indication of triggering by binaries.
Only one of the outbursting Class 0 protostars, HOPS 124, show potential
evidence for a close companion; in this case, substructure in its disk
resolved by ALMA (Sheehan et al., 2020). Furthermore, since Class 0 protostars
do not show a higher incidence of close binaries, triggering by binaries
cannot explain a higher rate of outbursts in the Class 0 phase (Tobin et al.,
2021).

##  5 Conclusions

1\. We compared Spitzer mid-IR photometry from 2004 to 2017 for 92 Class 0
Orion protostars in Orion and detected 3 bursts of $\geq 2$ mag in the 3.6 or
4.5 µm bands. The bursts of the prototstars HOPS 12 and HOPS 124 are new
detections. These join the known outbursts from the Class 0 protostar HOPS
383, and the Class I/FS protostars HOPS 41 and V2775 Ori.

2\. The median probability time interval between bursts is 199-438 yr for the
Class 0 protostars; this interval is a factor of 4 less than that estimated
for Class I/FS protostars. We cannot rule out, however, that all protostars
have a burst interval of 860 yr, as inferred for the combined sample of 5
protostars. Thus, we find suggestive but not conclusive evidence that the
burst interval is smaller for Class 0 than for Class I/FS protostars. This
high rate of bursts can be driven by the rapid infall onto to the disks of the
Class 0 protostars.

3\. The combined Spitzer, WISE, and NEOWISE photometry show that the bursts
last at least 9 yr. The entire HOPS 383 burst lasted 15 yr, the remaining 2
Class 0 bursts did not end in 2019. There is significant variability during
the bursts.

4\. With SOFIA, Spitzer, and Herschel data from 19 to 89 µm, we estimate burst
amplitudes for the three protostars. Together with the burst interval and
duration, these suggest a significant fraction of the mass during the Class 0
phase is accreted during bursts. The fractions vary from source to source, and
may be as low as 3% for HOPS 12 and as high as 100% for HOPS 383. While bursts
are not required to resolve the protostellar luminosity problem, they can play
a significant role in mass accretion onto low-mass protostars.

We thank the referee for their insightful comments. This work uses
observations from the Spitzer Space Telescope, operated by JPL/Caltech under a
contract with NASA. This paper also uses data from the Wide-field Infrared
Survey Explorer, a joint project of the University of California, Los Angeles,
and JPL/Caltech, funded by NASA. Observations were also made with the NASA/DLR
Stratospheric Observatory for Infrared Astronomy (SOFIA). SOFIA is jointly
operated by the Universities Space Research Association, Inc. (USRA), under
NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR
contract 50 OK 0901 to the University of Stuttgart. Herschel is an ESA space
observatory with science instruments provided by European-led Principal
Investigator consortia and with important participation from NASA. Finally,
this work makes use of the NASA/IPAC Infrared Science Archive, operated by
JPL/Caltech under a contract with NASA. STM and RAG were supported by the NASA
ADAP grant 80NSSC19K0591, and STM was supported by the NASA ADAP grant
80NSSC20K0454. RP was supported by the NASA ADAP grant 80NSSC18K1564. Support
for WJF was provided by NASA through award #07_0200 issued by USRA. AS
gratefully acknowledges funding support through Fondecyt Regular (project code
1180350), from the ANID BASAL project FB210003, and from the Chilean Centro de
Excelencia en Astrofísica y Tecnologías Afines (CATA) BASAL grant AFB-170002.
M.O. acknowledges support from the Spanish MINECO/AEI AYA2017-84390-C2-1-R
(co-funded by FEDER) and MCIN/AEI/10.13039/501100011033 through the
PID2020-114461GB-I00 grant, and from the State Agency for Research of the
Spanish MCIU through the “Center of Excellence Severo Ochoa” award for the
Instituto de Astrofísica de Andalucía (SEV-2017-0709). This work was completed
while STM was a Fulbright Scholar hosted by AS at the Universidad de
Concepcíon. The National Radio Astronomy Observatory is a facility of the
National Science Foundation operated under cooperative agreement by Associated
Universities, Inc.

##  Appendix A Color Corrections

Figure 5: The transmission profiles for the FORCAST 19.7 $\mu$m filter,
FORCAST 25.3 $\mu$m filter, FORCAST 31.7 $\mu$m filter and MIPS 24 $\mu$m
filter. The IRS spectrum of HOPS 124 is also displayed. All have been
normalized so that the peak value is unity.

To compare the FORCAST 19.7 $\mu$m, 25.3 $\mu$m and 31.7 $\mu$m photometry to
the MIPS 24 $\mu$m photometry, it is essential to correct for the different
wavelength ranges of the three bandpasses, all of which overlap with the MIPS
bandpass (Figure 5). This must take into account the spectral shape of the
observed protostars, the adopted spectral shapes for the calibration of the
FORCAST and MIPS photometry, and the nominal wavelengths adopted for each
filter. We estimate the correction by assuming that the burst does not change
the spectral shape of the protostar. We use the IRS spectra of HOPS 12 and
HOPS 124 (Figure 5; Furlan et al., 2016). For HOPS 383, which does not have an
IRS spectrum, we use the best fit radiative transfer model of the SED from
Furlan et al. (2016). We start by determining the flux density in each of the
bandpasses based upon the IRS spectrum. In the FORCAST photometry, the
calibration of the delivered data assumes a flat-spectrum source. Since the
units of flux density for the calibrated FORCAST photometry and the IRS
spectrum are in Janskys (i.e., they are in $F_{\nu}$ instead of
$F_{\lambda}$), yet the IRS spectra are parameterized by wavelength, we use
the following equation for the FORCAST flux density

$F^{\rm IRS/FOR}_{\nu}(\lambda_{0})=\frac{\int(1/\lambda)F^{\rm
IRS}_{\nu}(\lambda)T(\lambda)d\lambda}{\int(1/\lambda_{0})T(\lambda)d\lambda},$
(A1)

where $F^{\rm IRS/FOR}_{\nu}(\lambda_{0})$ is the inferred flux density for
the FORCAST filter, $F^{\rm IRS}_{\nu}(\lambda)$ is the flux density from the
IRS spectra, $T(\lambda)$ is the transmission in electrons per photon, and
$\lambda_{0}$ is the adopted nominal wavelength for the filter (i.e., 19.7,
25.3 or 31.7 $\mu$m). The MIPS calibration assumes a source with a Rayleigh-
Jeans spectrum (Engelbracht et al., 2007). The resulting flux density inferred
for MIPS at 23.675 $\mu$m (the nominal wavelength) is

$F^{\rm IRS/MIPS}_{\nu}(\lambda_{0})=\frac{\int(1/\lambda^{2})F^{\rm
IRS}_{\nu}(\lambda)R(\lambda)d\lambda}{\int(\lambda_{0}^{2}/\lambda^{4})R(\lambda)d\lambda},$
(A2)

where $F^{\rm IRS/MIPS}_{\nu}$ is the flux density from the MIPS observations
and $R(\lambda)$ is the response in electrons per energy.

Assuming the shape given by the IRS spectrum remains constant, then the ratio
of the two flux densities can be used to estimate the contemporaneous MIPS
flux densities from the actual FORCAST photometry. The estimated MIPS flux
density is given by

$F^{\rm FOR/MIPS}_{\nu}(23.675~{}\mu m)=\frac{F^{\rm
IRS/MIPS}_{\nu}(23.675~{}\mu m)}{F^{\rm IRS/FOR}_{\nu}(25.3~{}\mu m)}F^{\rm
FOR}_{\nu}(25.3~{}\mu m),$ (A3)

where $F^{\rm FOR}_{\nu}(25.3~{}\mu m)$ is the observed flux density. The
burst amplitude is then given by

$A=\frac{F^{\rm FOR/MIPS}_{\nu}(23.675~{}\mu m)}{F^{\rm
MIPS}_{\nu}(23.675~{}\mu m)}.$ (A4)

Similar equations can be used with the observed flux densities in the two
other FORCAST filters.

Using the flux densities in Table 2 for the FORCAST and MIPS data, we estimate
the burst amplitudes in Table 4. While the amplitudes determined from the 25.3
$\mu$m and 31.7 $\mu$m photometry are consistent, the amplitudes derived from
the 19.7 $\mu$m filter are higher. This is likely due to the breakdown of our
assumption that the pre-outburst spectrum approximates the spectrum during the
burst at this wavelength. Given that the 25.3 $\mu$m band is closest to the
the nominal wavelength of the MIPS 24 $\mu$m filter, we adopt these values in
our analysis of the bursts.

Table 4: FORCAST Burst Amplitudes Filter | 19.7 $\mu$m | 25.3 $\mu$m | 31.7 $\mu$m
---|---|---|---
HOPS 12 | 3.8 | 2.5 | 2.5
HOPS 124 | 8.5 | 5.4 | 5.1
HOPS 383 |  | 40 |

##  Appendix B Model and Posterior Probability Distribution

Table 5: Burst Intervals as a Function of Prior Sample | Uniform (yrs)aamedian value and 95% confidence limits in parentheses | Logarithmic (yrs)aamedian value and 95% confidence limits in parentheses | Jeffries (yrs)aamedian value and 95% confidence limits in parentheses
---|---|---|---
Class 0 | 701 (248 - 3310) | 438 (161-1884) | 371 (167 - 1086)
Class I/FS | 4203 (979 - 51391) | 1741 (527-12011) | 1350 (531 - 5130)
Class 0 in ONC | 487 (113 - 6502) | 199 (60-1393) | 155 (61 - 590)
Class I/FS in ONC | 22128 (721 - 366828) | 823 (155-21694) | 485 (148 - 3287)
All protostars | 1120 (531 - 3011) | 879 (400-2515) | 796 (418 - 1799)
All protostars in ONC | 546 (193 - 2579) | 341 (126-1468) | 289 (130 - 846)

Fischer et al. (2019) presented a model for the number of outbursts where the
probability of the number of bursts, $n_{b}$ is given by

$P(n_{b}|I_{b},\Delta
T,N_{p})=\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})}.$ (B1)

In this model, $I_{b}$ is the interval between bursts ($I_{b}=1/R_{b}$ where
$R_{b}$ is the rate of bursts), $N_{P}$ is the number of protostars, and
$\Delta T$ is the time interval over which bursts were monitored. The
probability of detecting one or more bursts for a given protostar is
$(1-e^{-\Delta T/I_{b}})$; this value goes to one as $\Delta T/I_{b}$ goes to
infinity. In the limit $\Delta T/I_{b}<<1$, this model can be approximated by
a binomial distribution (Hillenbrand & Findeisen, 2015) or Poisson
distribution (Contreras Peña et al., 2019). We have determined the posterior
probability distribution for $P(n_{b})$ for three different uninformative
priors: a uniform prior, a logarithimic prior, $P(I_{b})\propto 1/I_{b}$, and
a Jeffries prior, $P(I_{b})\propto 1/I_{b}^{1.5}$ when $\Delta T/I_{b}<<1$. In
all cases, the prior probability drops to 0 when $I_{b}>0.5$ Myr, the lifetime
of a protostar.

In Table 5, we compare the median value and 95% confidence intervals obtained
from the 2.5% and 97.5% points of the cumulative distributions. The uniform
weighting increases the median value of $I_{b}$; this prior increases by a
factor of ten for each successive log10 interval of $I_{b}$, and hence it is
weighted to longer values of $I_{B}$. In contrast, the Jeffries prior gives
the shortest values of $I_{b}$. Since it is appropriate to use a logarithmic
prior for determining the order of magnitude of the burst interval, we adopt
that prior in the paper. We find that the logarithmic and Jeffries priors
produce similar posterior probabilty distributions. Furthermore, the median
$I_{b}$ values for the uniform prior are within the logarithmic prior’s
confidence intervals. Only for the Class I/FS burst interval in the ONC, which
is dependent on one detected burst, is the posterior probability distribution
driven by the prior. In this case, $I_{b}$ is poorly constrained.

## References

  * Allen & Davis (2008) Allen, L. E., & Davis, C. J. 2008, Low Mass Star Formation in the Lynds 1641 Molecular Cloud, ed. B. Reipurth, Vol. 4, 621
  * Andre et al. (1993) Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122, doi: 10.1086/172425
  * André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102, doi: 10.1051/0004-6361/201014666
  * Bae et al. (2014) Bae, J., Hartmann, L., Zhu, Z., & Nelson, R. P. 2014, ApJ, 795, 61, doi: 10.1088/0004-637X/795/1/61
  * Bell & Lin (1994) Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987, doi: 10.1086/174206
  * Bell et al. (1995) Bell, K. R., Lin, D. N. C., Hartmann, L. W., & Kenyon, S. J. 1995, ApJ, 444, 376, doi: 10.1086/175612
  * Caratti o Garatti et al. (2011) Caratti o Garatti, A., Garcia Lopez, R., Scholz, A., et al. 2011, A&A, 526, L1, doi: 10.1051/0004-6361/201016146
  * Contreras Peña et al. (2019) Contreras Peña, C., Naylor, T., & Morrell, S. 2019, MNRAS, 486, 4590, doi: 10.1093/mnras/stz1019
  * Contreras Peña et al. (2017) Contreras Peña, C., Lucas, P. W., Minniti, D., et al. 2017, MNRAS, 465, 3011, doi: 10.1093/mnras/stw2801
  * Dunham et al. (2008) Dunham, M. M., Crapsi, A., Evans, Neal J., I., et al. 2008, ApJS, 179, 249, doi: 10.1086/591085
  * Dunham & Vorobyov (2012) Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52, doi: 10.1088/0004-637X/747/1/52
  * Engelbracht et al. (2007) Engelbracht, C. W., Blaylock, M., Su, K. Y. L., et al. 2007, PASP, 119, 994, doi: 10.1086/521881
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10, doi: 10.1086/422843
  * Fischer et al. (2019) Fischer, W. J., Safron, E., & Megeath, S. T. 2019, ApJ, 872, 183, doi: 10.3847/1538-4357/ab01dc
  * Fischer et al. (2012) Fischer, W. J., Megeath, S. T., Tobin, J. J., et al. 2012, ApJ, 756, 99, doi: 10.1088/0004-637X/756/1/99
  * Fischer et al. (2017) Fischer, W. J., Megeath, S. T., Furlan, E., et al. 2017, ApJ, 840, 69, doi: 10.3847/1538-4357/aa6d69
  * Fischer et al. (2020) —. 2020, ApJ, 905, 119, doi: 10.3847/1538-4357/abc7cb
  * Froebrich & Makin (2016) Froebrich, D., & Makin, S. V. 2016, MNRAS, 462, 1444, doi: 10.1093/mnras/stw1766
  * Furlan et al. (2016) Furlan, E., Fischer, W. J., Ali, B., et al. 2016, ApJS, 224, 5, doi: 10.3847/0067-0049/224/1/5
  * Green et al. (2016) Green, J. D., Kraus, A. L., Rizzuto, A. C., et al. 2016, ApJ, 830, 29, doi: 10.3847/0004-637X/830/1/29
  * Grosso et al. (2020) Grosso, N., Hamaguchi, K., Principe, D. A., & Kastner, J. H. 2020, A&A, 638, L4, doi: 10.1051/0004-6361/202038185
  * Gutermuth et al. (2009) Gutermuth, R. A., Megeath, S. T., Myers, P. C., et al. 2009, ApJS, 184, 18, doi: 10.1088/0067-0049/184/1/18
  * Gutermuth et al. (2008) Gutermuth, R. A., Myers, P. C., Megeath, S. T., et al. 2008, ApJ, 674, 336, doi: 10.1086/524722
  * Habel et al. (2021) Habel, N. M., Megeath, S. T., Booker, J. J., et al. 2021, ApJ, 911, 153, doi: 10.3847/1538-4357/abded8
  * Harper et al. (2018) Harper, D. A., Runyan, M. C., Dowell, C. D., et al. 2018, Journal of Astronomical Instrumentation, 7, 1840008, doi: 10.1142/S2251171718400081
  * Hillenbrand & Findeisen (2015) Hillenbrand, L. A., & Findeisen, K. P. 2015, ApJ, 808, 68, doi: 10.1088/0004-637X/808/1/68
  * Hsieh et al. (2019) Hsieh, T.-H., Murillo, N. M., Belloche, A., et al. 2019, ApJ, 884, 149, doi: 10.3847/1538-4357/ab425a
  * Ioannidis & Froebrich (2012) Ioannidis, G., & Froebrich, D. 2012, MNRAS, 425, 1380, doi: 10.1111/j.1365-2966.2012.21556.x
  * Kenyon et al. (1993) Kenyon, S. J., Calvet, N., & Hartmann, L. 1993, ApJ, 414, 676, doi: 10.1086/173114
  * Kenyon et al. (1994) Kenyon, S. J., Gomez, M., Marzke, R. O., & Hartmann, L. 1994, AJ, 108, 251, doi: 10.1086/117064
  * Kryukova et al. (2012) Kryukova, E., Megeath, S. T., Gutermuth, R. A., et al. 2012, AJ, 144, 31, doi: 10.1088/0004-6256/144/2/31
  * Landsman (1993) Landsman, W. B. 1993, in Astronomical Society of the Pacific Conference Series, Vol. 52, Astronomical Data Analysis Software and Systems II, ed. R. J. Hanisch, R. J. V. Brissenden, & J. Barnes, 246
  * Lang (2014) Lang, D. 2014, AJ, 147, 108, doi: 10.1088/0004-6256/147/5/108
  * Lee et al. (2021) Lee, Y.-H., Johnstone, D., Lee, J.-E., et al. 2021, ApJ, 920, 119, doi: 10.3847/1538-4357/ac1679
  * Mainzer et al. (2011) Mainzer, A., Bauer, J., Grav, T., et al. 2011, ApJ, 731, 53, doi: 10.1088/0004-637X/731/1/53
  * Megeath et al. (2012) Megeath, S. T., Gutermuth, R., Muzerolle, J., et al. 2012, AJ, 144, 192, doi: 10.1088/0004-6256/144/6/192
  * Megeath et al. (2016) —. 2016, AJ, 151, 5, doi: 10.3847/0004-6256/151/1/5
  * Meisner et al. (2017) Meisner, A. M., Lang, D., & Schlegel, D. J. 2017, AJ, 154, 161, doi: 10.3847/1538-3881/aa894e
  * Morales-Calderón et al. (2011) Morales-Calderón, M., Stauffer, J. R., Hillenbrand , L. A., et al. 2011, ApJ, 733, 50, doi: 10.1088/0004-637X/733/1/50
  * Muzerolle et al. (2013) Muzerolle, J., Furlan, E., Flaherty, K., Balog, Z., & Gutermuth, R. 2013, Nature, 493, 378, doi: 10.1038/nature11746
  * Offner & McKee (2011) Offner, S. S. R., & McKee, C. F. 2011, ApJ, 736, 53, doi: 10.1088/0004-637X/736/1/53
  * Park et al. (2021) Park, W., Lee, J.-E., Contreras Peña, C., et al. 2021, ApJ, 920, 132, doi: 10.3847/1538-4357/ac1745
  * Pilbratt et al. (2010) Pilbratt, G. L., Riedinger, J. R., Passvogel, T., et al. 2010, A&A, 518, L1, doi: 10.1051/0004-6361/201014759
  * Plunkett et al. (2015) Plunkett, A. L., Arce, H. G., Mardones, D., et al. 2015, Nature, 527, 70, doi: 10.1038/nature15702
  * Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2, doi: 10.1051/0004-6361/201014535
  * Reach et al. (2005) Reach, W. T., Megeath, S. T., Cohen, M., et al. 2005, PASP, 117, 978, doi: 10.1086/432670
  * Rebull et al. (2014) Rebull, L. M., Cody, A. M., Covey, K. R., et al. 2014, AJ, 148, 92, doi: 10.1088/0004-6256/148/5/92
  * Reipurth & Aspin (2004) Reipurth, B., & Aspin, C. 2004, ApJ, 608, L65, doi: 10.1086/422250
  * Riaz et al. (2021) Riaz, R., Schleicher, D. R. G., Vanaverbeke, S., & Klessen, R. S. 2021, MNRAS, 507, 6061, doi: 10.1093/mnras/stab2489
  * Rieke et al. (2004) Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25, doi: 10.1086/422717
  * Safron et al. (2015) Safron, E. J., Fischer, W. J., Megeath, S. T., et al. 2015, ApJ, 800, L5, doi: 10.1088/2041-8205/800/1/L5
  * Scholz et al. (2013) Scholz, A., Froebrich, D., & Wood, K. 2013, MNRAS, 430, 2910, doi: 10.1093/mnras/stt091
  * Sharma et al. (2020) Sharma, R., Tobin, J. J., Sheehan, P. D., et al. 2020, ApJ, 904, 78, doi: 10.3847/1538-4357/abbdf4
  * Sheehan et al. (2020) Sheehan, P. D., Tobin, J. J., Federman, S., Megeath, S. T., & Looney, L. W. 2020, ApJ, 902, 141, doi: 10.3847/1538-4357/abbad5
  * Stutz & Kainulainen (2015) Stutz, A. M., & Kainulainen, J. 2015, A&A, 577, L6, doi: 10.1051/0004-6361/201526243
  * Stutz et al. (2013) Stutz, A. M., Tobin, J. J., Stanke, T., et al. 2013, ApJ, 767, 36, doi: 10.1088/0004-637X/767/1/36
  * Tobin et al. (2020) Tobin, J. J., Sheehan, P. D., Megeath, S. T., et al. 2020, ApJ, 890, 130, doi: 10.3847/1538-4357/ab6f64
  * Tobin et al. (2021) Tobin, J. J., Offner, S. S. R., Kratter, K. M., et al. 2021, arXiv e-prints, arXiv:2111.05801.  https://arxiv.org/abs/2111.05801
  * Visser et al. (2015) Visser, R., Bergin, E. A., & Jørgensen, J. K. 2015, A&A, 577, A102, doi: 10.1051/0004-6361/201425365
  * Vorobyov & Basu (2015) Vorobyov, E. I., & Basu, S. 2015, ApJ, 805, 115, doi: 10.1088/0004-637X/805/2/115
  * Werner et al. (2004) Werner, M. W., Roellig, T. L., Low, F. J., et al. 2004, ApJS, 154, 1, doi: 10.1086/422992
  * Whitney et al. (2003) Whitney, B. A., Wood, K., Bjorkman, J. E., & Cohen, M. 2003, ApJ, 598, 1079, doi: 10.1086/379068
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868, doi: 10.1088/0004-6256/140/6/1868
  * Zhu et al. (2009) Zhu, Z., Hartmann, L., Gammie, C., & McKinney, J. C. 2009, ApJ, 701, 620, doi: 10.1088/0004-637X/701/1/620