Far-Infrared Luminosity Bursts Trace Mass Accretion onto Protostars

Evidence abounds that young stellar objects undergo luminous bursts of intense accretion that are short compared to the time it takes to form a star. It remains unclear how much these events contribute to the main-sequence masses of the stars. We demonstrate the power of time-series far-infrared (far-IR) photometry to answer this question compared to similar observations at shorter and longer wavelengths. We start with model spectral energy distributions that have been fit to 86 Class 0 protostars in the Orion molecular clouds. The protostars sample a broad range of envelope densities, cavity geometries, and viewing angles. We then increase the luminosity of each model by factors of 10, 50, and 100 and assess how these luminosity increases manifest in the form of flux increases over wavelength ranges of interest. We find that the fractional change in the far-IR luminosity during a burst more closely traces the change in the accretion rate than photometric diagnostics at mid-infrared and submillimeter wavelengths. We also show that observations at far-IR and longer wavelengths reliably track accretion changes without confusion from large, variable circumstellar and interstellar extinction that plague studies at shorter wavelengths. We close by discussing the ability of a proposed far-IR surveyor for the 2030s to enable improvements in our understanding of the role of accretion bursts in mass assembly.


INTRODUCTION
In young stellar objects (YSOs), stochastic enhancements of the accretion rate from the protoplanetary disk onto the star seem to play an important role in building up the stellar mass (Kenyon et al. 1990;Offner & Mc-Kee 2011;Dunham & Vorobyov 2012;Fischer et al. 2019;Lee et al. 2021;Zakri et al. 2022;Wang et al. 2023).To some degree, these events supplement the consistent but slowly declining flow of disk mass that persists for a few million years.The stochastic enhancements, known as bursts or outbursts, are observed as dramatic increases in the bolometric luminosities of the objects.
Historically, these bursts have been divided into two categories.The FU Orionis events (FUors) feature increases of 4 to 5 mag at optical wavelengths over a few months, and this state persists for many years.The first known example began its burst in 1936 (Herbig 1966) and still has not subsided.The EX Lupi events (sometimes referred to as EXors) feature increases of about 1 mag at optical wavelengths, persist for months to a few years, and may recur in any given star (Herbig 2007).In optical and near-infrared (near-IR) spectroscopy, EX Lup stars resemble young stellar objects with atypically high accretion rates.FU Ori stars, on the other hand, have very different spectra that are dominated by a self-luminous accretion disk.See Reipurth & Aspin (2010), Audard et al. (2014), Section 2.5 of Hartmann et al. (2016), andFischer et al. (2023) for reviews of the burst phenomenon.
Descriptions of these two classes are necessarily somewhat vague due to the diversity of light curves and spectra seen within each.In recent years, two developments have complicated the picture further.First, many outbursts with intermediate photometric and spectroscopic characteristics have been observed.V1647 Ori has been considered emblematic of these (Contreras Peña et al. 2017).It began an outburst in 2003, faded, and then recovered to its outburst state in 2008.Spectroscopically, it resembles an FU Ori star except for H I and CO emission that are more typical of EX Lup stars (Connelley & Reipurth 2018).
Second, infrared surveys of star-forming regions or the Galactic plane (Scholz et al. 2013;Caratti o Garatti et al. 2017;Contreras Peña et al. 2017;Lucas et al. 2017;Fischer et al. 2019;Lucas et al. 2020;Park et al. 2021;Zakri et al. 2022) and targeted submillimeter observations (Herczeg et al. 2017;Lee et al. 2021) have detected bursts from objects with SEDs consistent with protostars that are still deeply embedded in their natal envelopes.Mid-IR searches frequently define bursts as brightenings in excess of 1 or 2 mag, finding that they occur in a given protostar on ∼ 1000 yr timescales (Park et al. 2021;Zakri et al. 2022).These are less dramatic and more frequent than the classical FU Ori bursts, which brighten in the optical by ∼ 5 mag and occur on ∼ 10,000 yr timescales (Hillenbrand & Findeisen 2015).Therefore it is crucial to continue investigating these protostellar bursts to constrain the responsible physics.Safron et al. (2015) reported the first known burst of a Class 0 protostar, in which the infalling envelope is inferred to still contain more mass than the star itself.The burst, of HOPS 383, was discovered in an analysis of 24 µm photometry.This protostar brightened in 2004-2006 and had faded by 2017 (Grosso et al. 2020), a light curve suggestive of V1647 Ori.Due to its embedded nature, no optical/near-IR spectroscopy was possible to determine whether it had features similar to EX Lup, FU Ori, or something else.
The case of HOPS 383 is instructive because it tests the limits of what we can determine without spectroscopy and with only sporadic photometry at the necessary wavelengths.While the onset and end of the burst could be dated with reasonable precision, it was difficult to monitor the evolution of the light curve at intermediate times, and it was impossible to grasp the physical insight into the nature of the burst that would have been possible with optical and near-IR spectroscopy.
Due to these observational limitations, it is at present difficult to answer a fundamental question about the assembly of low-mass stars: What fraction of their mainsequence masses is assembled in stochastic bursts?To understand this, we need to determine, at least in a statistical sense, the amplitudes and durations of bursts.Recently, Wang et al. (2023) presented a framework for using burst durations and amplitudes to infer what fraction of a star's mass is accreted during each mode and applied it to 70 years of photometry of EX Lup.They concluded that the largest bursts are responsible for about half of the accreted mass during the monitoring time, with the other half arising from the combination of small bursts and quiescent accretion.
It is particularly important to measure the burst durations and amplitudes for Class 0 protostars, before the majority of the stellar mass has been assembled.This paper demonstrates that a far-infrared (far-IR) survey of Galactic star-forming regions is essential for determining these quantities for statistically significant numbers of Class 0 protostars.Protostellar SEDs for a wide range of evolutionary states and viewing angles peak in the far-IR (Whitney et al. 2003), so they are most easily detected at such wavelengths.
For the protostellar mass assembly question specifically, far-IR monitoring has two additional benefits.First, this wavelength range gives us the most reliable estimates of the true burst amplitudes.In the far IR, the spectral energy distribution (SED) is most directly responsive to changes in the bolometric luminosity and therefore the accretion rate.Second, the far IR allows us to circumvent the high circumstellar and interstellar extinction that cause incompleteness in mid-infrared (mid-IR) surveys and that can itself be time dependent, producing spurious burst-like signals.
In this work, we start with 1.2-870 µm SEDs that were fit by Furlan et al. (2016) to 92 Class 0 protostars, the largest population of this class in a single molecular cloud complex within 500 pc, the Orion Molecular Clouds.This allows us to sample the range of protostellar luminosities and evolutionary states present in a real cloud complex.For each SED, we increase the luminosity in the associated model by fixed factors and consider how effectively those factors are recovered by photometry in different wavelength regimes.In a companion paper, Lee et al. (in preparation) consider the observational cadence needed for a far-IR survey to quantify burst rates, amplitudes, and durations.
In Section 2, we introduce the set of model SEDs used to fit the Class 0 Orion protostars and how these models are modified to account for bursts.In Section 3 we examine how burst amplitudes depend on wavelength.In Section 4 we show the benefits of the far IR for circumventing foreground extinction.In Section 5 we discuss how future far-IR missions such as the proposed Probe Far-Infrared Mission for Astrophysics (PRIMA) will be useful in addressing the issues presented here.Section 6 presents our conclusions.Here we describe how model SEDs are used to infer the effects of bursts on photometry of Class 0 protostars at different wavelengths.We start with a set of model SEDs used to fit the multiwavelength SEDs of Class 0 Orion protostars by Furlan et al. (2016) in the context of the Herschel Orion Protostar Survey (HOPS), a key program of the Herschel Space Observatory.Led by PI S. T. Megeath, HOPS is described by Fischer et al. (2020).We then modify these SEDs to account for bursts of varying intensity and examine how these bursts are reflected in photometry at different wavelengths.

The HOPS Model Grid
Furlan et al. ( 2016) presented a grid of 30,400 model SEDs designed to infer physical properties of protostars.The grid contains 3040 model protostars, and each has ten SEDs as viewed from different angles.These model SEDs were then fit to the observed SEDs of 319 Orion protostars, constructed with data from the Two-Micron All-Sky Survey (2MASS; Skrutskie et al. 2006), the Spitzer Space Telescope (Werner et al. 2004), the Herschel Space Observatory (Pilbratt et al. 2010), and the Atacama Pathfinder Experiment (APEX; Siringo et al. 2009Siringo et al. , 2010)).Table 1 gives the wavelengths and approximate angular resolutions of datasets used for the SEDs.The techniques used to obtain the photometry are explained in Section 3.1 of Furlan et al. (2016).
Of the 319 protostars, Furlan et al. (2016) determined 92 to be of Class 0, the least evolved class, via their mid-IR spectral slopes and their bolometric temperatures.In Section 3, we analyze potential bursts from 86 of these Class 0 protostars, dropping six sources where the model There is a central star, not visible at this scale, and a flared disk (salmon).The radius of the disk can vary.For visibility, here it is shown at twice the size of the largest model.Additionally there is an envelope (light blue) of fixed radius 10,000 au.The envelope density profile follows the Terebey et al. (1984) rotating collapse prescription and is normalized to one of 19 values.The envelope cavity has an opening angle that ranges from 5 • to 45 • in steps of 10 • .All five possibilities are shown with dotted lines, and the 15 • case is shown as evacuated.Finally, the system can be viewed over a range of angles from nearly pole-on (18 • ) to nearly edge-on (87 • ).
SED fits either lack sufficient mid-IR flux for the analysis or are already near the maximum luminosity in the grid.
The grid is described in full detail by Furlan et al. (2016); here we summarize the important properties.Models are calculated with the 2008 version of the HOCHUNK radiative transfer code (Whitney et al. 2003).Each model contains a central protostar, a circumstellar disk, and an envelope with radius 10,000 au.The envelope density is that of a rotating, collapsing core with a fixed infall rate (Terebey et al. 1984), plus a bipolar cavity evacuated by an outflow.Figure 1 illustrates the axisymmetric geometry of the model.
The disk and envelope consist of dust with opacities from Ormel et al. (2011) that account for grain growth at an age of 3 × 10 5 yr and include the effects of ices.Properties relevant to embedded protostars are varied, as described in the next two paragraphs, while others are held constant.Constant parameters of note include the stellar mass (0.5 M ⊙ ), stellar temperature (4000 K), disk mass (0.05 M ⊙ ), and envelope radius (10,000 au); see Table 3 of Furlan et al. (2016) for the full list.This is an important distinction from larger, more widely used grids that vary stellar properties and other features that have little effect on SEDs when there is a dense envelope present.
Four parameters are varied in the grid.First, the total luminosity of the protostar can take on one of eight values extending from 0.1 L ⊙ to 303 L ⊙ in roughly equal logarithmic steps.Because the stellar and accretion luminosity are both reprocessed by the disk and envelope, protostellar SEDs depend only weakly on the relative contributions of these inputs.Therefore the luminosity is modified by adjusting the radius of the star or the accretion rate onto the star.To improve the fits by allowing for a continuous range of luminosities, a limitedrange scaling factor is applied, based on the finding that the shape of an SED is only weakly dependent on small changes in the luminosity (Kenyon et al. 1993).Table 2 shows how the required luminosities are obtained for different choices of stellar radius, accretion rate, and scaling factor.
Second, the envelope density can take one of 19 values.The envelope density at 1000 au can be 0 or can range from 1.19 × 10 −20 g cm −3 to 1.78 × 10 −16 g cm −3 in roughly equal logarithmic steps.The envelope density elsewhere is set by the rotating collapse model mentioned above.Third, the outer radius of the disk, which is also the centrifugal radius of the envelope, can be 5, 50, 100, or 500 au.Fourth, the cavity opening angle can be 5, 15, 25, 35, or 45 • .
To account for inclination (i) dependence, each model can be viewed from one of ten angles ranging from 18 • (almost face-on) to 87 • (almost edge-on) in equal steps of cos i. Model fluxes are determined for 24 apertures that run from 420 au to 10,080 au in steps of 420 au.For each wavelength, an aperture is chosen that represents the angular resolution of the relevant instrument (Table 1), assuming a source distance of 420 au.Finally, foreground extinction is applied with laws from Mathis (1990) if A J < 0.76 or McClure (2009) otherwise.For each protostar, the best model is the one that has the least average weighted logarithmic deviation between the observed and model SED (Fischer et al. 2012).

Adding Bursts to the HOPS Model Grid
The best-fit model for each protostar has a luminosity associated with it; see Column 8 of Table 1 in Furlan et al. (2016).To simulate bursts, we increase each of these luminosities by factors of 10, 50, and 100 and choose a model from the existing grid with this lumi-  2, the given stellar radius and mass accretion rate are chosen, and a scaling factor within the given range is chosen to reproduce the burst luminosity.As with the initial fitting exercise, we again assume that the SED does not depend significantly on the division between stellar and accretion luminosity or on luminosity changes covered by the range of scaling factors.These assumptions were validated by Furlan et al. (2016).
Figure 2 shows the original, factor-of-10 burst, and factor-of-100 burst models for the nearly face-on protostar HOPS 256 and the nearly edge-on protostar HOPS 397.The upper panels show how the SEDs change, brightening and shifting to bluer wavelengths as the bursts intensify.The lower panels show how the ratios of the two burst SEDs to their quiescent counterpart vary with wavelength in different ways depending on the details of the best-fit model.

BURST AMPLITUDES WITH WAVELENGTH
For each model, we calculate flux densities (for single wavelengths) or integrated luminosities (for a range of adjacent, simultaneously observed wavelengths) that would be observed by several existing and planned telescopes.To simulate WISE (Wright et al. 2010) or NE-OWISE (Mainzer et al. 2014) observations, we calculate the 3.4 and 4.6 µm flux densities.To simulate observations such as those obtained by the James Clerk Maxwell Telescope transient survey (Herczeg et al. 2017), we cal-  2016) with a nearly face-on viewing angle of 18 • , and HOPS 397 (right) was modeled by the same authors with a nearly edge-on viewing angle of 87 • .In each upper panel, black symbols mark photometry (+ signs) or limits (triangles) from Spitzer, Herschel, and APEX.The gray spectrum for HOPS 256 is from Spitzer.The red curve is the best-fit SED from a radiative transfer model, and the blue and magenta curves show the SEDs for the same models, scaled up by factors of 10 or 100 in luminosity, respectively.Rectangles show the wavelength regions of photometry discussed in Section 3: purple for the WISE 3.4 and 4.6 µm bands, green for one of the far-IR space mission concepts (see Section 5), and orange for the JCMT 450 and 850 µm bands.In each lower panel, the ratios of the factor-of-10 and factor-of-100 burst SEDs to the best-fit SEDs are plotted with solid blue and magenta lines, while dotted blue and magenta lines mark fixed ratios of 10 and 100.In the mid IR, the observed flux changes are greater than the luminosity changes for the face-on protostar and less than the luminosity changes for the edge-on protostar.In the far IR, the integrated flux change gives a robust estimate of the luminosity change.In the submillimeter, the flux changes are significantly less than the luminosity changes.
culate the 450 and 850 µm flux densities.Finally, we calculate the far-IR luminosity obtained by integrating under the 25 to 235 µm SED as observed with 16 separate bandpasses placed over this range at roughly equal logarithmic intervals.
We define the inferred burst amplitude in each wavelength range as the flux density (in the mid IR and submillimeter) or luminosity (in the far IR) observed in the burst SED to that of the pre-burst, quiescent SED.We define the actual burst amplitude as the change in the true luminosity of the protostar; i.e., a factor of 10, 50, or 100.We include 86 Class 0 protostars in the sample.Four of the initial 92 are excluded because their pre-burst luminosities exceed 202 L ⊙ .This places them in the bottom row of Table 2; in this case, simulating a burst involves multiplying the SED by a constant, which is not instructive.Two additional protostars are excluded because their pre-burst models have no detectable 3.4 µm flux, leaving their post-to-pre-burst ratios poorly constrained.
Figure 3 shows how the ratio of the inferred burst amplitude to the actual burst amplitude depends on wavelength.For example, if the model luminosity increased by a factor of 10 but the 3.4 µm flux increased by a factor of 12, the value on the vertical axis of the figure would be 1.2.
In the mid IR, the ratio covers a broad range and is centered slightly above unity.Mid-IR fluxes are heavily dependent on geometry.In protostars that are viewed through their outflow cavities, a disproportionately large amount of the mid-IR flux is emitted along the line of sight, and the fractional increase in the mid-IR flux is large compared to the fractional increase in luminosity.In protostars that are viewed through their disks, very little flux reaches the line of sight, and the opposite occurs.Furthermore, when the envelope is heated by the burst, the location of the effective surface seen by the observer changes and the local temperature associated with this surface also changes such that the emission can vary considerably.For the majority of the envelope, the mid IR is on the Wien side of the SED, so small changes in optical depth or temperature can lead to large changes in flux density.The potential for large changes makes the mid-IR a reasonable wavelength range in which to search for bursts, even if the change in flux density there is likely to overestimate the change in accretion rate.
At submillimeter wavelengths, the SED responds to the temperature change in the disk and envelope, the magnitude of which is significantly reduced compared to the luminosity change (Johnstone et al. 2013).
The far IR traces the luminosity change best, in that it is centered just below unity and has a small dispersion.This is for two reasons.First, these changes are determined by integrating over a broad range of wavelengths, from 25 to 235 µm, capturing the peak of the SED for a broad range of evolutionary states (Whitney et al. 2003).Second, the far IR probes luminosity in a geometry-independent way, responding to the heating of the inner envelope, which is fairly uniform regardless of geometry and viewing angle.
An additional advantage to observing with a large number of far-IR bandpasses is that flux ratios among these bandpasses change as outbursts propagate through and heat the protostellar envelope.Figure 2 shows how the SED peak shifts to shorter wavelengths as the bursts intensify.While any single bandpass in this range would reflect the actual burst amplitude reasonably well, inte-grating over several bandpasses results in a more robust measurement.
The choice in the far IR of 16 separate bandpasses extending from 25 to 235 µm is motivated by plans for a proposed space mission concept (see Section 5).To estimate the extent to which the result depends on the specific number of bandpasses and their wavelengths, we repeated the exercise for two previous far-IR instruments.The Photodetector Array Camera and Spectrometer (PACS; Poglitsch et al. 2010) aboard the Herschel Space Observatory had three bandpasses centered at 70, 100, and 160 µm.The Far-Infrared Surveyor (FIS; Kawada et al. 2007) aboard the AKARI space mission (Murakami et al. 2007) had four bandpasses centered at 65, 90, 140, and 160 µm.
The far-IR symbol in Figure 3 has a median of 0.82, and the central half of the distribution spans a range of 0.26.The ratios for the Herschel simulation have a median of 0.68, and the central half of the distribution spans a range of 0.31.The ratios for the AKARI simulation have a median of 0.84, and the central half of the distribution spans a range of 0.39.The results are similar for these three collections of far-IR bandpasses, but the actual burst amplitudes are more closely recovered, with a smaller dispersion, when many bandpasses are available across the full far-IR range.2020), who used similar methods except that we present the effect of bursts on SEDs that have been fit to an ensemble of 86 well-characterized Class 0 protostars in a single molecular cloud complex.Similar findings regarding the wavelength dependence of flux changes were reported by Stecklum et al. (2021) in an analysis of the burst of G353.93−0.03, a massive YSO.Our results also build on the mid-IR versus submillimeter multiwavelength analysis of observed protostellar variability by Contreras Peña et al. (2020).That study revealed both a clear correlation in brightness variations between the two widely separated wavelengths and a clear indication that neither wavelength responded in a manner directly proportional to the changing central source luminosity.

AVOIDING EXTINCTION WITH FAR-IR SURVEYS
Another reason to survey for bursts in the far IR is to avoid the complications of large, variable circumstellar and interstellar extinction that affect observations at shorter wavelengths.This causes three distinct problems in assessing the rates, amplitudes, and durations of accretion bursts.
First, it can be ambiguous whether brightness changes in the mid IR are due to accretion or to changing foreground extinction.Redder mid-IR colors during out-burst can be used to rule out a dust clearing event, since this would yield bluer colors (e.g., Zakri et al. 2022), but the interpretation is more difficult in the event of bluer mid-IR colors.The far IR is only mildly sensitive to changes in extinction, making it a more robust probe.
We demonstrate the robustness of far-IR and submillimeter photometry against variable extinction in Fig- ure 4. For each of the 86 Class 0 protostars, we consider the foreground A V modeled by Furlan et al. (2016) and change this by a randomly selected value between −30 and 30 mag, enforcing a minimum A V of 0. The threshold of 30 is based on the maximum change seen in V2492 Cyg, a star with a history of extreme variability due to both accretion and extinction changes (Hillenbrand et al. 2013).The mid-IR fluxes vary by a large range, while the far-IR fluxes vary only slightly, and the variation of the submillimeter fluxes is imperceptible.The distributions are peaked slightly below 1, since increasing A V by up to 30 is always possible, while a decrease in A V is capped by our prohibition against negative A V .
Second, point sources detected in mid-IR imaging surveys are often knots of shock emission (Gutermuth et al. 2009;Koenig & Leisawitz 2014) or scattering surfaces (Yoon et al. 2022) within outflows.These may be several arcseconds from the protostar, while the protostar itself is hidden behind intervening dust.Such locations could be directly illuminated by an accretion burst in the central protostar or indirectly heated through shocks associ-ated with varying conditions within the jet-like outflow.Either of these events may occur significantly after the onset of the burst, complicating efforts to understand a given protostar's accretion history.
Third, many Class 0 protostars are too deeply embedded to be detected before a luminosity outburst.The classification of HOPS 383 as a protostar was tenuous before its outburst (Safron et al. 2015).Of the 92 Class 0 protostars in Orion (Furlan et al. 2016), 16 (17%) were not classified as protostars based on mid-IR Spitzer data alone and required Herschel far-IR photometry to be classified as such (Stutz et al. 2013).In regions that are less well characterized than Orion, far-IR surveys will be needed to recover the full population of young Class 0 protostars.

PROSPECTS FOR FAR-IR MONITORING OF PROTOSTELLAR OUTBURSTS
Far-IR monitoring of enough protostars to understand the role of luminosity bursts in stellar mass assembly will require a new mission with the requisite sensitivity, mapping speed, and wavelength coverage.The As-tro2020 Decadal Report recommended "a far-IR imaging or spectroscopy probe mission," prompting NASA to issue an Announcement of Opportunity for an Astrophysics Probe Explorer in July 2023.One of the mission concepts developed in response to this call is the Probe Far-Infrared Mission for Astrophysics (PRIMA; Glenn 2023), with science instruments well suited for far-IR monitoring of protostellar outbursts.
PRIMA is planned to be a cryogenically cooled far-IR observatory for the 2030s, with a nominal mission lifetime of five years.It will improve mapping speeds by up to four orders of magnitude with respect to its far-IR predecessors and is expected to open wide discovery space.
The two science instruments planned for PRIMA are PRIMAger, an imager (Burgarella et al. 2023), and FIRESS, a spectrometer (Bradford et al. 2023).PRIMAger is a sensitive multi-band spectrophotometric imager.It offers hyperspectral narrow-band imaging (with spectral resolving power R ∼ 10) from 25 to 80 µm and polarimetric capabilities in four broadband filters from 80 to 260 µm.FIRESS offers coverage from 24 to 235 µm, in either low spectral resolution (R ∼ 130) or moderate spectral resolution (≳ a few × 1000).
The PRIMA instruments are well suited to monitoring protostellar SEDs over time and performing followup observations of spectral lines to investigate changes in chemistry and physical properties during accretion events.The wide bandwidth of PRIMAger makes it highly effective for monitoring protostellar SEDs, and surveying can be efficient since protostars are highly clustered in the sky.FIRESS can be used to monitor changing physical conditions via analysis of far-IR spectral lines including rotational transitions of H 2 O, high-J CO, OH, [O i] and [C ii] (e.g., Manoj et al. 2013;Karska et al. 2018).Lee et al. (in preparation) argue that a monitoring survey of 2000 protostars is necessary and sufficient to robustly answer the question of whether protostars gain the majority of their mass through accretion bursts or through steady-state processes.They discuss monitoring 2000 protostellar SEDs with PRIMA, tracking their brightness over the five-year mission, and comparing the fluxes to those obtained twenty years earlier with Herschel.An initial spectral scan of all protostars with FIRESS, followed by triggered FIRESS observations during and after observed bursts, will enable detailed monitoring of the changing physical conditions.These observations are also critical for measuring protostellar burst light curves to constrain both the burst mechanisms and the event durations.

CONCLUSIONS
With model SEDs that were fit to the observed SEDs of 86 Class 0 protostars in the Orion molecular clouds by Furlan et al. (2016), we explored how major luminosity outbursts due to accretion manifest as a function of wavelength.The protostars sample the ranges of envelope density, cavity opening angle, and viewing angle encountered in the largest star-forming region in the neareast 500 pc.
We find that photometry over several distinct bands in the far-IR range, e.g., 25 to 235 µm, is essential for three reasons.First, the change in luminosity seen over this range, measured as a fraction of the total luminosity change, has a small dispersion and is close to unity, unlike that seen at shorter or longer wavelengths.Second, the location of the SED peak shifts to shorter wavelengths during the burst, and this can be tracked with multiple bandpasses over this wavelength range.Third, this wavelength range is less susceptible to the confusing effects of variable circumstellar and interstellar extinction than shorter wavelengths.
Due to these advantages, a far-IR time-domain survey of nearby star-forming regions as outlined in Lee et al. (in preparation) may lead to a breakthrough in our understanding of the importance of outbursts in assembling the main-sequence masses of stars.

Figure 1 .
Figure1.Radiative transfer model geometry.There is a central star, not visible at this scale, and a flared disk (salmon).The radius of the disk can vary.For visibility, here it is shown at twice the size of the largest model.Additionally there is an envelope (light blue) of fixed radius 10,000 au.The envelope density profile follows theTerebey et al. (1984) rotating collapse prescription and is normalized to one of 19 values.The envelope cavity has an opening angle that ranges from 5 • to 45 • in steps of 10 • .All five possibilities are shown with dotted lines, and the 15 • case is shown as evacuated.Finally, the system can be viewed over a range of angles from nearly pole-on (18 • ) to nearly edge-on (87 • ).

Figure 2 .
Figure2.Demonstration of the burst procedure for two Class 0 protostars.HOPS 256 (left) was modeled byFurlan et al. (2016) with a nearly face-on viewing angle of 18 • , and HOPS 397 (right) was modeled by the same authors with a nearly edge-on viewing angle of 87 • .In each upper panel, black symbols mark photometry (+ signs) or limits (triangles) from Spitzer, Herschel, and APEX.The gray spectrum for HOPS 256 is from Spitzer.The red curve is the best-fit SED from a radiative transfer model, and the blue and magenta curves show the SEDs for the same models, scaled up by factors of 10 or 100 in luminosity, respectively.Rectangles show the wavelength regions of photometry discussed in Section 3: purple for the WISE 3.4 and 4.6 µm bands, green for one of the far-IR space mission concepts (see Section 5), and orange for the JCMT 450 and 850 µm bands.In each lower panel, the ratios of the factor-of-10 and factor-of-100 burst SEDs to the best-fit SEDs are plotted with solid blue and magenta lines, while dotted blue and magenta lines mark fixed ratios of 10 and 100.In the mid IR, the observed flux changes are greater than the luminosity changes for the face-on protostar and less than the luminosity changes for the edge-on protostar.In the far IR, the integrated flux change gives a robust estimate of the luminosity change.In the submillimeter, the flux changes are significantly less than the luminosity changes.

Figure 3 .
Figure 3. Ratio of inferred to actual burst amplitude for bursts in 86 Class 0 protostars.The 25-235 µm point refers to the result of integrating the SED over far-IR photometric bands from 25 to 235 µm.In these violin plots, the shape along the vertical axis shows the continuous distribution of points obtained via kernel density estimation from the discrete measurements.Horizontal lines mark the first quartile (bottom), median (central), and third quartile (top) of each distribution.

Figure 4 .
Figure 4. Flux ratios for the ensemble of protostars if the V band extinction AV toward each protostar is changed by a randomly selected value between −30 and 30 mag, enforcing a minimum AV of 0. The submillimeter ratios are uniformly too close to unity to be discerned well in the plot.These theoretical results are similar to those found by MacFarlane et al. (2019) and Baek et al. (2020), who used similar methods except that we present the effect of bursts on SEDs that have been fit to an ensemble of 86 well-characterized Class 0 protostars in a single molecular cloud complex.Similar findings regarding the wavelength dependence of flux changes were reported byStecklum et al. (2021) in an analysis of the burst of G353.93−0.03, a massive YSO.Our results also build on the mid-IR versus submillimeter multiwavelength analysis of observed protostellar variability by ContrerasPeña et al. (2020).That study revealed both a clear correlation in brightness variations between the two widely separated wavelengths and a clear indication that neither wavelength responded in a manner directly proportional to the changing central source luminosity.

Table 1 .
Data Used in Orion SEDs

Table 2 .
Dependence of Luminosity on Model Parameters nosity, holding other parameters (envelope density, disk radius, cavity opening angle, inclination, and foreground extinction) fixed.Within each of the luminosity ranges shown in Table