A study of 90 GHz dust emissivity on molecular cloud and filament scales

Recent observations from the MUSTANG2 instrument on the Green Bank Telescope have revealed evidence of enhanced long-wavelength emission in the dust spectral energy distribution (SED) in the Orion Molecular Cloud (OMC) 2/3 filament on 25"(0.1 pc) scales. Here we present a measurement of the SED on larger spatial scales (map size 0.5-3 degrees or 3-20 pc), at somewhat lower resolution (120", corresponding to 0.25 pc at 400 pc) using data from the Herschel satellite and Atacama Cosmology Telescope (ACT). We then extend the 120"-scale investigation to other regions covered in the Herschel Gould Belt Survey (HGBS) specifically: the dense filaments in the southerly regions of Orion A; Orion B; and Serpens-S. Our dataset in aggregate covers approximately 10 square degrees, with continuum photometry spanning from 160um to 3mm. These OMC 2/3 data display excess emission at 3mm, though less (10.9% excess) than what is seen at higher resolution. Strikingly, we find that the enhancement is present even more strongly in the other filaments we targeted, with an average excess of 42.4% and 30/46 slices showing an inconsistency with the modified blackbody to at least 4{\sigma}. Applying this analysis to the other targeted regions, we lay the groundwork for future high-resolution analyses. Additionally, we also consider a two-component dust model motivated by Planck results and an amorphous grain dust model. While both of these have been proposed to explain deviations in emission from a generic modified blackbody (MBB), we find that they have significant drawbacks, requiring many spectral points or lacking experimental data coverage.

Star forming regions within the Milky Way represent an enormous variation in scale from the size of the clouds (∼ 10 18 m) down to the individual stars formed within (∼ 10 9 m). Giant molecular clouds (GMC) are composed of gas and dust which collapse gravitationally to form denser regions which eventually form individual stars. This process is defined by the hierarchical nature of the structure formation as the clouds form filamentary structures (Hennebelle & Falgarone 2012) during their collapse into dense cores. These filaments, especially the smallest and densest ones (proposed characteristic size ∼0.1 pc, Arzoumanian et al. (2019)), are a particularly active area of investigation (Peretto et al. 2013;André et al. 2014). These studies include cataloguing of individual stars in the process of forming (Fiorellino et al. 2020;Könyves et al. 2020;Polychroni et al. 2013;Stutz et al. 2013), mapping the magnetic field strength and orientation within the clouds (Chen et al. 2019;Fissel et al. 2019), and the measuring of the thermal SEDs (Schnee et al. 2014;Sadavoy et al. 2016;Mason et al. 2020).
We focus on the submillimeter and millimeter wave emission from the clouds and internal filamentary structure. This regime is dominated by thermal emission from dust grains which are generally assumed to emit as a modified blackbody (MBB), given by where T d is the temperature of the dust grain determined by the local radiation environment, and β is the spectral index of the dust, set by the physical properties of the grains such as their composition and size. These MBB models of dust emission agree well with observations at the submillimeter and millimeter wavelengths, with dust temperatures varying across the cloud and the spectral index typically taking on a value of 1.5 < β < 2.5 Sadavoy et al. 2013;Aghanim et al. 2016).
In addition to a generic modified blackbody spectrum, there have been recent studies of more specific models which treat the emission more in-depth than a single power-law spectrum. Models from Draine & Lazarian (1999) and Draine & Hensley (2013) take into consideration the ferromagnetic resonance properties of ironcontaining grains which produced large magnetic dipole cross-sections at frequencies of 50-100 GHz. Another such model is the two-component dust model of Meisner & Finkbeiner (2015) which was designed to describe the emission structure of the galactic cirrus and comprises a cold, shallow β(∼ 1.6) component in combination with a hot, sharp β(∼ 2.7) component to describe the diffuse emission. A final model which seeks to describe the behavior of the dust through a description of the particles is the amorphous grain model of (Paradis et al. 2011). This spectrum relies on model of the emissivity of the dielectric grains over the frequency band and at different temperatures.
In the north of Orion A, we find the high line-mass Integral Shaped Filament (ISF) out of which the Orion Nebula Cluster is forming (e.g. Stutz & Kainulainen (2015); Stutz & Gould (2016)). The northern portion of the ISF is commonly denoted OMC2/3, and was the subject of study at high resolution (Schnee et al. (2014), hereafter S14, Sadavoy et al. (2016), hereafter S16). Due to its high line mass and proximity, the region is bright and observationally accessible. Meanwhile, its status as a high mass star and cluster forming filament, with the accompanying advanced evolutionary stage and elevated UV radiation field, motivate the extension of studies of the thermal dust emission to other potentially less disturbed filaments.
The emission by dust in OMC 2/3 has a rich history of investigation. S14 took the initial measurements from the MUSTANG experiment (3.3 mm) and combined them with observations from the MAMBO (1.2 mm) instrument on the IRAM (Institut de Radioastronomie Millimetrique) 30 m telescope and NH 3 -derived temperature maps to model the spectral index β of the thermal dust emission. This analysis showed high emission in the 3.3 mm data and produced surprisingly shallow SEDs with β ∼ 0.9, which was attributed to the potential growth of mm-sized grains in the regions of interest. S16 expanded on this analysis, including 160-500 µm data from the PACS (Photodetector Array Camera & Spectrometer) and SPIRE (Spectral and Photometric Imaging Receiver) instruments as well as 2 mm data from the GISMO (Goddard-Iram Superconducting 2-Millimeter Observer) instrument to further constrain the dust SED. These more constrained SEDs displayed a β ∼ 1.7 − 1.8, in line with expectations (S16). The shallower β of S14 is then attributed to enhanced emission at 3.3 mm (90 GHz).
The most recent development comes from (Mason et al. 2020), hereafter M20, which investigated the enhanced emission found in S16. M20 includes data from the commissioning of the MUSTANG2 (Multiplexed Squid-TES Array at Ninety Gigahertz 2) instrument as well as the GBT (Green Bank Telescope) Ka band receiver at 31 GHz (∼1 cm) to further extend the low frequency arm of the SED. Additionally, this analysis includes SCUBA-2 (Submillimetre Common-User Bolometer Array 2) observations at 450 and 850 µm which serve to further constrain the low-frequency tail of the SED. Results from M20 indicate that the enhanced emission seen at 3.3 mm is confirmed by the 1 cm data point and that it is inconsistent with both anomalous microwave emission (AME) and spinning dust emission. Such flattening of the SED has also been predicted theoretically from amorphous dust models (Meny et al. 2007;Paradis et al. 2011;Coupeaud et al. 2011;Nashimoto et al. 2020), but these predictions have not yet been confirmed observationally.
Understanding the emission properties of astrophysical dust and its spectrum is important for much of astrophysics and cosmology. Dust emission is often used as a mass tracer and estimator (Eales et al. 2012;Groves et al. 2015;Paradis et al. 2019) and a systematic disagreement with models would significantly affect dust mass estimates. In addition, polarized emission is used to trace the plane-of-the-sky projected magnetic field and, under specific assumptions, estimate the field strengths (Heitsch et al. 2001;Fissel et al. 2019;Crutcher 2012). Pertaining to cosmology, the cosmic microwave background (CMB) is measured through sightlines that are contaminated by the polarized emission of the dusty interstellar medium of the Milky Way. These foregrounds are modelled (Vansyngel et al. 2018;Hensley & Bull 2018) and removed to get at the underlying CMB and its polarization spectrum. These important uses of polarized dust emission and the models upon which they are based must, however, also describe the unpolarized emission, in which we have begun to see a break in the spectral properties between 150 and 90 GHz. To this end, we have expanded the studies of S14, S16, and M20 of OMC 2/3 to additional clouds within the Herschel Gould Belt Survey (HGBS) to garner an understanding of whether this break in the spectral properties is widespread in nearby molecular clouds. In this paper we use the methods of M20 extended to larger angular scales and other regions using data from the Atacama Cosmology Telescope (ACT), Herschel PACS, and Herschel SPIRE.
The structure of the paper is as follows. In Section 2 we review the data used in each analysis mode (25 and 120 ), discuss the regions we have targeted in this work, and lay out the map production pipeline. Section 3 details the method of SED extraction and the models used to fit the data. In Section 4 we compare the results of the pipeline at 25 to the previous analysis as well as discuss potential sources of contamination within the data. We discuss the results of this study and the potential for follow-up in Section 5. Finally, we conclude in Section 6. In this paper, we refer to the difference between a model-predicted brightness as enhanced emission or the elevation of a data point and all error bars represent a 1-σ uncertainty (68% confidence) on the associated quantity.

DATASETS
The data that were used for this study varied between the high-and low-resolution analyses. In the initial high-resolution verification and calibration of the pipeline, we included data from the PACS, SPIRE, SCUBA-2, MAMBO (Max-Planck Millimeter Bolometer Array), GISMO, and MUSTANG2 instruments to cover a spectral range from 160µm to 3mm at a resolution of 25 matching the 350µm SPIRE array. The study was extended to low-resolution with the wealth of survey data provided by the ACT instrument. For the low-resolution analysis we retain the Herschel survey data on the high frequency end and combine it with the ACT full-sky data on the low-frequency side. The enormous footprint of the ACT map affords coverage of almost any HGBS region we aim to target. Within the HGBS we targeted a handful of star forming regions, particularly chosen to span a range in density and star formation activity level (isolated vs clustered) while retaining a molecular hydrogen column density in excess of 1e22 cm −2 (Polychroni et al. 2013;Könyves et al. 2020;Fiorellino et al. 2020). For these low-resolution analyses, we match the ACT 90 GHz instrument which has a beam size of 120 and cover the same spectral range of 160µm to 3 mm. At 90 GHz, our high-resolution data comes from the MUSTANG2 instrument on the Robert C. Byrd Green Bank Telescope (GBT). MUSTANG2 features 215 horncoupled transition edge sensor (TES) bolometers, a bandpass of 75-105 GHz, a field of view of 4.25 , and a beam size of 10 giving excellent resolution and mapping power at 90 GHz. For more information about the MUSTANG2 instrument, see Dicker et al. (2014) and Stanchfield et al. (2016).
In addition to producing the initial OMC 2/3 90 GHz map, this instrument provides an excellent option for high-resolution followup mapping of interesting regions seen in the larger survey analysis. As the excess signal is seen in the low-frequency (90 GHz and longer λ) tail, the availability of a high-resolution and on-sky resource for mapping these regions is important to the understanding of this behavior.

Legacy Datasets
Our data in the 160, 250, and 350 µm bands come from the HGBS observations with the PACS and SPIRE instruments. These bolometer cameras operated on the Herschel satellite (Pilbratt et al. 2010). The 350µm band sets the limiting resolution of 25 to which we match all of the maps. At 450 and 850 µm our data comes from the SCUBA-2 instrument (Holland et al. 2013) which is a dual-band bolometer receiver on the James Clerk Maxwell Telescope. The remaining data at 1.2mm and 2mm come from the GISMO (Staguhn et al. 2006) and MAMBO (Bertoldi et al. 2003) instruments which operated on the IRAM 30m telescope.

ACT
The Advanced ACTPol instrument is a cosmic microwave background (CMB) camera on the Atacama Cosmology Telescope that observes roughly half of the sky, mapping the polarization spectrum of the CMB and measuring the Sunyaev-Zel'dovich Effect (SZE) of galaxy clusters. The experiment features more than 5500 polarized transition edge sensors (TES) in five bands, centered at nominal frequencies of 28, 41, 90, 150, and 230 GHz. Further information about the AdvACT experiment can be found in Henderson et al. (2016). We use the 90-230 GHz data as the low-frequency data has not yet been analyzed, which sets a resolution limit of 120 due to the beam size of the 90 GHz band. For this analysis, we used the data release which included observations through the 2017-2018 season (Naess et al. 2020) which has been combined with the Planck maps (Ade et al. 2014a) to recover large scale flux around bright sources such as OMC 1. The product used here is cutouts of the HGBS regions from the full ACT map. Notably, this 3 mm data product is calibrated using the CMB, rather than the planets used for MUSTANG2 observations, which provides an additional layer of insulation from systematics that might contaminate our conclusions.

SPIRE
Extending the scope of this analysis to 120 we give up access to the 450 and 850 µm spectral coverage provided by the SCUBA-2 instrument due to severe filtering effects. We recover this region of the spectrum through the inclusion of the SPIRE 500 µm HGBS maps, which, with a native resolution of 36.7 , no longer exceed the resolution limit of 25 set by the 350 µm maps in the high-resolution analysis. This coverage at 500 µm provides an important data point in the center of the SEDs between the 350 µm SPIRE and 1.4 mm (220 GHz) ACT data in the transition into the low frequency tail.

Instrumental Bandpasses
To ensure our data on the spectral emission curves are properly located in frequency space, we calculate the response-weighted average frequency of each instrument. For the Herschel instruments, we use the frequencies calculated in M20, and, for the ACT instrument, we obtained the spectral response curves R(ν) and calculated the band-center as The response-weighted band centers can be found in Table 1. These bandpass-corrected frequencies are used in the analyses performed in Section 3.2 and amount to small shifts (typically ∼ 1%) in the spectral location of the data points.

DATA REDUCTION
In this section we detail the methods and tools used extract the SEDs from the original maps. We present the pipeline that we have developed to bring these disparate data products into a usable and common format. The discussion of the pipeline includes also the testing and verification steps taken to ensure that the data are treated properly. Throughout this analysis, all flux data are associated with the errors found in 2, which includes calibration error (Herschel - Bendo et al. (2013), ACT -Naess et al. (2020)), typical noise levels in the maps, and the noise introduced through the map processing steps.

Pipeline
A goal of this analysis was to recreate and verify the data reduction pipeline in python based on the IDLproduced results from M20. As such, we started with a goal of accurately reproducing the results seen in the 25 analysis of the OMC 2/3 region, while creating a scalable framework that could be expanded to include other data, such as the 120 data presented in this paper, with minimal modification. We break the process down into multiple steps, each of which we discuss below.   Table 2. Associated errors for each instrument and, in the case of the specific map noise, each specific band. The map noise levels are the rms noise of the signal levels seen in a roughly 82 kilopixel region in each map multiple beam distances from the filamentary structures as compared to the typical signals at each peak.

Initial maps
The first step was to extract sub-maps around the regions of interest from much larger maps, such as the full-sky ACT maps or large Herschel regions. The rough postage-stamping is achieved through a script which uses the functionality of the astropy.io library to load the ACT maps and another map as a set of reference coordinates, which are then used to create a preliminary, large rectangular cutout of the region.
Once the maps have been postage stamped, the fits files are then scanned using a wrapper function which parses all maps' bunit keywords, compares them to a dictionary of known units, and converts them to the desired output unit (MJy/sr in this analysis). Once all maps have been properly resized and the values correctly mapped in the common and desired units, they are saved in preparation for the next step.

Reprojected maps
The next stage in the pipeline is to convert all maps to the same world coordinate system (WCS), pixelization, and desired region. Our output WCS is defined by first selecting the targeted region in a map in DS9 and determining the size and central pixel coordinates from the initial maps. The original map coordinates are then scaled to the new pixelization and an output shape along with a WCS header is created for this selected region. Once these objects are created, we resample the maps to the new coordinate system using the reproject package which is an extension of the utilities found in the astropy library. Specifically, we use the reproject interp function, which performs a bilinear resampling of the original fits image and WCS onto the newly defined WCS. For these maps, we chose to use a final pixelization of 2 to match that of the original MUSTANG2 map.
We performed extensive verification of these functions in order to verify that the output data is both properly scaled and in the correct locations. The first order test mapped the original data to the same WCS, which showed that the manipulations did not shift the locations of the data or appreciably change the values of the pixels. After verifying this, we next tested the flux conservation of the reproject functions, we expect that where p 0 is the pixel values in the original map and p is the pixel values in the new map. The right hand side is also scaled by the relative pixel sizes to conserve flux. This test was applied to both the entire map and small test regions around bright sources in different maps to ensure that there was no bias affecting larger flux values.
Results of this test showed that the effect was 0.05-0.5% in magnitude and that the exact value depended on the structure and change in pixelization of the map itself, which is included our assumed 5% map error.

Resolution matching
The final stage in the map production pipeline is to match the resolution of the maps to that of the lowest resolution instrument. Consider a data-product map, which consists of a grid of sampled data points and associated values (data) and convolved with the instrument point-source response kernel (B 0 ), creating the initial map (map 0 ). To degrade the resolution of this map to that of an instrument with a worse point source response (B f ), a second convolution with another response kernel (b) must be applied such that For the simplest case of a Gaussian beam, which we assume here, this can be solved for the characteristic convolution function b as a Gaussian smoothing function obeying, Where these σ are the standard deviations of the Gaussian smoothing functions. Once these values have been obtained, we create a set of 2D Gaussian convolution kernels and apply them to the maps using the astropy.convolution library and associated functions. As the standard convolution kernel application sets the values beyond the map boundary to 0 for computational purposes, the maps were significantly oversized during development in order to prevent artefacts from appearing in the slice datasets. Once these convolutions have been applied, the maps are in their final states and are fully prepared for data extraction and SED production. We verified the flux conservation of this method across the entire map as, excluding edge effects, the sum of the pixel values in the map should be static. For the initial 25 analysis, the fractional difference in total pixel values ranged from 10 −8 − 10 −17 , effectively zero, with the largest difference in the most drastically changed maps and no values causing appreciable change in the maps. At 120 , the edge effects are more important due to the larger convolution kernel and the boundary conditions. Because of these effects, we see flux conservation errors of up to a 0.5% level across the entire map, though the level returns to the virtually undetectable of the 25 analysis in the central, targeted region of the map where the SED data are extracted.

Slices
To extract our SED datapoints, we used a slicing and deramping method. The concept behind this method is to remove any linear drifts in the background as well as the mean levels that may enter from map making or diffuse emission not associated with the filament itself. For the 25 analysis of OMC 2/3 (verification of the pipeline), the slices were in predetermined locations from the previous analysis in M20. For all other analyses (120 investigations), the location of the slices was determined by hand and covered the high SNR filamentary structures in each map. These slices can be found for all regions in Figure 1. We perform the slicing and deramping operation as the maps are zero-median and this allows us to subtract a linear drift and remove the background mean level near the clouds. The 500 µm maps were chosen to be the representative maps for SNR cuts as they represented both the center of the spectral band and a middle resolution.
In order to choose slices that minimized the noise contribution of looking at low-signal areas, we first made rough signal-to-noise cuts of the regions in each 500 µm map. To do this, we selected a low emission region far away from the filaments of approximately 10 5 pixels in each map, calculated the mean and variance of the signal levels in these dark regions, and subtracted the mean from the map and divided by RMS noise level. With these new SNR maps, we masked all regions below a cutoff level of 5-σ as being regions considered "low-signal". An example of this procedure applied to the OrionA-S map can be seen in Figure 2. We then imported these maps into DS9 and manually created a series of slices that extended far off of the filament into the voided areas at either end and were oriented to slice through regions of varying signal level. Once created, these slices were converted into pairs of endpoints which were fed into the data extraction algorithm detailed below.

Data extraction
The final step in the data processing pipeline is the extraction of the individual SED datapoints from each map. The data are extracted using an algorithm that removes a linear trend from the data along the slice and finds the location of the highest emission within each region. This process serves to create a slice which has the background level removed as well as any linear drifts that may occur across the map from mapmaking artifacts such as bowling, which were particularly evident in the JCMT and MUSTANG2 data. Once the slice has been detrended, the signal at the peak intensity that remains is the SED datapoint.
The method to detrend the data is as follows. First, we apply the slices to the maps. These are defined by their endpoints, between which we draw a line on the map, creating a masked region. Second, at each end point we create a circular mask of radius 20 in which we calculate the mean value, representing the noisesuppressed value at each end of the slice. Third, we create a linear regression model between the two endpoints of the slice to capture any linear changes across the map as well as the background. Fourth, at each point along the slice we calculate the distance from the endpoint and subtract the linear regression model value. Fifth, we restrict our search to the central 60% of the slice (targeted on the filament itself) and find the local maximum, this is saved as the location chosen for a given slice and frequency. Once this process has been repeated for all frequencies in a given slice, the SED data locations are extracted by taking the median of the lo- cations chosen in the previous steps to minimize bias introduced by the noise inherent to the data. Once the median locations are calculated, the data are extracted for each frequnecy at that location in the slice and the process is repeated for all slices, producing the initial SEDs used in this analysis. This technique is partic-ularly adept at removing large-scale backgrounds such as the CO contamination in the 90 GHz ACT+Planck maps versus the ACT-only maps, returning datapoints that agree to sub-noise and calibration levels. These data are discussed further in Section 4.3.

SED fitting
With the data extracted from the maps, we apply models to determine the presence of excess emission at low frequencies as well as verify the pipeline consistency with the 25 data. The dust is modeled to emit as a Modified Blackbody (MBB), which is characterized by free parameters for the dust temperature T d and spectral index, β (see Equation 1). In order to quantify the level of excess emission at low frequencies, we characterize the behavior of three different fits to the data. They are as follows: • Full fit: We fit the entire SED from 160 µm to 3.3 mm with a modified blackbody and characterize the quality of the fit. This is referred to as the MBBall in the SEDs produced.
• Short λ: We fit the SED from 160 µm to 2 mm, omitting the 90 GHz point, with a modified blackbody and characterize the quality of the fit to this data. This fit is referred to as the MBBno90 fit in the produced SEDs.
• Short λ extrapolated: Not a truly unique fit, but an extrapolation of the previous down to the 3.3 mm data, allowing us to characterize the quality of that fit to the entire dataset. This serves as a check as to whether the 90 GHz data is significantly different from what the rest of the data would predict it to be. This is referred to as the extrapolated MBBno90 fit throughout the text.
We fit the data through the use of the scipy.optimize library function curve fit, which, in combination with the MBB model we provide, performs a least squares fit taking into account the uncertainties in the data provided. The SEDs and both fits for each slice can be found in Figures 3-8. We characterize the quality of the fits through the χ 2 per degree of freedom for each variant, the lists of which can be found in Tables 4 and 5. As discussed in Section 2.3, we use the response-weighted average frequency ν 0 for the purposes of fitting, residual calculations, and reduced χ 2 values. Verification of the pipeline is discussed in Section 4 with a comparison of M20 and this analysis fit values found in Figure 9.

SED Monte-Carlo simulations
The initial SEDs are constructed from the data that were extracted using the methods described in Section 3.1.5. To quantify the impact of calibration and measurement errors (see Table 2), as well as any potential biases inherent in our analysis, we performed a suite of Monte-Carlo simulations. The procedure was as follows. For each slice, the original MBBno90 fit values were extracted at the frequencies corresponding to the instruments' bandpass-averaged frequency. Next, a set of 10 6 Gaussian random values were generated, with a mean value of 1 and a σ draw corresponding to the calibration error for each instrument, with the maximally-pessimistic assumption that each instrument (ACT, SPIRE, and PACS) had fully correlated calibration uncertainty. These draws were then applied as multiplicative factors to the original predicted data points and the resulting SED refit, producing a varying suite of parameters (T, β, and A) as well as new predictions for 90 GHz brightness. A pair of corner plots showing the parameter posteriors is displayed in Figure 10, corresponding to a case with no excess emission (top) and a case with a significant 90 GHz excess (bottom). The parameter posteriors are qualitatively similar in both cases.
This suite of simulations provided two key pieces of information: an assessment of the bias of the initial SED fits as compared to the noisy realizations and an understanding of the range of predicted values associated with each fit. Notably, for 53 slices and 10 6 realizations, we find that the initial fits are remarkably unbiased (∆ 90 < 0.2%, see Tables 4 and 5), with typical deviations from the final mean parameter values shown in Table 3. Finally, we use the distribution of predicted values to characterize the disagreement statistics between the SED predictions and the measured map values, the results of which are discussed in Section 5.

90 GHz Monte-Carlo simulations
In addition to using this procedure to model the effects of high-frequency calibration errors on the SEDs, we also applied this technique to determine the range of possible 90 GHz values. These simulations were constructed in a similar manner -all slices shared the same set of 10 6 draws for the calibration errors. The key difference comes in adding the map noise values to each slice data point. In this case, the slices were separated into groups corresponding to the particular map they were taken from and the "dark noise" levels used to create a per-slice additional noise prescription. The resulting collection of data points represent a large sampling of the range of possible values for the 90 GHz emission   and providing insight into the significance of the excesses seen in the data.

Full data fits (MBBall)
We first attempt a simple modified blackbody fit to cover the entire spectrum of the data, similar to the MBBall performed in the M20 analysis in the OMC 2/3 region. After analyzing the initial OMC 2/3 test region (Figure 3, 120 resolution) which was observed in M20, we looked at the three other regions in Serpens, Orion B, and Orion A-S. In Section 5 we will examine the residuals of this fit more closely region by region. With the entire spectrum modeled, we follow the methods of M20 and consider how a spectrum omitting the 3 mm point fits the data. We first restrict ourselves to a MBB fit which covers the range from 160 µm to 2 mm and consider the quality of this fit to the data. Figure 9. Plots of the calculated parameters T d and β by slice for the M20 and our new pipeline analysis methods for the MBBall fit. Error bars for the red points represent the 1-σ uncertainties extracted from the output covariance matrix of the fitting software. For the blue points, the errors are those reported in M20. Of note is that slice 16 comprises poor data which fails to converge to a solution in the new pipeline, and the results of the temperature fit were not reported in the M20 analysis. For these reasons we leave this slice out of the verification analysis.

Parameter
Mean Deviation (nσ) Temperature (T) 0.04 Spectral Index (β) -0.007 Amplitude at 250µm (A) 0.06 Table 3. The mean deviation of the initial fit parameters across all slices as compared to the 10 6 Monte-Carlo realizations. As is shown here, the initial fits were slightly biased high on the amplitude and temperature and slightly low on the spectral index. To correct for this, the bootstrap realization values and their associated distributions have been used in lieu of the initial fits throughout this analysis. Of note, these biases reported here would correspond to values of ∆T = 0.01K and ∆β = −0.0004 on average across the entire sample.  Table 5).
(Bottom) The SED parameter corner plot for Orion A-S slice 9, which showed > 6σ excess emission (see Table 5). While the SED parameters vary between the two slices, the fitted parameters remain robust in each case.
The results of this model are covered in-depth in Section 5. Restricting the data to higher frequencies improves the quality of the fits, which are shown as the green curves in Figures 3-8. If the emission is elevated as we seem to see, we should expect a trend of increased deviation from the model at 3 mm relative to that seen  Table 4. Reduced χ 2 for the Serpens regions in the two different cases as well as the OMC 2/3 region. The first case is the modified blackbody fit to all data point and the χ 2 n calculated from these data. The second case is the blackbody spectrum fit to the λ ≤ 2mm datapoints and the χ 2 n calculated from the λ ≤ 2mm dataset as well. The number of degrees of freedom are shown next to each dataset in parentheses. We have additionally included the fit parameters from the bootstrap analysis of the data and the excess at 90 GHz compared to the prediction from the model fit using only the higher frequencies (short λ model) is shown as a fractional excess ("Excess %"). Slices with bolded "Excess (σ)" line lie above the 4-σ significance level of excess inconsistent with 0. Slice 7 of Serpens 3 was flagged during data quality control due to contamination of the data from a Galactic HII region.
with the full fit, as well as an increased χ 2 n as compared to both other models. Not only do the 90 GHz data lie above the modeled SED curves, but when the fits exclude the 90 GHz data, there is a systematic and large increase in the amount by which the 90 GHz points lie above the curves, clearly showing an excess emission level compared to expectations. As with the previous analyses, the results of this are presented region by region in Section 5.

Further modeling
We also consider further approaches to fitting these spectra from a "quality of fit" perspective. The first approach is to consider a two-component dust population which models the SED as a combination of two components at different temperatures that emit with different spectral characteristics. We follow Meisner & Finkbeiner (2015), hereafter M15, and consider the χ 2 /n (n), no90 params, excess emission Location Full fit (n=4) Table 5. Reduced χ 2 for the OrionA-S and OrionB regions in the two different cases. The first case is the modified blackbody fit to all data point and the χ 2 n calculated from these data. The second case is the blackbody spectrum fit to the λ ≤ 2mm datapoints and the χ 2 n calculated from the λ ≤ 2mm dataset as well. The number of degrees of freedom are shown next to each dataset in parentheses. We have additionally included the fit parameters from the bootstrap analysis of the data and the excess at 90 GHz compared to the prediction from the model fit using only the higher frequencies (short λ model) is shown as a fractional excess ("Excess %"). Slices with bolded "Excess (σ)" line lie above the 4-σ significance level of excess inconsistent with 0. Note that slice 13 for OrionA was added as a "sanity check" and includes a strong 90 and 150 GHz source that breaks the blackbody curve. It was not included in any calculations due to the known contamination.

model
In this model we have wrapped together the fraction and optical efficiency terms f n and q n as well as the overall proportionality constant into a pair of terms A and B which function as the amplitudes of the two dust populations. With our relatively limited number of data points, we opt to use the best-fit parameters found in M15 for β 1 = 1.67 and T 1 = 9.15 to model a second, cool component of the dust emitting in tandem. We then fit the amplitude of this cool component and a full, second modified blackbody to the SEDs. The results of these fits are discussed in Section 5.6. The second model we consider is that of Paradis et al. (2011), in which the grains are treated as an amorphous solid, with an emissivity that differs from the standard treatment. This SED is modeled as Where they have calculated the greybody term as a function of frequencies and temperatures for a variety of datasets which cover the spectrum that we are modeling within. As was tested in M20, we find that the T 0 = 17.5K model provides a reasonable estimate of the temperatures of the dust we model (∼15K average temperature) and use the opacities that correspond to that dust temperature as that measured constant is the closest fit to our data. That said, we do not constrain the model to assume a 17.5 K dust temperature, but allow that to be fit as well. The results of this model are discussed in Section 5.7.

Verifcation of pipeline at 25
The fitting method that is applied here, as well as the data extraction methods are slightly different from those found in M20, and as such we expect a slight variance in the MBB parameters for identical slices at 25 . In order to verify that our pipeline produces reasonable results, we fit the data for the 24 slices extracted with our pipeline using the datasets at 25 from M20. The resulting parameters are compared with those from M20 to ensure that the two methods are in good agreement. A graphical comparison of the fitted parameters with their associated errors can be seen in Figure 9. When we compare the solutions to the least-squares fit that each algorithm finds, we see that there exists a small, but consistent difference in the fit parameters across the majority of the slices. In the case of the fit temperatures T d in each slice, we find that our pipeline predicts a higher temperature across all slices, with a median difference and error of ∆T d /σ T d = -0.190±0.095. For the spectral indices, we find that this trend is reversed, with our pipeline solving for a slightly lower spectral index across most slices, for median difference and error of ∆β/σ β = 0.13±0.11. We accept these small variances as reasonable and associate them with the well-known degeneracy between β and T (Shetty et al. 2009), the different implementations of the slicing and deramping method, and the least-squares fitting algorithms.

Free-Free considerations
Of the most concern for contamination of the SEDs at 90 GHz is free-free emission, which has a ν −0.1 spectrum. In both S16 and M20, this was treated through the examination of serendipitous VLA X-band data which overlapped with the OMC 2/3 region. Of the 24 slices examined in M20, only one was found to have any nonnegligible flux associated with free-free emission. To ameliorate concerns that this contamination might provide false positives for enhanced emission in our maps, we used the NRAO VLA Sky Survey (NVSS) catalog and postage-stamp database (Condon et al. 1998) to identify objects which were potential sources of contamination near or on the data taken from each slice. This was then cross-referenced with the Wide-field Infrared Survey Explorer (WISE) catalog of HII regions (Anderson et al. 2014). This allowed me to determine if the sources were simply 1.4 GHz with a quickly falling off (synchrotron) or slowly dropping (free-free) spectrum. Of the sources identified as potentially problematic (pixel flux assuming free-free spectrum ≥ 5% of the flux at 90 GHz) only slice 7 of the Serpens 3 region was associated with an HII region (G029.825+02.232) within a distance of 15 . Of note, Slice 13 of the Orion A-S region was also flagged and coincided with an HII region; however, that data was already excised as the ACT 90 and 150 maps had shown a strong source not seen in other frequencies and was originally included solely to test the code. We flagged the Serpens 3 Slice 7 region as contaminated at a > 5σ level and removed it from the calculations.

Spectral line contamination
In addition to the free-free considerations discussed above, we considered the potential for map contamination through molecular emission lines. Of particular interest are the coupling of the 115 GHz CO(J=1-0) line to the Planck 100 GHz maps which are used in the production of the ACT dataset and the HCN emission line at 88.6 GHz which sits within all bandpasses of the ACT 90 GHz maps. To evaluate this systematic concern, we reproduced the analysis with the ACT-only maps, which are not sensitive to the CO transition as it lies outside the bandpass, and found that the results agreed with those of the ACT+Planck maps to within a maximum discrepancy of 3%. The background subtraction tests in the Orion A-S area corresponded to values of 60-80 K km/s CO(J=1-0) brightness, in good agreement with the survey data shown in Figure 11 with typical values of 50-80 K km/s (Dame et al. 2001). Our procedure is robust to CO contamination in the Planck passband because of the size of the Planck beam (9.65 , Ade et al. (2014b)), which is the same size scale as the slices and the CO emission turns up overwhelmingly as a constant background level. This background is removed during the slicing and deramping process. Figure 11. Measured CO(J=1-0) emission in the Orion A molecular cloud complex from Dame et al. (2001). This analysis avoids the bright HII region in between seen on the right side of the image, and is focused on the trailing tail, where the emission levels are between 50-80 K km/s. The remaining problematic line is the HCN line, which has potential for contamination in all instruments used in these analyses. To estimate the level of contamination from this line, we scaled the conversion factor found in Ade et al. (2014c) for CO by the bandpass ratio and the emission frequencies for HCN vs CO and compared the expected signal from HCN in these regions to the excess seen in the 90 GHz maps. Estimates of the typical emission levels were found in Rydbeck et al. (1981) and Ohashi et al. (2014) which showed peaks in the densest regions of ∼ 60 K km/s and more typical levels throughout the cloud of ∼5-6. In order to account for the excess signal seen at 90 GHz, all regions would need to be in excess of 65 K km/s HCN brightness, which is more than 10x the typical levels seen across the region. In sum, we did not find any evidence that CO or HCN emission can account for the excess signal level seen in this work.

OMC 2/3
The first region we consider in detail is the OMC 2/3 region that was examined in M20. For the MBBall fit model and 120 resolution, we find that when we apply this method to the 5 slices in the OMC 2/3 region we see an excellent agreement with the model with a median χ 2 n = χ 2 /n = 1.00 for n=4 degrees of freedom. In this region the minimum was χ 2 n = 0.66 and the maximum was χ 2 n = 1.79. As did M20, we examined the 3 mm datapoints and find that, while they lie above the SED in all cases, they are much more in line with the model with an average excess of 4.9%. When examined at 25 we saw a mean brightness excess of 32%, a far more noticeable discrepancy than at 120 . This smaller difference between 90 GHz and the model contributes to the significantly better agreement and fit seen in our analysis.
We next consider the results of the MBBno90 fit, and find that the median was χ 2 n = 0.85 for n=3 degrees of freedom, the minimum was χ 2 n = 0.05, and the maximum was χ 2 n = 1.74 across the five slices in this map. The drop in the median relative to the full data fit comes as a result of slices 3 and 4, both of which are nearly perfectly fit in this regime by the model. The final model is the MBBno90 fit extrapolated down to 90 GHz which gives us a median χ 2 n = 1.65 for n=4 degrees of freedom, a minimum of χ 2 n = 0.72, and a maximum of χ 2 n = 2.77 across 5 slices with a notable increase in all of the slices. Looking at the deviation of the 90 GHz datapoint we see a mean excess of 10.9%. It appears that there is slight evidence that the 90 GHz data is elevated with respect to the model when it is constrained by that data, and that the evidence becomes significant when the model is not constrained by that data. Unfortunately, in this region, none of excesses seen broke the 4-σ threshold we have set, though 2/5 were above a 2-σ level. Even though less excess is seen at 120 as compared to the large amount seen at 25 , we do not speculate as to the cause of this as the conclusions come from disparate datasets, the slices do not coincide with each other, and we lack additional data for comparison in other regions.

Serpens
The Serpens region from the HGBS contains a total of 20 slices, of which 19 are uncontaminated (Section 4.3), split into three sub regions for data processing. Using the MBBall fit method and combining the data into one set, we obtain a median χ 2 n = 6.75 with a minimum of χ 2 n = 0.45 and a maximum of χ 2 n = 17.40. We attribute this degraded consistency with the MBB to the fact that in 17/19 slices, we find that the 3 mm data point lies above the SED curve, with an average excess of 17.1%. When we apply the MBBno90 fit to this region, we find a median χ 2 n = 1.23 for n=3 degrees of freedom, a minimum of χ 2 n = 0.41, and a maximum of χ 2 n = 5.76. In this case, we see a reduction in χ 2 n in 17/19 slices, which tracks with the excess emission seen at 90 GHz in 17/19 slices. We then applied the extrapolation to 90 GHz and found that a median χ 2 n = 11.17 for n=4 degrees of freedom, a minimum of χ 2 n = 0.67, and a maximum of χ 2 n = 24.42 with an increase compared to both other models visible in all slices. For this region, we also see an average excess emission to 34.7% relative to the model. Additionally, we found that for 13/19 slices, the excess at 90 GHz is significant to at least a 4-σ level.

Orion B
In the Orion B region we measured a total of 11 slices. The MBBall fit model applied to the data returned a median χ 2 n = 6.48 with a minimum of χ 2 n = 0.24 and a maximum of χ 2 n = 41.20. As with the Serpens region, we find that this reduced consistency with the model as compared to the OMC 2/3 region stems from the fact that most (9/11) of the 3 mm measurements are far above the SED, displaying a mean excess of 30.7%. The results of the MBBno90 fit are in line with the other regions', showing a drop to a median χ 2 n = 1.26 for n=3 degrees of freedom with a minimum of χ 2 n = 0.17 and a maximum of χ 2 n = 7.65. We also see that again, most slices (7/11) have a reduction in the χ 2 n value while the others remain nearly constant, which we expect given the excess emission seen in 9/11 slices. When extrapolated down to 90 GHz, the fit degrades to a median χ 2 n = 9.41 for n=4 degrees of freedom, a minimum of χ 2 n = 0.26, and a maximum of χ 2 n = 48.90. Additionally, the mean excess seen in the 90 GHz data points within Orion-B is increased to 56.2% and for 5/11 slices the excess is significant to at least a 4-σ level.

Orion A-S
The final region is the Orion A-S region which contains a total of 17 slices, of which 16 are not contaminated by an HII region (Section 4.3). We applied the MBBall fit to this region and found a median χ 2 n = 6.79 with a minimum of χ 2 n = 1.36 and a maximum of χ 2 n = 62.43. Again, we find that the reduced consistency with an MBB model is expected, as 15/16 of the slices show an elevated emission level at 3 mm, with an average excess of 27.1%. Applying the MBBno90 fit next, we see a median χ 2 n = 3.97 for n=3 degrees of freedom, a minimum of χ 2 n = 0.11, and a maximum of χ 2 n = 15.83. As with the rest of the regions, we see a majority of slices (13/16) with a reduction in the χ 2 n value as the model is able to better fit the short wavelength data. With these fits examined, we looked at the extrapolation of the MBBno90 fit to 90 GHz and found that this new fit had a median χ 2 n = 10.79 for n=4 degrees of freedom, a minimum of χ 2 n = 2.27, and a maximum of χ 2 n = 68.78 with an increase versus both models in 14/16 slices. Additionally, the average elevation of the 90 GHz data point had increased to 52.0% with a 4-σ or higher significance in 12/16 slices.

Overall results
In this section, we describe the trends seen across all regions for the MBBall fit, MBBno90 fit and extrapolated MBBno90 fit. For the overall MBBall fit behavior we saw a median χ 2 n = 6.19, enhanced emission in 45/51 slices, and a mean elevation of the 90 GHz data of 21.9%. When we restricted the SED to λ ≤ 2 mm the average quality of the fit was dramatically improved,  with a median χ 2 n = 1.57. The MBBno90 fit was then extrapolated down to the 90 GHz data, which resulted in a reduced quality of the fits versus both previous analyses with a median χ 2 n = 9.96 and enhanced emission in 45/51 slices. Figure 12 shows the per-slice residuals for the no90 model and, in particular, the residuals at 90 GHz are clear outliers compared to all other frequencies.
To explore this spread further and understand the significance of the result, we present the results of the Monte-Carlo noise and calibration error simulations for the 90 GHz values in particular in Figures 13 and 14. The two distributions from Figure 13 are subtracted Figure 14. Distribution of fractional excess emission relative to the model as predicted by the simulations in Figure  13. In this particular slice the mean excess is ∼50% and is quite inconsistent with zero. Figure 15. Visualization of the 90 GHz residuals by region with the associated 1-σ uncertainty in the %-excess emission relative to the model predictions. Two slices in the Orion B and one in the Orion S region are not displayed on this plot in order to preserve the readability. These slices display large excess emission values ranging from 147-307%.
from each other and scaled by the mean value of the model prediction to create the distribution shown in Figure 14. The mean and standard deviation of this new distribution are then extracted and used to calculate the tension with the null hypothesis of 0 excess emission at 90 GHz, which are shown as the "Excess σ" column for each slice in Tables 4 and 5. The results sorted by region and value, along with their associated 1-σ uncertainties are displayed in Figure 15. We find that, at the 4-σ or greater significance level, nearly 60% (30/51) of the slices are inconsistent with the modified blackbody model, which is in otherwise good agreement with the rest of our spectral bands. These probability distributions based on MC noise simulations provide clear evidence that the excesses seen in the measured data points are a real and significant effect at 90 GHz.
In addition to the above, we also examined the correlation between the β and T d values versus the excess seen at 90 GHz. We find a correlation r of 0.133 for β and the excess and a value of r=-0.010 for T d and the excess, both of which imply a lack of correlation between the variable. Thus, this analysis shows that there exists a significant and systematic elevation of the emission at 90 GHz on 120 scales with respect to the standard modified blackbody models that is inconsistent with both random fluctuations and calibration errors and is uncorrelated with the fit parameters.

Two-component model
Our first additional model to consider is the twocomponent model (see Section 3.3) that was described in M15. An example of this fit is shown in Figure 16. At 90GHz, the fits produced a lower residual than the generic modified blackbody with a mean excess of 17.0% across all regions. The high-frequency end of the spectrum, however, was not as well fit, which led to an overall median χ 2 n = 5.51 for n=2 degrees of freedom, a minimum of χ 2 n = 0.20, and a maximum of χ 2 n = 42.54. As it stands, this model seems to provide a somewhat reasonable description of the SEDs, which is to be expected for a model with more independent parameters. Unfortunately, while it does better than the single modified blackbody at 90 GHz, it performs worse over the entire spectrum on average while still under-predicting the emission at 90 GHz. In addition, the model predicts a second component with unusually high βs (mean ∼2.7) and low temperatures (9-11 K) inconsistent with both the expectation of a warm second component and the NH 3 temperatures seen in OMC 2/3 (Friesen et al. 2017). As it stands, this model, while reasonably motivated for the cold, high Galactic latitude cirrus, does not provide an adequate description of the emission within these warm and dense star-forming regions where less spectral coverage is available at high-resolution to balance the number of required parameters.

Amorphous dust model
We finally consider the results of fitting the amorphous grain dust model. This model relies on a standard blackbody with an emissivity term that is a function of temperature and frequency. For our fits, we allowed the dust temperature and the normalization to be free parameters while fixing the emissivity to scale as the 17.5 K measurements from Paradis et al. (2011), the closest match to the mean 15.2 K temperatures we found to fit the SEDs. An example fit is shown in Figure 17. This model provided a significantly worse fit to the data than the two-component model with a median Figure 16. An example of the two-component dust model fit. Overall the fit provides a reasonable estimate of many flux density measurements, but the constraint set by the 90 GHz data causes it to over-predict the 150 GHz and still under-predict the 90 GHz values. The slice shown here is arbitrary, but represents the sample well. Figure 17. Example fit using the amorphous dust model on a slice SED. In this slice, the model over-predicts the 90 GHz emission while under-predicting the 250-500µm emission, leaving the errors spread out over the SED rather than centralized at 90 GHz as with the MBB. χ 2 n = 8.13 for n=5 degrees of freedom, a minimum of χ 2 n = 0.37, and a maximum of χ 2 n = 115.1. Interestingly, this is the sole model which displays large errors in predictions across the spectrum, rather than failing at any one specific frequency repeatedly. This model fares better than any other model at predicting the surface brightness with an average excess of -1.04%, meaning it slightly over predicted the 90 GHz flux on average. The overall quality of this fit, however, is bad. We find that is it comparable to that of the 90 GHz constrained MBBall fit and is significantly worse than that of the short wavelength MBBno90 fit. There exists hope for this model, however, as more detailed (T, ν) data would allow for us to include the -temperature in the model in the future.

DISCUSSION AND CONCLUSIONS
Here we have presented a new study with significant evidence for enhanced emission, relative to a modified blackbody model, at 90 GHz in the nearby molecular cloud regions of Orion A, Orion B, and Serpens at 120 resolution. Considering first the analysis of only the OMC 2/3 region at 120 resolution and the search for evidence of enhanced emission using the ACT and Herschel data we found evidence for enhanced emission at 90 GHz, with all slices showing enhanced emission above the MBB. The data as compared to the MBBall fit model has a modest 4.9% average elevation and when compared to the extrapolated model it rose to 10.9%, though no slices within this region were inconsistent with 0 excess to our 4-σ cutoff. In this region we also saw that no90 fit produced the best average agreement with the model and that fit, when extrapolated to 90 GHz, significantly degraded the quality of the fit relative to both other methods.
With the prototype analysis of the OMC 2/3 region completed and showing tentative evidence for emission above an MBB spectrum, we extended our study to other molecular clouds within the HGBS. We chose to focus on three targets at similar distances (∼400pc) to the OMC 2/3 filament, the Orion A-S, Orion B, and Serpens regions. Within these regions we found an even higher fraction (41 of 46) of the slices showing elevated emission and a striking departure in the behavior of the 90 GHz emission, with 30/46 slices displaying a 4-σ or greater inconsistency with the modified blackbody at 90G Hz. These regions have provided much stronger evidence for enhanced emission than even the OMC 2/3 region which prompted this study and was the backbone of the results in M20.
In addition to the results of the modified blackbody SED models, we have considered the possible effects of free-free and line emission contamination as well as twocomponent and amorphous dust models to explain the 90 GHz discrepancies. In the case of free-free contamination, we found only two slices (one of which was selected due to a low-frequency source for code testing) out of the 53 total which were identified as co-located with an HII region from the NVSS+WISE catalogs. We elected to remove these slices from consideration due to the likely presence of contamination. For the case of contamination by molecular line emission, we found that the Planck data were indeed contaminated by the CO line. The Planck beam was large compared to our slices, however, so the CO line emission appeared as an approximately constant background term which was removed via the deramping. The results were in line with the ACT only maps. The other potential cause of contami-nation is the HCN line, which we can ignore given how weak the emission was. Turning now to the additional models, we found that, while they resulted in similar reduced χ 2 n values as the MBBall fit model, they each produced results which were concerning, either unreasonable parameters, poor fits, or both. In the case of the two-component model, we restricted ourselves to a fixed cold component model from M15 in combination with a free, warm component. This resulted in overall slightly improved (but still poor) fits from the MBBall fit (χ 2 n = 5.51 vs χ 2 n = 6.19), a systematic under-prediction of the 90 GHz flux (17.0% low), and temperatures inconsistent with the NH 3 in OMC 2/3 and βs averaging over 2.6, which are unusually high. On the other hand, the amorphous dust grain model provided a drastically worse fit to the curves (χ 2 n = 8.81 vs χ 2 n = 6.19) than the MBBall fit and good agreement with the model at 90 GHz (1.04% over prediction). However, this model fails to describe the high-frequency end of the spectrum well, as is particularly evident in Figure 17. Both of these models have the potential to provide a better description of the SED than a generic modified blackbody; however, in the case of the amorphous dust, more measurements of the behavior at different temperatures would be required, and for the two-component model, we would need more spectral coverage to comfortably use a six-parameter model to fit the data.
When compared to the results of M20, we found that our data reduction pipeline produced results in excellent agreement with M20, having a median difference and standard deviation of ∆T d /σ T d = -0.190±0.095 and ∆β/σ β = 0.13±0.11. The method that we used in this paper has the benefits of being computationally lightweight and having the ability to extract information on these clouds solely through existing survey and archival data. This study has shown that our pipeline is a viable method through which we can combine the data at 120 resolution and measure the SEDs across a variety of sources of varying brightness. By extracting information about the spectral energy distribution in these regions, and especially about the prevalence of enhancement at 90 GHz, we have paved the way for complementary observations and analyses to examine the emission on smaller angular scales and lower frequencies.
With our concerns about contamination assuaged and additional models fit and ruled out, we consider possible reasons for this deviation from the modified blackbody. While we are not sure what the source of this emission is, whether it is dust-based or from some other source, analyses have so far ruled out contributions from spectral line emission of CO and HCN, AME, and spinning dust variants. This leaves a gap in our understanding of the dust, though we should not be surprised that such a simple MBB model fails to encapsulate the complex behavior of the emission from these clouds. It is known already that low-frequency emission does not extrapolate well to higher frequencies (Adam et al. 2016), and perhaps the spectral range we consider is simply too wide for a single temperature or fixed spectral index modified blackbody to accurately describe. This seemingly common elevation of emission at 90 GHz is troubling, as it means that using standard dust models for core envelopes (e.g., Ossenkopf & Henning (1994)) to convert thermal dust emission at 3mm to dust masses may overestimate the measured masses significantly. Whatever the source of this emission may be, it represents a significant deviation from the generic model and therefore warrants further investigation.
In order to better understand the properties of this enhanced emission, there are three obvious paths for future investigation. The first path is to consider these regions in polarization, which would trace the dust emission in particular and could not be contaminated by free-free or other unpolarized emission. This would give us an understanding of the association with dust or a alternative contributing source. This would require sensitive surveys of the molecular clouds in polarization, on wellmatched angular scales to the currently existing ACT polarization data. The second is to further extend this study to additional clouds and regions. We have so far concerned ourselves with only a small fraction of the available HGBS and ACT data. There exist data on tens of clouds available to study with a wealth of information about this behavior to extract. Our analysis has shown a systematic increase in the surface brightness at 90 GHz, with more than 90% of the > 50 slices showing an increase even when the model is constrained by that data. We have studied only 3 relatively small regions within these cloud, however, leaving potentially hundreds of slices for future analyses. These future works would provide significant spatial coverage to test the ubiquity of the 90 GHz excess, and a much larger sample size. The third path is to delve further into the behavior of the clouds studied in this paper with the methods of M20 and data available at higher resolution. We have provided evidence that there exists enhanced emission in these regions for 120 resolution at 90 GHz, which makes them prime targets for followup observations with high resolution on instruments such as MUS-TANG2 and the GBT Ka-band receiver, the latter of which would provide the very low frequency coverage that showed the largest enhancements in M20. Additionally, where legacy spectral coverage of these clouds lapses, there exist instruments such as NIKA2 (Calvo et al. 2016), at 1.25 and 2 mm, on the IRAM 30m telescope which provide high angular resolution coverage in the portion spectrum between the GBT (3 mm) instruments and the Herschel data (160-500µm), crucial for bridging the gap in the tail of the SED between 600 GHz and 90 GHz and determining the presence of enhanced emission. With this suite of available instruments, there exists significant opportunity for followup observation and analysis of these regions to determine the nature of dust emission in star-forming molecular clouds.
Acknowledgements: IL would like to acknowledge the assistance of Marius Lungu from the ACT collaboration in finding and sharing the ACT bandpass data used within this paper. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (Grant agreement No. 851435). The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. the authors would also like to acknowledge the thoughtful and constructive suggestions of the referee which helped to improve the manuscript.