The Power of SOFIA/FORCAST in Estimating Internal Luminosities of Low Mass Class 0/I Protostars

With the Stratospheric Observatory for Infrared Astronomy (SOFIA) routinely operating science flights, we demonstrate that observations with the Faint Object infraRed CAmera for the SOFIA Telescope (FORCAST) can provide reliable estimates of the internal luminosities, $L_{\rm int}$, of protostars. We have developed a technique to estimate $L_{\rm int}$ using a pair of FORCAST filters: one"short-wavelength"filter centered within 19.7-25.3 $\mu$m, and one"long-wavelength"filter within 31.5-37.1 $\mu$m. These $L_{\rm int}$ estimates are reliable to within 30-40% for 67% of protostars and to within a factor of 2.3-2.6 for 99% of protostars. The filter pair comprised of F25.3$\mu$m and F37.1$\mu$m achieves the best sensitivity and most constrained results. We evaluate several assumptions that could lead to systematic uncertainties. The OH5 dust opacity matches observational constraints for protostellar environments best, though not perfectly; we find that any improved dust model will have a small impact of 5-10% on the $L_{\rm int}$ estimates. For protostellar envelopes, the TSC84 model yields masses that are twice those of the Ulrich model, but we conclude this mass difference does not significantly impact results at the mid-infrared wavelengths probed by FORCAST. Thus, FORCAST is a powerful instrument for luminosity studies targeting newly discovered protostars or suspected protostars lacking detections longward of 24 $\mu$m. Furthermore, with its dynamic range and greater angular resolution, FORCAST may be used to characterize protostars that were either saturated or merged with other sources in previous surveys using the Spitzer Space Telescope or Herschel Space Observatory.


INTRODUCTION
The Spitzer Space Telescope enabled large infrared surveys of nearby star-forming molecular clouds yielding a census of young stellar objects (YSOs) in each cloud. In particular, two Spitzer legacy projects, "From Molecular Cores to Planet-Forming Disks" (c2d; Evans et al. 2003) and "Gould's Belt" (GB), observed star-forming regions in 18 molecular clouds, resulting in the identification of 2966 YSO candidates, including 326 protostellar (Class 0/I) candidates (Dunham et al. 2015). Two other Spitzer legacy projects were focused on the large star-forming regions of the Taurus (Rebull et al. 2010) and Orion ) molecular clouds, within which more than 3800 YSO candidates, including at least 500 protostellar candidates, were identified.
Among the most straightforward observational characteristics of protostars to derive is the bolometric luminosity, provided the spectral energy distributions (SEDs) are sufficiently covered, especially in the farinfrared and submillimeter regimes that dominate the emission. However, many protostars have not been observed at these wavelengths and, if they have, the observations may lack the angular resolution necessary to reliably characterize the thermal emission from dust in the protostellar envelope. Furthermore, the bolometric luminosity is "contaminated" by external heating by the interstellar radiation field; the internal (photospheric and accretion) luminosity, L int , better represents an intrinsic property of the protostar. Differences between bolometric and intrinsic luminosities tend not to be significant for typical or high luminosity protostars; those with luminosities < ∼ 1.0 L are most affected by external heating (e.g., Whitney et al. 2013;Dunham et al. 2008;Evans et al. 2001). Dunham et al. (2008) found that fluxes at 70 µm alone were reliable indicators of L int .
Spitzer and Herschel surveys provided 70 µm fluxes for protostars, which may be used to estimate their internal luminosities. However, many protostars either lack 70 µm observations, or these observations suffer from insufficient dynamic range or angular resolution. With Spitzer and Herschel no longer obtaining such observations, a different approach is necessary to derive these estimates. We therefore use radiative transfer models to investigate, in a manner similar to that of Dunham et al. (2008), the relationships between internal luminosities and FORCAST mid-infrared fluxes, which provide better dynamic range and angular resolution. We demonstrate that FORCAST observations are sufficient to estimate internal luminosities of protostars with reliability comparable to that achieved by 70 µm observations. In §2, we summarize the protostar models used in this study. We discuss in §3 the relevant characteristics of FORCAST imaging observations adopted to survey these models. We present in §4 results from these models, which confirm consistency with previous studies; we characterize relationships between observed FORCAST fluxes and internal luminosities of protostars. In §5, we discuss the applicability and limitations of our results, and how these results may be used to further investigate low-mass protostars in nearby star-forming environments. We summarize our findings in §6.
Following Dunham et al. (2008), who used the RadMC code (Dullemond & Dominik 2004) to model protostars observed with Spitzer IRAC (3-8 µm; Fazio et al. 2004) and MIPS (24, 70 µm;Rieke et al. 2004), we considered 350 models of typical protostars and flared disks within rotationally flattened protostellar envelopes, heated by external interstellar radiation fields (ISRFs), with assumed properties as summarized in this section. For each model, we obtained results for 10 inclinations, i, uniformly spaced between cos i of 0 (edge-on disk) and 1 (face-on disk), or cos i = [0.05, 0.15, 0.25, ..., 0.95]; thus, 3500 SEDs were constructed with a distribution of inclinations reflecting that expected for real protostars randomly oriented. To limit statistical variations in the emergent fluxes, each model followed 10, 40, or 160 million photons, whichever was sufficient to yield signal-to-noise ratios (SNRs) of at least 5 at all inclinations and wavebands considered in this study, where SNRs were computed by Hochunk3d following Wood et al. (1996).
The protostars emit as blackbodies at temperature 3000 K with randomly selected (uniformly, in log space) luminosities in the range 0.03-30 L , extending to more luminous protostars than Dunham et al. (2008). As mentioned in Crapsi et al. (2008), the precise temperature assumed for the protostars is not critical since all of the emission is reprocessed by the disks and envelopes.
The flared protostellar disks have a density structure, ρ disk , that decreases as a power law in the midplane radially ( ) while decreasing exponentially perpendicular to the midplane (z) according to (e.g., Shakura & Sunyaev 1973;Lazareff et al. 1990;Pringle 1981;Bjorkman 1997;Hartmann 1998;Whitney et al. 2003b): (1) where R * is radius of the protostar, and the scale height increases as a power law, h( ) ∝ ( /R * ) β with β = 9/7, which is consistent with a self-irradiated passive disk (Chiang & Goldreich 1997). No accretion energy is considered in the disks. The disks have inner radii given by the sublimation temperature and outer radii of 100 AU, where the scale height is 20 AU. The disk masses, which set the overall density normalizations, ρ d0 , are randomly selected (uniformly, in log space) in the range 10 −5 − 10 −3 M .
The rotationally flattened envelopes have density profiles, ρ env , that may be parameterized in terms of the centrifugal radius, R c , and the polar angle, θ 0 , of the streamline of infalling material at large radial distances, r (Cassen & Moosman 1981;Ulrich 1976): where µ ≡ cos θ and µ 0 ≡ cos θ 0 . The constant ρ o is defined by whereṀ env is the mass infall rate in the envelope, M * is the mass of the central protostar, M disk is the mass of the disk, and M disk M * . The Ulrich profile assumes the gas is in free-fall towards a fixed central mass. While the Ulrich profile is typically adopted for the entire envelope, as we have also done in our study, we remind readers that it most accurately reflects free-fall envelope densities at radial distances, r, within which the mass is dominated by the central protostar rather than the disk or envelope. Thus, the Ulrich profile deviates from an accurate collapse profile as the envelope mass interior to r increases appreciably relative to M * (e.g., Shu 1977), which is likely the case in real protostellar envelopes (see §5).
In this formulation, the three input parameterṡ M env , M * , and R c suffice to specify the envelope density profile. The parametersṀ env and M * are related to the density normalization and collapse timescale. The parameter R c is related to rotation and is often set equal to the disk radius. A fourth (less important) parameter arises because of the necessity to set a maximum cloud envelope radius, R env , in order to compute a model. Other formulations of the Ulrich profile are present in the literature; for example, Furlan et al. (2016) preferred to express the profile in terms of the envelope density at a fiducial radial distance (1000 AU), assuming M * = 0.5 M . We instead have recast Equation 2 in terms of M env , which is often used in the literature, and use it to set the normalization instead ofṀ env , by noting that the streamlines become radial at large r, so that in the limit r/R c → ∞ and µ/µ 0 → 1 the envelope density simplifies to become The mass of the envelope is then given by In order to facilitate comparison with the Dunham et al. (2008) and Crapsi et al. (2008) results, we adopt M * = 0.5 M for the model suite. Similarly, we consider envelopes with outer radii of 14,000 AU, envelope masses randomly selected (uniformly, in log space) in the range 1-10 M , R c randomly selected in the range 100-900 AU, and bipolar cavities (created from protostellar outflows) with shape following the streamline with opening angle of 15 • ; the density within each cavity is set to the density of the outermost region of the envelope (Dunham et al. 2008). As discussed in Whitney et al. (2013), the external ISRF adopted by default in Hochunk3d is that found by Mathis et al. (1983) for the solar neighborhood, while Dunham et al. (2008) adopted that of Black (1994) modified at ultraviolet wavelengths for consistency with Draine (1978). Evans et al. (2001) discuss differences between these ISRFs, though we note that the default ISRF in Hochunk3d does not include the cosmic background component dominating at millimeter wavelengths. For consistency with Dunham et al. (2008), we adapted Hochunk3d to use the "Black-Draine" ISRF in our current study. To account for environmental differences among protostellar envelopes, including differing amounts of dust in the molecular clouds surrounding these envelopes, the strength of the ISRF was adjusted by a scale factor and then attenuated and reddened. For each envelope, this scale factor is ran-domly selected in the range 1/3 − 3, distributed logarithmically about unity, and the dust visual extinction is randomly selected in the range 1-5 magnitudes.

Dust Grain Properties
The optical properties of the envelope dust adopted by Dunham et al. (2008) were not available; therefore, we experimented with different dust grain populations available in the literature. The first three grain populations that we considered were readily available in the Hochunk3d distribution. The first population, which we refer to as "KMH-ice" dust, was that found by Kim et al. (1994) for the average Galactic interstellar medium, except with water-ice mantles making up the outer 5% (in radius) of the grains. The second population that we tried was the "molecular cloud model" (hereafter, referred to as "MCM") dust appropriate for protostellar envelopes and described in Whitney et al. (2013). The third dust population was "model 1" dust, which we refer to as "WM1" dust, used by Wood et al. (2002) to model the disk of the classical T Tauri star HH 30 IRS. A fourth population was the thinly ice-mantled, coagulated dust of Ossenkopf & Henning (1994), often referred to as "OH5" grains in the literature (e.g., Evans et al. 2001;Shirley et al. 2005), augmented by the opacities of Pollack et al. (1994) at wavelengths shorter than 1.25 µm, as described in Dunham et al. (2010). The last population was that of Ormel et al. (2011) adopted by Furlan et al. (2016), which includes a mixture of ice-coated silicate and bare graphite grains of radii 0.1-3 µm. The OH5 and Ormel populations were not available in Hochunk3d, but we included them for this study. For reference, the opacities and albedos for the five considered grain populations are shown in Figure 1.
Observationally derived, infrared and submillimeter dust opacities, for protostellar environments, relative to the opacity at 2.2 µm are shown in Figure 2. A comparison with the relative opacities from grain populations considered in this study suggests that the OH5 grains best reproduce these observations. For this reason, we adopt the OH5 population in this study.
As evident in Figure 2, none of the grain populations yield relative opacities at 1.2-850 µm that are fully consistent with observations. Increasing the relative opacity of OH5 grains by 35% for λ ≥ 2.5 µm yields better agreement with observations. In order to obtain some handle on how a grain population better constructed for protostellar environments may affect our results, we rerun the models with these "revised OH5" opacities. We stress, however, that artificially increasing the midinfrared and submillimeter relative opacities of the OH5 grains is not consistent with element abundance constraints of grain populations; such opacities would result from larger grains, and inclusion of these larger grains would necessarily come at the expense of smaller grains to conserve element abundances. Constructing a pro- tostellar grain population is beyond the scope of this study.

FORCAST FILTERS AND SENSITIVITIES
The FORCAST instrument (Herter et al. 2012;Adams et al. 2010) on SOFIA (Young et al. 2012) obtains mid-infrared images and spectra at 5.4-37.1 µm on two detectors: the short-wavelength channel (SWC) and the long-wavelength channel (LWC). Using a dichroic, these channels simultaneously image two wavebands; alternatively, a single channel may be used to directly image one waveband. For our study, we consider only FORCAST images using the seven filters, listed in Table 1, in the range 19.7-37.1 µm for typical observing conditions: specifically, an altitude of 41,000 feet, 7.1 µm of precipitable water vapor at the zenith, and telescope pointings at 50 • from the zenith (e.g., Horn & Becklin 2001).
For purposes of discussion, we adopt fiducial sensitivity limits as those point source flux densities associated with SNR=3 after an hour exposure time. Most FOR-CAST surveys of star-forming regions are likely to require greater SNRs achieved in reasonable times; thus, we expect most studies will focus on sources brighter than given by these limits. Using the online SOFIA Instrument Time Estimator 1 , we determined these fiducial sensitivity limits, in typical observing conditions, for 1 https://dcs.sofia.usra.edu/proposalDevelopment/SITE/ the FORCAST filters operating in direct and dichroic modes, as listed in Table 1. Figure 3 illustrates the effective transmissions of four of these FORCAST filters, accounting for the atmosphere, optics (e.g., the filter itself, optical blockers, and dichroic, as appropriate), and detector response. Except for the F 24.2 µm filter, which operates only in direct mode, we included in Hochunk3d the dichroic transmission functions of the filters listed in Table 1 in order to derive FORCAST flux densities of protostellar models. For the F 24.2 µm filter, we included the direct transmission function. For each filter, the shapes of the direct and dichroic transmission functions are similar; the primary difference is in the overall scale factor of Figure 3. The effective transmission functions of four of the seven FORCAST filters considered in this study, assuming typical observing conditions. The functions associated with observations obtained in direct and dichroic modes are plotted by the black and red curves, respectively, and account for absorption by the atmosphere and optical elements as well as the detector response. the transmission. Therefore, no significant difference is expected in flux densities derived from dichroic and direct transmission functions, only in the observing time required to detect them, particularly for filters with effective wavelengths greater than 30 µm.

RESULTS
Following the approach of Dunham et al. (2008) and adopting OH5 dust, as discussed in §2.1, our results for FORCAST and MIPS fluxes, at a distance of 140 pc, as a function of L int , are shown in Figure 4, demonstrating that L int is best indicated by the 70 µm flux. This figure illustrates increased scatter, particularly at smaller wavelengths and for lower luminosity protostars. The increased scatter is primarily due to geometric effects of inclination. Scatter introduced from inclination may be understood by referring to the SEDs of a standard Class I protostar shown in Figure 14 of Whitney et al. (2013). The flux at 70 µm is relatively unchanged with inclination, while fluxes at shorter wavelengths, particularly for λ < 40 µm, vary considerably for the same protostar observed at different inclinations.
Like Dunham et al. (2008), we derive fluxes within 6 -radius (840 AU at 140 pc) apertures, which is the Spitzer resolution (full width at half maximum; FWHM) at 24 µm, for λ < 40 µm and within 20 -radius (2800 AU at 140 pc) apertures at 70 µm. Typical resolutions achieved by FORCAST filters in 19.7-37.1 µm are 2.1-3.4 ; thus, in principle, our aperture fluxes for FOR-CAST filters will capture a greater fraction of the total fluxes. At least for point sources, the 6 -radius aperture already captured most of the flux in Spitzer observations; any difference between 6 -radius aperture fluxes derived from Spitzer observations and those derived from FORCAST observations is expected to be negligible.
Using a linear least-squares fitting method, we determined the best-fit parameter values, m and b, characterizing the dependence of flux, F ν (measured in erg s −1 cm −2 ) at a distance 140 pc, on L int (measured in L ): where we note that observed photometry is typically given in terms of flux density, S ν = F ν ν −1 . The best fits are illustrated in Figure 4, and the associated parameter values are listed in Table 2, which also lists the standard reduced chi-squared χ 2 red statistics for assessing the quality of these fits. More directly meaningful is Column 5 of Table 2, which lists values for the dispersion (σ) between best-fit and input log L int values where the best-fit values can be explicitly written, for clarity, as using the values for m and b listed in Table 2. These dispersions, σ L , provide a direct means for quantifying the reliability of L int estimates based on these fits. For example, σ L = 0.12 when using 70 µm fluxes; thus, L int estimates based on these fluxes are reliable to within a factor of 1.3 for 67% of the models (i.e., 1σ) and 2.3 for 99% of the models (i.e., 3σ), assuming normal distributions. In contrast, the 19.7 µm fluxes, which yield σ L = 0.55, result in luminosities reliable only to within a factor of 45 (3σ; factor of 3.5 for 1σ). Clearly, the capability that Spitzer and Herschel had in obtaining 70 µm fluxes was critical in characterizing protostars. Equation 10, with fluxes at an adopted distance of 140 pc, may be converted to a form directly applicable to observations, for which S ν is typically given in Jy and valid for any distance d, as where ν is the effective frequency, given in Hz, of the filter, and the best-fit parameter values m and b may be obtained from Table 2. Focusing on 70 µm, for example, L int may be estimated using Thus, a protostar observed at 70 µm to be 1 Jy at a distance of 140 pc suggests that L int is ∼0.1 L , reliable to within a factor of 2.3 (3σ), as previously discussed.
With the scatter in the correlations between FOR-CAST fluxes and L int being primarily a function of inclination, we explored whether utilizing two FORCAST fluxes may improve estimates of L int . Again referring to Figure 14 of Whitney et al. (2013), the slopes of the SEDs in the 20-40 µm regime appear to be correlated with inclination, suggesting that two FORCAST fluxes would in principle provide a first-order luminosity estimate from the average flux level and second-order correction to the estimate from the slope. For example, Kryukova et al. (2012) found that, for protostars lacking 70 µm fluxes, better luminosity estimates could be achieved by considering both the Spitzer 24 µm fluxes and slopes of available 3.6-24 µm SEDs than by considering only 24 µm fluxes alone. Such luminosity estimates were reliable to within a factor of ∼11 (3σ), compared to a factor of 48 (3σ; from σ L = 0.56 in Table 2) based on 24 µm fluxes alone, representing a marked improvement.
Our approach to use two FORCAST fluxes is similar to that by Kryukova et al. (2012) to use Spitzer 3.6-24 µm SEDs, but we might expect a greater improvement since FORCAST extends to longer mid-infrared wavelengths. Toward this end, we considered pairs of FORCAST filters, where the first filter was one of longer wavelengths (i.e., 31.5-37.1 µm) and the second filter  Table 1. The black and red points represent our models with cos i ≥ 0.5 and cos i < 0.5, respectively. Those models dominated by photons originating from the ISRF (see §5) are identified by plus signs; FORCAST is not sensitive to such models. Our best-fit lines to all models, excluding the ISRF-dominated models, are shown as black lines. The MIPS 24 µm and 70 µm panels include the fits from Dunham et al. (2008) as green lines, for reference.
was one of shorter wavelengths (i.e., 19.7-25.3 µm). Linear regression was then used to determine the best-fit coefficients to where the fluxes at 140 pc associated with Filter 1 and Filter 2 are denoted as F ν1 and F ν2 , respectively. Table 3 lists these coefficients for the different filter pairs, and Figure 5 compares L int (fit) with those input into the model. In general, there is reasonable agreement for all models, particularly those detectable by FORCAST, with increased dispersion for intrinsically fainter protostars. Similar to Table 2, Table 3 also lists values for σ L to quantify the reliability of luminosity estimates based on these fits for all models detectable by FORCAST.
Regardless of the FORCAST filter combination, the two-filter fits provide luminosity estimates reliable at least to within a factor of 2.6 (3σ). Filter combinations including F 25.3 µm generally provide the most constrained estimates. The best FORCAST filter combination is F 37.1 µm with F 25.3 µm, which yields σ L = 0.12, providing luminosities reliable to within a factor of 2.3 (3σ), comparable to that achieved from 70 µm observations.
Our fits to Equation 13 may be recast in a form more directly applicable to observations of a source at distance d, in general, as  where S ν1 and S ν2 are the observed flux densities in Jy in Filters 1 and 2, respectively, and Λ is the coefficient accounting for the conversion of units and overall normalization given by where ν 1 and ν 2 are the effective frequencies, in Hz, associated with Filters 1 and 2, respectively. For example, focusing explicitly on 37.1 µm and 25.3 µm, L int may be estimated using where the flux densities S ν,37.1 and S ν,25.3 are given in Jy.
While FORCAST filter pairs yield L int estimates with reliability comparable to those achieved previously with 70 µm observations, observational biases are evident in Figure 5 and depend on the specific filter pair, protostellar luminosity, and sensitivity of the FORCAST observations. Brighter protostars are preferentially detected, resulting in fitted luminosity estimates that systematically overestimate L int , especially evident for low luminosity protostars L int < 0.3 L observed with a filter pair that includes F 19.7 µm. The filter combination F 37.1 µm with F 25.3 µm shows the least bias, though the luminosities are still overestimated for the lowest luminosity protostars. The degree to which L int estimates are biased increase for relatively low luminosity protostars and for less sensitive observations. Figure 5 also shows that nearly all FORCASTdetectable models lie within the 3σ L ranges, extrapolated from computed σ L dispersions listed in Table 3, suggesting that these ranges overestimate the ranges associated with 99% of FORCAST-detectable models. For example, a careful analysis that accounts for the asymmetric and non-normal distribution suggests that Figure 5. Comparison of Lint estimates, derived from mid-infrared fluxes, with input Lint for half of the two-filter combinations considered in this study. Plots associated with Lint estimates using long-wavelength filters F 31.5 µm, F 34.8 µm, and F 37.1 µm are plotted in the left, center, and right panels; those estimates using short-wavelength filters F 19.7 µm and F 25.3 µm are plotted in the top and bottom panels. These plots are similar to those utilizing filters F 24.2 µm and F 33.6 µm, which are not included in this figure. The black points represent models with a flux in at least one of the two relevant FORCAST bands less than the dichroic fiducial sensitivity limit; the red points represent models detectable within the fiducial 1-hour exposures. ISRF-dominated models are not included. The solid lines represent perfect agreement between the estimates and model input values, while the dashed lines represent estimates within 3σL ranges.
L int (fit) is consistent with L int to within a factor of 1.9-2.1 for 99% of models detectable by F 37.1 µm and F 25.3 µm, slightly smaller than the factor of 2.3 extrapolated from σ L . (A difference in σ L of only 0.01-0.02 accounts for this effect.) In other words, σ L slightly understates the reliability of L int estimates derived from FORCAST filter pairs. Given different systematic uncertainties, such as discussed in §5.2 and §5.3, we continue to adopt σ L from Table 3 as they are conservative measures of the reliability of L int estimates. Figures 4 and 5, FORCAST is not sensitive to the faintest of protostars. While protostars at 140 pc with L int > ∼ 0.2 L are detectable at 37.1 µm, only those more luminous than ∼0.7 L are detectable at 19.7 µm. Our method to use a pair of FORCAST filters to determine internal luminosities of protostars is viable for protostars detectable in those filters, which is driven primarily by the sensitivity at shorter wavelengths. Thus, the choice of filter pair is important. Furthermore, the viability of this method is not solely a function of the internal luminosity, but inclination and other properties play a role.

DISCUSSION As illustrated in
In Figure 6, we explore the interplay of internal luminosity, envelope mass, and inclination in determining the detectability of protostars with two different filter pairs: F 37.1 µm with F 19.7 µm; and F 37.1 µm with Figure 6. Comparison of (Lint, Menv) parameter space probed by our models and those detectable by FORCAST, for select pairs of filters and inclinations. Top panels are models observed with filters F 37.1 µm and F 19.7 µm; bottom panels are models observed with filters F 37.1 µm and F 25.3 µm. Left panels include models for all inclinations, while the right panels exclude models with nearly edge-on inclinations (specifically, with cos i < 0.2). First, models that are undetectable by FORCAST are plotted with larger and thicker black symbols, and then models that are detectable by FORCAST are overplotted as smaller and thinner red symbols. The symbol is a dot for a model dominated by protostellar radiation; it is a plus sign for a model dominated by scattered or reprocessed radiation from the ISRF. Note that there are no red plus signs plotted since all models dominated by ISRF are undetectable by FORCAST and therefore appear as black plus signs. Since there are 10 models (one for each inclination) for each (Lint, Menv) probed, only red dots appear for those models detectable for all inclinations; red dots on top of larger black dots appear for models detectable for only some inclinations; and only larger black dots appear for models undetectable for all inclinations. Finally, the dotted and solid curves correspond to the region where 50% and 25%, respectively, of the models are detectable by FORCAST. F 25.3 µm. Comparing the two left panels, we see that protostars of lower luminosities are detectable when using F 25.3 µm rather than F 19.7 µm, as expected, especially for greater envelope masses. While envelope mass affects detectability, it has less impact (i.e., the solid and dotted curves exhibit steeper slopes) when using F 25.3 µm. Comparing each of the right panels with its adjacent left panel, we see that a greater fraction of lower luminosity protostars are detectable, if those models with nearly edge-on protostars are excluded.
The sensitivity of F 25.3 µm, over that of F 19.7 µm, to more models is a compelling reason to favor it. As previously mentioned, Table 3 demonstrates that F 25.3 µm paired with F 37.1 µm yields L int estimates with less uncertainty. Given these considerations of completeness and precision, observations of protostars with F 37.1 µm and F 25.3 µm are likely best for the purpose of determining L int .
We find models for which more than half of the radiation within 20 -radius (2800 AU at 140 pc) apertures are reprocessed or scattered photons originating from the external ISRF rather than from the protostellar system. While such observed ISRF photons are dependent on the strength of field, extinction from the parental molecular cloud, and properties of the protostellar envelope, they are not tied to the internal protostellar luminosity. Thus, contribution (or "contamination") from the ISRF primarily serves to add scatter in the relationships between F ν and L int , and it enables a guide to the level of precision possible on estimates of L int for protostellar systems in a typical range of environments. ISRF-dominated models are more prevalent at 19.7 µm than at 25.3 µm and are found more for nearly edge-on, less luminous protostars. While Hochunk3d enables tracking of the sources of the photons imaged from each system, an observer, in general, does not know a priori the relative contribution of the ISRF. It may be possible to estimate this contamination based on the observed radial profile, enabling results with better L int precision. For our study, we did not pursue such an investigation since we were able to estimate L int with comparable precision to that provided by Spitzer and Herschel. Furthermore, models dominated by the ISRF in the mid-infrared FORCAST bands are at least two orders of magnitude below the fiducial FORCAST sensitivity limits.

Inclination and External Heating
Traditionally, the luminosity of a protostar has been determined by integrating an SED, which requires sufficient spectral coverage, especially from the infrared to submillimeter regimes. We refer to this luminosity as the observed bolometric luminosity, L bol (SED). The evacuated cavity and protostellar disk primarily, and envelope density profile secondarily, result in a non-uniform escape of infrared photons. Light detected from a poleon protostar suffers less extinction relative to the same protostar observed edge-on. Therefore, L bol (SED) will overestimate the true bolometric luminosity in the case of the pole-on protostar and underestimate it for the edge-on protostar. This effect, which we refer to as the "flashlight effect," has been documented in the literature (e.g., Zhang et al. 2013;Whitney et al. 2003b;Yorke & Bodenheimer 1999). For example, in their Figure 10, Whitney et al. (2003b) demonstrated that L bol (SED) overestimated the true bolometric luminosity by about a factor of 2 for their pole-on protostars with bipolar cavities and underestimated it by 50% for the same protostars observed edge-on. Our models show a similar trend.
The bolometric luminosity includes the internal luminosity of the protostar as well as a component, or "contamination," due to external heating by the ISRF. For our models, this contamination is typically ∼0.3 L , consistent with previous studies (e.g., Evans et al. 2001), and reaches as high as ∼1 L in some cases. Thus, not only does L bol (SED) suffer from the flashlight effect, but the contamination from external heating can be significant, particularly for protostars with L int < ∼ 1 L . In our method to estimate protostellar luminosities from FORCAST fluxes, the fluxes were empirically fit to L int (Equation 13); thus, it calibrates out the effect of external heating, in a statistical sense. But, does our method suffer from the flashlight effect? Utilizing a pair of FORCAST filters was intended to account for inclination, the primary factor in the large scatter in correlations between FORCAST fluxes and L int shown in Figure 4. In Figure 7, we plot L int relative to L int (fit) derived from F 37.1 µm and F 25.3 µm, as a function of inclination, demonstrating that our method results in reliable luminosity estimates that do not depend on inclination. Plots for other filter pairs show similar results. Thus, our method successfully uses pairs of FORCAST filters to estimate L int to better characterize protostars more efficiently than obtaining a full SED to determine L bol (SED).

Consideration of Aperture Sizes
Deriving FORCAST fluxes from 6 -radius apertures enabled us to compare our results directly with Dunham et al. 2008, who used the same aperture size for 10-40 µm. In principle, the flux-luminosity relationships derived by our study and by Dunham et al. 2008 apply strictly to fluxes derived with the same physical size of the aperture -i.e., radius of 840 AU. In practice, however, because these apertures include most (typically > ∼ 90%) of the mid-infrared emission from the protostars, it is not important to adhere to the same physical aperture size when deriving fluxes for protostars at different distances, if the apertures include a larger physical area. For example, one could simply use 6 -radius apertures to measure fluxes for all protostars in the Gould's Belt molecular clouds, with distances ∼140-500 pc. The relatively small amount of flux added by including a region of 3000 AU for protostars at 500 pc compared to a region of 840 AU for those at 140 pc potentially introduces a systematic error that is insignificant compared to other dominant systematic errors. Furthermore, since our relationships involve pairs of FORCAST filters, such a systematic error is expected to be even more muted since the measured mid-infrared FORCAST fluxes will be increased similarly in both filters when increasing the physical aperture size. Figure 7. Plot of Lint/Lint(fit) as a function of cos i for models detectable by FORCAST filters F 37.1 µm and F 25.3 µm, demonstrating that our method using this filter pair yields Lint estimates that do not depend on inclination. Individual models are represented by red points (the same models plotted by red points in the bottom right panel of Figure 5), while the median and dispersion (computed in log space) at each inclination are represented by the black filled circles and error bars. The horizontal solid line represents the median value for all detectable models, and dashed lines represent the 3σL range. Plots for other filter pairs show similar results.
We tested these expectations explicitly by using total flux densities (obtained from large 100 -radius apertures capturing emission from the entire protostellar envelope of outer radius 14,000 AU) in the relationships given by Equation 14 to derive L int estimates that were typically within 5% of those obtained from the 6 -radius apertures from which the relationships were derived. We do not recommend using Equation 14 (or Equation 16) with flux densities obtained from apertures much smaller than 840 AU since the more extended size of the protostar (from scattered light) at shorter FORCAST wavelengths relative to that at longer wavelengths may result in greater systematic errors.

Impact of Dust Grain Population
While OH5 grain opacities are most consistent with observational constraints, there are discrepancies, as discussed in §2.1. To quantify the impact of the shortcoming of OH5 grains on L int estimates derived from a pair of FORCAST filters, we reran our protostellar models using "revised OH5" grains, with 35% greater opacity for λ ≥ 2.5 µm compared to OH5 grains. The flux densities at 37.1 µm and 25.3 µm for all FORCASTdetectable models were used to obtained L int estimates using Equation 16. These estimates were typically ∼5% less than those obtained using OH5 grains; most mod-els with revised OH5 grains yielded L int estimates that were within 5-10% of those obtained with OH5 grains. We therefore expect any improved dust model to have a relatively small effect on our results.

Applicability of our Results
We stress that our results apply to those embedded protostellar sources at an early evolutionary stage exhibiting a protostellar envelope, as described in §2. Historically, such protostars have been observationally identified as Class 0 or Class I (hereafter, Class 0/I) sources, based on thermal dust emission or the slopes, α, of the infrared SEDs (Lada 1987, Andre et al. 1993. Class 0/I sources are commonly defined as those with α ≥ 0.3, while Flat sources are those with 0.3 > α ≥ −0.3 and Class II sources, primarily representing evolved YSOs, are those with −0.3 > α ≥ −1.6 (Greene et al. 1994).
While sources observationally identified as Class 0/I are likely bona fide protostars with envelopes, it is possible that some of these sources instead are the more evolved Class II sources obscured by sufficient molecular cloud material such that their SEDs mimic those of protostars. Such "contamination" is most prevalent in embedded young clusters, where the intracluster material may provide significant extinction along the lines of sight. Based on a Spitzer study of the NGC 2264 and IC 348 clusters (Forbrich et al. 2010), up to a third of sources previously identified as Class 0/I were found to be consistent with extincted Class II sources, though half of these possible extincted Class II sources are also still consistent with being Class 0/I protostars. More recently, Carney et al. (2016) use HCO + J=3-2, C 18 O J=3-2, and 850 µm observations to distinguish Class 0/I protostars from extincted Class II sources in Perseus and Taurus; they found ∼30% of sources classified as Class 0/I based on their SEDs were likely extincted Class II sources. Thus, most sources classified as Class 0/I based on their infrared SEDs are bona fide Class 0/I protostars, especially in relatively isolated star-forming regions, with no more than ∼20-30% expected to be extincted Class II sources in regions of embedded clusters. If a source exhibits α ≥ 0.3, it is most likely a protostar, and our results most likely apply.
A source exhibiting α < 0.3 may also be a protostar described by our modeling; in fact, the distribution of α exhibited by our models, shown in Figure 8, peaks at α ∼ 0.1 and shows an extended tail for α ≤ −0.3, slopes that are more indicative of Flat-spectrum and Class II sources. Since most decreasing infrared SEDs are associated with Class II sources, additional evidence beyond the 2-24 µm SED would be necessary to believe reasonably that a particular source with such an SED is a protostar, rather than an evolved Class II source.
The case for identifying a source as a protostar may be bolstered by FORCAST 19-37 µm observations. Figure 8 shows that all protostellar models exhibit midinfrared SED slopes, α(19-37 µm) > 0.5, in the FOR- Figure 8. Slopes of the SEDs characteristic of our protostellar models. The mid-infrared slopes α (19-37 µm), derived from FORCAST bands, are plotted relative to the infrared slopes α traditionally defined by the 2-24 µm bands. Those models detectable by FORCAST are plotted as red points, while those models not detectable are plotted as black points. The histograms along the left represent the distributions of α (19-37 µm), with black showing the distribution for all models and red showing that for only the detectable models. The histograms along the top represent the distributions of α, with the same color scheme. Focusing on the FORCASTdetectable models, the median slopes are α ∼ 0.3 and α(19-37 µm) ∼ 3.5, with 95% of these models falling into the following ranges: α = [−0.3, 0.2] and α(19-37 µm) = [1, 6.5].
The identification and classification of protostars is beyond the scope of this paper, but we have provided some considerations based on the traditionally defined α and how FORCAST observations to derive α(19-37 µm) may help. If the sources are indeed protostars, those same FORCAST observations can be used to estimate their internal luminosities by Equation 14 (or, Equation 16 specifically for 25.3 µm and 37.1 µm observations).

Envelope Mass
The assumptions and applicability of the Ulrich envelope density profile are important to consider, especially in the context of high envelope masses relevant for protostar models. The assumptions include: 1) free-fall toward a central mass, M * +M disk , at the free-fall velocity given by setting r = R env to focus on effects near the adopted cloud boundary and 2) pressure terms that are small compared to kinetic energy, which can be written as the where a s is the thermal sound speed of the gas. For example, the fiducial case of R env = 14, 000 AU and M * + M disk ≈ M * = 0.5 M , as assumed in our modeling, results in v ff = 0.25 km s −1 , which is on order of the thermal sound speed of a s = 0.19 km s −1 for T = 10 K gas (e.g., Terebey et al. 1984). Thus, pressure terms are important near the adopted edge of the envelope, with the result that the density distributions in real protostellar envelopes will deviate from that given by the Ulrich profile near the envelope boundaries. The first assumption, as expressed in Equation 17, applies if the central mass dominates the gravitational potential. However, the envelope masses considered here and in previous studies (e.g., Dunham et al. 2008), which also use the Ulrich envelope model, typically exceed the central masses. To evaluate the size of the effect, we compare the Ulrich envelope mass computed from Equation 7, which further assumes the central mass is dominated by the star (i.e., M disk M * ), with that of the TSC84 (Terebey et al. 1984) self-consistent cloud collapse model. This comparison is appropriate because the TSC84 model includes pressure effects and asymptotically matches the Ulrich model at small radii.
In the TSC84 model, outside the collapsing region of the envelope and when rotational effects are small, a simple formula (involving the leading term) gives the total mass interior to r (e.g., Equation 3 of TSC84; Shu 1977): where r > R exp (t), the radius of the expansion wave representing the boundary of the collapsing region at time t and given by R exp (t) = a s t. In terms of the mass infall rate, given byṀ where m 0 = 0.975 (Shu 1977), this total mass may be expressed as The difference between this total mass and the central mass, which has already collapsed to form the protostar and disk, is the desired envelope mass interior to r > R exp in the TSC84 model: Using Equation 21 to substitute for M tot , and noting that the central mass is the envelope mass interior to r > R exp is then given by The factor in parenthesis has a value of order unity at r = R exp and of order 2 at r R exp . Thus, the factor in parenthesis ranges from about 1 to 2 over the valid r ≥ R exp regime. Note that computing the envelope mass at smaller r requires using the full collapse solution, and can be done numerically, but it is not necessary for our purpose.
For direct comparison with the envelope mass in the Ulrich free-fall model, we rewrite Equation 7 in terms of v ff , using Equation 17, to obtain which becomes the inequality, if indeed the pressure terms are small compared to kinetic energy (i.e., Equation 18). This expression is similar in form to that of the TSC84 model in Equation 24 above and demonstrates that the Ulrich profile underestimates the envelope mass, compared with a realistic envelope model having both gravity and pressure terms. For the adopted case of M * = 0.5 M , Figure 9 compares the envelope mass at different accretion rates and at fixed r = R env = 14, 000 AU envelope radius. The mass difference is about a factor of two, large enough to be important at millimeter wavelengths where the observations are sensitive to cloud (i.e. envelope) mass. However, the mass difference is less important at the mid-infrared wavelengths relevant to this study, where the fluxes generated are not sensitive to the treatment of the outer boundary (e.g. Whitney & Hartmann 1993).
6. SUMMARY In this study, we have established an approach whereby a pair of FORCAST filters may be used to estimate the luminosities of protostars. Empirical relationships are derived for different combinations of a long-wavelength filter (F 31.5 µm, F 33.6 µm, F 34.8 µm, F 37.1 µm) paired with a short-wavelength filter (F 19.7 µm, F 24.2 µm, F 25.3 µm). We find that the best pairing is F 37.1 µm with F 25.3 µm, resulting in luminosity estimates reliable to within a factor of 2.3 for 99% of protostars, which is comparable to the precision achievable in previous studies utilizing Spitzer or Herschel 70-µm data. The luminosity is estimated using Equation 16 once the flux densities S ν,37.1 and S ν,25.3 in Jy are known. Table 3 gives results for other FORCAST filter pairs, which may be used with Equations 14 and 15 to estimate luminosities.
With many protostars lacking data at wavelengths 70 µm or longer, obtaining FORCAST observations and applying our results may be the best approach currently to determine their luminosities. Furthermore, the higher angular resolution achievable by FORCAST enables partitioning of emission among sources blended in previous observations and better constraints on the SEDs and luminosities of the components. Our approach requires data using only a pair of FORCAST filters, not a well covered SED in the infrared and submillimeter regimes, and is independent of the inclination of the protostar.
In §2.1 we consider available dust model opacities. Figure 2 shows that the OH5 opacity (augmented with Pollack optical constants) fits the observational data best, though not perfectly. We find that an improved dust model would affect the luminosity estimates by only 5-10%, thus supporting the choice of OH5 dust for protostellar envelopes.
We also compare ( §5.5) the commonly assumed Ulrich density profile for protostellar envelopes with a more realistic profile for envelope masses comparable to, or greater than, the embedded source. Real protostellar envelopes likely have material collapsing slower than assumed free fall velocities, and pressure terms become appreciable in the outer regions of the envelopes. Such considerations suggest that the Ulrich profile underestimates total envelope masses by about a factor of two, but this deficiency has little effect on observed infrared emission from the protostar.