© ESO 2022

1 Introduction

Astronomers typically classify stars into two or three mass regimes for a number of reasons that can generally by summarised as their expected courses of evolution and environmental influences. High-mass stars are widely agreed to include any star of M ≳ 8 M_⊙ – the minimum initial mass of a supernova progenitor (Heger et al. 2003; Smartt 2009, and references therein) – and stars of lower mass sometimes are divided into intermediate-mass (1.5 M_⊙ < M < 8 M_⊙), and low-mass (M ≲ 1.5 M_⊙) depending on whether or not the stars are small and cool enough to sustain a convection in their outer layers (Hansen et al. 2004, Sect. 2.2). Low-mass stars make up the overwhelming majority of stars in the Galaxy, and along with intermediate-mass stars, they are the hosts of planetary systems. High-mass stars are rare and may not have planets, but they produce most of a galaxy’s light, and are major drivers of the physical and chemical evolution of the interstellar medium (ISM) via fierce winds, ionising radiation, and supernovae. Most stars of all masses form in clusters (Lada & Lada 2003), but high-mass stars do so most exclusively. This, on top of their rarity and short lifetimes, makes high-mass stars difficult to study in the process of formation; high-mass star-forming regions (SFRs) are characteristically distant, heavily extincted, and very crowded.

The formation process for high-mass stars has long been a subject of heated debate, with the main source of contention being whether high-mass stars begin as high-mass cores (~0.1 pc, that is they are core-fed) or if they are made high-mass by starting in a particularly dense environment (that is they are clump-fed, where clumps are parsec-scale condensations). The former paradigm is encapsulated in the core accretion or turbulent core model (McKee & Tan 2003), essentially a scaled up, accelerated analogue to how low- and intermediate-mass stars are understood to form. The latter paradigm started with the so-called competitive accretion model proposed by Bonnell et al. (2001, 2004), where similar-sized condensations can gain more mass if they happen to start near the centre of the cloud’s potential. This has since been largely superseded by the hierarchical collapse model (Vázquez-Semadeni et al. 2009, 2019; Padoan et al. 2020), where the large-scale cloud collapses into filaments and nodes (or ridges and hubs, respectively) and concentrates gas into parsec-scale clumps. Therein, low- to intermediate-mass cores that are born from turbulent fragmentation (Padoan & Nordlund 2002) can rapidly accrete a massive envelope. High-mass prestellar core candidates exist (see for example Russeil et al. 2010; Tackenberg et al. 2012; Csengeri et al. 2017; Tigé et al. 2017; Motte et al. 2018), but none have been confirmed, and similar candidates from earlier studies have subsequently been found to contain deeply embedded protostars. On the other hand, high-mass prestellar clumps and high-mass clumps containing protostellar cores of a wide range of masses are abundant in the literature (Bontemps et al. 2010; Svoboda et al. 2016; Sanhueza et al. 2019; Kong et al. 2021; Pitts & Barnes 2021, to name a few). This evidence would seem to favour the clump-fed scenario, but only one high-mass prestellar core would need to be confirmed to prove the core-fed scenario is also viable.

Another potential way to distinguish the above scenarios is to look at their predictions for the distribution of matter in protostellar envelopes. The radial mass density distributions for protostellar envelopes of all masses are typically modelled as a power law of the form ρ(r) ∝ r^p, where the index p is between 0 and − 2. Interpretations may vary depending on whether one subscribes to the historically favoured inside-out collapse model of Shu (1977), or the outside-in mode that hierarchical collapse appears to produce (Gómez et al. 2021). Under the former scenario, a slope of −2 is indicative of an isothermal or Bonnor-Ebert sphere, −1.5 indicates an isolated sphere in free-fall, and −1 is the slope for the density distribution of gas in virial equilibrium to where the inside-out collapse has not yet propagated (McLaughlin & Pudritz 1997). In the hierarchical collapse scenario, solving the continuity equation for a spherical envelope reveals p = − 2 to be an attractor towards which any other value of p will evolve (Vázquez-Semadeni et al. 2019). Low- to intermediate-mass protostars, tending as they do to form in less dense environments, should have envelopes well-described by a free-falling spherical envelope model. If we start with that assumption and if the core accretion model of massive SF is true, massive protostars should have similar values of p. If a hierarchical collapse model suits massive stars better, one of three possibilities may be observed:

• Evolution towards p = −2 will be quicker, so most observed massive protostars will have approximately that value of p.

• Turbulent fragmentation will cause the spherical envelope approximation to break down, and p will be flatter than observed for lower-mass stars as different cores pull the mass in different directions.

• The ridge-hub morphology of the larger cloud in the hierarchical collapse scenario will lead to accretion flows that strongly violate the assumption of spherical symmetry, which may lead to either flatter values of p or an increased dispersion in p.

Thus in most cases, if the transition from intermediate to high-mass stars is accompanied by a transition to formation by hierarchical collapse, we should observe a population-wide change in the distribution of p and any parameters derived from the density profile (for example mass).

To test the above prediction, we needed constraints provided by detailed models of the density (and temperature) distributions over a large sample of protostars of all masses, but especially in the poorly-sampled region between about 10 and 100 M_⊙ (~ 50–1000 L_⊙). Based on the relative scales of the stellar initial mass function (IMF) and the similar-shaped core mass function (CMF), it is expected that between 15 and 50%, typically around 30%, of the mass of a core will end up in the star (Alves et al. 2007; Jørgensen et al. 2007b; Benedettini et al. 2018; Könyves et al. 2020), which means that the sampling gap spans the range of core masses with the potential to form either intermediate- or high-mass stars. Maps of the high-mass SFR Cygnus-X (d ~ 1.4 kpc, Rygl et al. 2012) by Motte et al. (2007), hereafter abbreviated as MBS07 revealed a wealth of protostellar cores with masses overlapping the high end of the desired mass range. Further, subsequent studies have since resolved some of these cores into fragments with masses suggestive of intermediate-mass star progenitors (see for example Bontemps et al. 2010; Duarte-Cabral et al. 2013), although continuum data from Herschel do not resolve these fragments. We selected ten of these sources for the following analysis in this article, and in anticipation of future chemical modelling.

In Sect. 2, we discuss how we obtained and pre-processed our continuum data to get fluxes for spectral energy distribution (SED) fitting; how we used SED fitting to estimate parameters, like envelope mass and luminosity, for bench-marking 1D envelopemodels; and how we then used those SED fits and derived parameters to model the temperature and density structures of our sources. In Sect. 3 we compile and compare the results from SED-fitting and 1D envelope modelling, and compare both sets of results to values that other studies derived for the same objects. In Sect. 4, we place our results in the context of individual source morphologies, and of comparable data on as many resolved protostellar sources as we could find spanning seven decades in protostellar luminosity, to identify trends and determine if they differ between low- and high-mass protostars. We conclude with a summary of our results in Sect. 5.

2 Methods

As discussed in Sect. 1, our science goal was to establish whether or not different physics govern the evolution of high- and low-mass protostellar envelopes by comparing trends in their structural parameters at roughly the same evolutionary phase, but over a wide range of masses and luminosities. Our technical goals include filling the sampling gap in the literature at intermediate masses and luminosities, which we were well-positioned to do thanks to concurrent work on 10 moderately-massive protostellar sources in Cygnus-X. However, extracting envelope parameters from available data is a multi-stage process; the temperature and density profiles of protostellar envelopes cannot be directly observed over much of a protostellar envelope. Resolution limits aside, the dominant gas and lightest molecule in any protostellar envelope, H₂, does not have a permanent dipole moment and therefore cannot emit at temperatures below about 500 K. Other tracers are needed to determine the density distribution and detect the protostellar emission as it emerges, reprocessed, from the cold envelope. Our tracer of choice is thermal dust emission, which has fairly large random uncertainties in its abundance relative to H₂ (as explained inSect. 2.3) but is robust to a much wider range of temperatures than most molecular tracers.

The procedure described in the following subsections has three main phases: photometry, SED-fitting, and spherical envelope modelling. The necessity of photometry becomes immediately evident in the phases after it, but the results of SED-fitting and spherical envelope modelling may at first glance look redundant. There are several reasons why we do both. Spherical envelope modelling incorporates more detailed physics through ray-tracing, and is able to recover radial structure that SED-fitting collapses to a point. However, the parameter space for spherical envelope modelling ismuch too large to optimise over without prior constraints. The parameters from SED fitting and their uncertainties provide those constraints. Dust SED fitting is also the first line of study for many pre- and protostellar objects, andis a way to derive mass-weighted average parameters that can be compared in lieu of radially resolved models of the envelope structure, if such models are not feasible with the available data. As a long-term goal, we want to make both sets of data available so that similar future comparisons of low and high-mass protostellar sources can include our data, regardless of whether or not the structure of individual protostellar envelopes can be resolved.

Fig. 1 SMA 868 μm (345 GHz) images of all 10 PILS-Cygnus sources (N53 and N54 are both featured in the middle right panel) contoured in red at 10, 40, and 70% of the maximum flux. The cyan ellipses in the upper right corners show the FWHMa of both axes and positionangle of the synthesised SMA beam, and the line segments alone the bottoms of each panel show how long 10 000 AU and its corresponding angular length appears at the adopted distance to each object. We note the change of scale for S26. All images are normalised to their peak flux densities, which are given in Table 1.

2.1 Observational data

The objects we contribute to this metastudy come from the PILS-Cygnus survey (PI: Lars E. Kristensen), for which 10 protostellar envelopes in Cygnus-X were selected based on brightness in continuum, SiO 2–1, and H₂ O 2₀₂ − 1₁₁ emission (Motte et al. 2007; Mottram et al. 2014). These were mapped with the SubMillimeter Array (SMA) telescope on Mauna Kea in the continuum and a variety of molecular species. The SMA continuum observations, taken in June and November of 2017, covered a 32 GHz-wide band centred at 345 GHz (868 μm) in both compact and extended array configurations for a combined spatial resolution of 0′′.7 and sensitivity of 0.15 Jy km s⁻¹ (see van der Walt et al. 2021, for detailed description of data reduction). Figure 1 shows the 868 μm (345 GHz) continuum images of all 10 sources, contoured at 10, 40, and 70% of the maximum flux in each image, with the synthesised SMA beam dimensions shown in the upper right of each panel and the image scales along the bottom of each panel. These images were not used directly in dust SED fitting because much of the outer envelope emission was resolved out. They were, however, subject to Gaussian decomposition to separately determine the number of core fragments inside each envelope as a check against Bontemps et al. (2010) and to compare core multiplicity to structural envelope parameters in search of trends.

Table 1 gives the sky coordinates and designations per MBS07 for all 10 objects, and Fig. 2 shows Herschel RGB cutouts around all of the sources at 70, 160, and 250 μm, with the sources and nearby objects from the DR catalogue (Downes & Rinehart 1966) labelled. Of those 10 sources, nine could be identified in archival mid-infrared (MIR) to submillimetre (submm) broadband continuum data, and eight could be separated from the extended filamentary emission and/or other sources in the DR 21 Ridge. N38 could not be conclusively identified at any wavelength, and N54 was too faint beyond 160 μm to construct a believable SED.

In most calculations, we adopted the average distance to Cygnus-X of 1.4 kpc from Rygl et al. (2012), except for S26, for which Rygl et al. (2012) finds a distance of 3.3 kpc. We used more precise values from Rygl et al. (2012) for individual objects, where available, to estimate the projected distances from Cygnus OB2 (hereafter abbreviated Cyg OB2). However, distance estimates to the centre of Cyg OB2 range from about 1.4 kpc (Hanson 2003) to about 1.7 kpc (Massey & Thompson 1991; Berlanas et al. 2019), and the recent work by Berlanas et al. (2019) using Gaia parallaxes suggest both estimates may be correct because Cyg-OB2 may be two OB associations, one each at 1.4 and 1.7 kpc, projected on top of each other. This revelation deeply complicates any interpretation of the PILS-Cygnus source morphologies based on deprojected distances to Cyg OB2. For now, we only report the projected distances in Table 1, which should be interpreted as lower limits on the absolute separations from Cyg OB2.

We concentrated mainly on the FIR/submm continuum data from Herschel/PACS (Poglitsch et al. 2010) and Herschel/SPIRE (Griffin et al. 2010) as part of the Hi-GAL survey (Molinari et al. 2016), and JCMT/SCUBA(-2) (Di Francesco et al. 2008; Holland et al. 2013) for SED computation. To compute bolometeric luminosities, we supplemented these data with photometry from 2MASS (Skrutskie et al. 2006), Spitzer/IRAC (Fazio et al. 2004; Beerer et al. 2010), Spitzer/MIPS (Rieke et al. 2004), MSX (Egan et al. 2003), and upper limits from IRAS (Neugebauer et al. 1984). We performed most of the photometry ourselves; for MSX and IRAS, whose large uncertainties limited their influence upon the fit, we used point-source catalogue fluxes. Table 2 shows the instrument and filter specifications for every data set we used for any purpose, photometry or otherwise, including the corresponding wavelengths, resolutions, and pixel scales. All fluxes and their associated uncertainties are listed in Appendix A, Table A.1.

For IRAC and 2MASS data, in which the sources were largely resolved, we performed aperture photometry using the Astropy-adjacent Photutils package with some adaptations for the highly crowded surroundings. Instead of background annuli, which could not be configured to avoid all neighbouring sources, we took the average of two background apertures of identical size to the target aperture, placed in uncontaminated areas on opposite sides of the target. Uncertainties in the fluxes listed in Appendix A are the statistical uncertainties in the differences between the target aperture flux and the background aperture flux with the uncertainties in each aperture flux propagated through. Systematic uncertainties due to the presence of filamentary emission or absorption are likely to be larger. S8 in particular has a complex extended morphology, like a comma embedded in a Y-shape, that made it difficult to decide what should be included or excluded from the source. The automatic photometry results were remarkably discontinuous with the rest of the SED, sowe redid the photometry for it manually in SAO DS9, following the same procedures and using the 90% contour levels as a guide for the aperture size.

2.2 Deconvolution of sources

Because many of the sources are unresolved and separated from each other or nearby objects by less than the beam FWHM of most Herschel bands, we developed a simple Gaussian deconvolution and fitting routine with all parameters free to vary within strict limits determined by eye (for example the coordinates of the centre of each Gaussian were allowed to fall anywhere within an area the size of each image’s resolution element, centred on the coordinates given in Motte et al. 2007). It fits up to four components and a flat background to cutouts with dimensions tailored to enclose no more than four targets, but also extend at least 10σ away from any source in at least two directions. The fluxes used to compute the SEDs and all derived parameters use the best-fit Gaussian integral of the deconvolved, background-less target flux.

Since the objects often did not separate cleanly at the longest wavelengths, and filamentary contamination remained an issue, RMS deviations from the observations were typically used in lieu of the unrealistically small calculated fitting uncertainties. We tried to estimate the uncertainties by distributing the RMS deviations amongst the modelled sources according to their fluxes, but this was not possible with N30, N38, and N54. With N30, it was unclear whether a one- or two-component Gaussian was preferable at short wavelengths, and at long wavelengths, most of the flux that should have been assigned to N30 was being assigned to a minor object identified as N31 in MBS07. MIR images from Spitzer suggest N31 is either a filamentary spur or some other distortion of N30, so we opted for the one-component model in the hopes that the inclusion of the spur and the relatively high fitted background would roughly balance out. In FIR/submm images where N38 would have been resolved if it were visible, the nearest source was 26′′ away, a much greater distance than the position uncertainties for any of the aforementioned observing instruments and only a couple arcseconds shy of the distance between N38 and N48. N54’s location suggested that it would be blended with N53 at all FIR/submm wavelengths, but at the wavelengths in which their centres should have been more than a FWHM apart (ergo, able to be deconvolved), N54 is entirely swamped by both N53 and another source on the opposite side that was not part of our sample, N52. The total flux from N54 was within the margins of error for N53 and N52 at PACS wavelengths, and at longer wavelengths the fit was sometimes negative before we explicitly restricted the Gaussian models to positive amplitudes. Like N38, the fits from N54 could not be trusted. At 70 and 160 μm, blending with N54 is not expected to contribute more flux to N53 than is accounted for by its large error margins. The same seems likely at longer wavelengths because N54 is quite blue in RGB images using 250, 160, and 70 μm data, but our data do not allow confirmation.

Fig. 2 Herschel three-colour cut-outs showing the PILS-Cygnus sources in their natal environments at 70, 160, and 250 μm, with the DR objects labelled. Herschel beams are enlarged by two-fold for visibility.

2.3 SED-fitting

The SED fitting routine used for the deconvolved data is described in detail in Pitts et al. (2019). The program uses MPFIT (Markwardt 2009) to derive an initial dust temperature (T_dust) and molecular hydrogen column density ($N_{{\text{H}}_2}$), interpolates on grids of colour correction factors¹, propagates the estimated uncertainty in the colour correction factor to the input uncertainty, and recomputes the SED from the colour-corrected data. The program normally assumes a single modified blackbody, or greybody, component along the line of sight following the prescription of Hildebrand (1983). Up to three additional components of the same functional form may be added, provided they have sufficiently distinct peak temperatures. The greybody equation takes the form $$I_{\nu }\approx B_{\nu }(T_{\text{dust}})\left(1-\text{exp}\left(-{\left(\frac{\nu }{{\nu }_0}\right)}^{\beta }\frac{{\kappa }_0{N}_{{\text{H}}_2}}{\gamma }\mu m_{\text{H}}\right)\right)$$(1)

where B_ν (T_dust) is the Planck function, β is the dust emissivity index (typically with 1 ≲ β ≲ 2), κ₀ is the dust opacity at reference frequency ν₀, γ is the gas-to-dust ratio (assumed to be 100), μ is the mean molecular weight per hydrogen molecule (2.8; see Appendix A of Kauffmann et al. 2008), and m_H is the mass of a hydrogen atom. For each object, we took the best-fitting value of κ₀ at 350 μm (ν₀ = 856 GHz) for grains with thin ice mantles from Ossenkopf & Henning (1994). Note that Eq. (2.3) does not assume optically thin dust emission, because for most PILS-Cygnus sources, optical depth at 70 μm is not much less than 1.

Ordinarily T_dust is somewhat degenerate with β and $N_{{\text{H}}_2}$ is degenerate with both γ and κ₀, so only T_dust and $N_{{\text{H}}_2}$ were initially left as free parameters. However, we found that we could get β to converge by disabling colour correction and fitting only data at wavelengths longer than either 70 or 160 μm depending on whether or not the source was warmer than about 30 K. For N30, N51, and S26, which are also prominent at shorter wavelengths, we first fit only data at λ ≳ 70 μm to derive β, then fixed β and fit a sum of greybodies with the understanding that the warmer components would only be first-order approximations of the luminosity. Since the uncertainty in β dominates over all other factors, we estimated the uncertainties in the temperature and column density parameters for the warm and/or hot components by varying β within its uncertainty margins from the single-component fit, and taking the range of output values as the margins of uncertainty for the other parameters. We emphasise that the warmer temperature components are probably not physical, and that the tendency of such components to separate around 10 μm is an artefact of the silicate absorption feature. Except for the purpose of calculating bolometric luminosities, we are only interested in the parameters of the cold thermal dust SED, and warmer gas accounts for ≪ 1% of the total mass.

We also note that some sources – N12, N51, N53, and N63 – have had emission at 70 μm excluded or assigned to a different temperature component of the SED fit. In these cases, the 70 μm emission may be inconsistent with the rest of the cold(est) dust SED component due to contamination by nonthermal emission from very small grains (see, for example Desert et al. 1990; Bernard et al. 2008 for more details). That the rest of the sources have 70 μm emission consistent with the SED fit to data at longer wavelengths does not indicate a lack of such contamination. It only indicates either that the amount of contamination is at the level of noise, or that the sum of thermal and nonthermal contributions at 70 μm lie within the expected range of values possible for a single cold greybody SED given the uncertainty in β.

As usual, γ and κ₀ are thought to be uncertain by at least a factor of about two and 2.5, respectively (Ossenkopf & Henning 1994; Beckwith et al. 1990; Motte et al. 2007; Reach et al. 2015, in particular suggest the gas-to-dust-ratio may drop by more than a factor of 3 in cold clumps/cores), so the uncertainties in T_dust and $N_{{\text{H}}_2}$ and all derived quantities refer to those contributed by the observational data, with the factor of two to three left implicit. This code was originally intended for parsec-scale prestellar clumps, not centrally-heated cores, so the temperature and density results may be best interpreted primarily as initial values for Transphere. To compute the masses and luminosities from SED-fitting, we used the 3σ mean radii from MBS07 corrected by the distances from Rygl et al. (2012).

Finally, we remark that this SED fitting routine does not model the broad 9.7 μm silicate absorption feature or the prominent emission lines attributed to polycyclic aromatic hydrocarbons (PAHs) between 3 and 20 μm. Modelling any of these features would have changed the total luminosity by much less than its margins of uncertainty, and no other parameter besides the total luminosity depends significantly on accurately characterising the emission at λ < 50 μm.

2.4 Modelling with Transphere

Transphere (Dullemond et al. 2002) can, among other things, calculate the dust radiative transfer in 1D through a spherical envelope using user-specified envelope masses, luminosities, and density power law indices, and output radial temperature profiles and SED fits at a reference radius. We excluded data at λ < 50 μm from fitting because at those wavelengths, the emission is dominated by the deep interior which is likely to be in a disk, ergo the approximation of spherical symmetry breaks down (Jørgensen et al. 2007a; Lommen et al. 2008; Kristensen et al. 2012, Appendix C). The important user-supplied parameters are listed, and their effects on the fitted SED and radial temperature profiles described, as follows:

Stellar luminosity L_TP controls the amplitude of the SED and the position of the peak in wavelength-space via the usual L_TP ∝ T⁴ and Wien’s Displacement Law. Increasing L_TP also boosts the temperature profile at all radii. Instead of L_bol, we used only the luminosity calculated from Herschel data and longer wavelengths as initial values of L_TP because of the aforementioned breakdown of spherical symmetry at shorter wavelengths.

Stellar surface temperature T⋆ has no noticeable effects on the envelope SED if it is initialised within a range that allows the temperature profile to converge (typically between 3000 and 10 000 K, endpoints not always inclusive). Indeed, we should not expect stellar surface temperature to matter because, for a protostar embedded in a spherical envelope, all surface emission should be reprocessed by dust in the envelope before it escapes. Increasing T_⋆ does slightly steepen the radial temperatureprofile at small (r ≲ 100 AU) radii, but this is not relevant to our studies because the assumption of spherical symmetry breaks down at these radii. We fixed T_⋆ = 5000 K for all PILS-Cygnus sources.

Envelope mass M_env has a rather complex effect on the SED. Increasing M_env, all else equal, pushes the SED towards longer wavelengths and makes the SED both slightly taller and narrower due to increased reprocessing of photons from the central source. The radial temperature profile also steepens with increasing M_env, all else equal, reaching higher temperatures at small radii and lower temperatures at large radii (greater than a few × 10³ AU). This change in temperature profile reflects the radiative trapping effects of the increased opacity required at small radii (Hartmann 1998). Initial envelope mass values were set to the masses calculated from dust SED fitting and varied iteratively with luminosity.

Power-law index p does affect the envelope SED, but the variation is within the data uncertainties for the entire plausible range of p. Increasing the magnitude of p (making it more negative) with all else equal will steepen the radial temperature profile for radii less than a few × 10³ AU, but temperature profiles for different values of p converge at larger radii. We initialised p to − 1.8 for all PILS-Cygnus sources, and then adjusted it iteratively.

Inner envelope radius r_in dictates how close to the protostar we truncate the envelope, and determines whether or how much of the inner envelope is sampled with more finely-spaced bins in radius to resolve optically thick gas. By default the program switches to a finer radial grid at 50 AU; temperature profiles failed to converge for any value of r_in within about 20 AU of this radius, and contained one or more discontinuities for r_in between about 80 and 200 AU. We fixed r_in at 20 AU to avoid this issue and err on the side of including material that a more sophisticated model might place in an inner disk.

Outer envelope radius r_out has a similar effect on the SED as M_env, but not identical: increasing r_out makes the SED narrower and very slightly taller, but does not significantly change the wavelength of the peak if the envelope mass and all other variables are held constant. On the radial temperature profile, varying r_out has the effect of changing the concavity in log-log space. We initialised r_out to 20 000 AU and fit it iterativelyalong with p.

Interstellar Radiation Field (ISRF) strength affects the height of the SED at λ ≲ 50 μm and the temperature profile in the outermost few × 10³ AU (out of 1to a few × 10⁴ AU). Increasing the ISRF strength can boost the SED at λ ≲ 50 μm, but never enough to explain the observed fluxes at those wavelengths without convergence errors. We fixed the ISRF strength at 1 G₀ with G₀ in Habing (1968) units.

To find the optimal envelope masses, luminosities, and density profiles of the PILS-Cygnus sources, we took an iterative χ² -minimisation approach using models computed with Transphere. First, we computed a large grid of model SEDs over a wide range of masses, luminosities, and density power law indices for a handful of possible outer envelope radii (r_out). Second, we minimised χ² between the model SEDs and data with λ > 50 μm by interpolating on a grid of χ² as a function ofluminosity (L_TP) and envelope mass (M_env) only, given afixed p. Third, we created smaller grids of model intensity images for each object, given L_TP and M_env, for all reasonable values of p and r_out. Fourth, we interpolated on a grid of χ² as a function of p and r_out to minimise χ² between the model images and close-cropped cutouts from the observational data. Finally, we took the resulting p and r_out, and repeated the second through fourth steps until L_TP, M_env, p, and r_out converged. Convergence typically occurred by the third iteration.

Data at λ ≲ 30 μm were excluded from the fit because these data mostly trace warm (T_dust ≳ 100 K) dust and gas close to the protostar where it is likely to have collapsed into a disk. The disk region around a protostar is not well-approximated by a 1D spherical model (see, for example Jørgensen et al. 2002 or Appendix C of Kristensen et al. 2012). Accordingly, because we have not attempted to model and subtract off a disk component, the envelope masses from Transphere modelling are likely to be upper limits, especially for more evolved sources with brighter MIR emission components. As discussed in Appendix C of Kristensen et al. (2012), not subtracting the warm inner disk effectively boosts the steepness of a protostar’s flux profile by convolving the (marginally) resolved outer envelope with the unresolved inner disk. This in turn is expected to inflate the power-law index of the density profile, so those too are expected to be somewhat overestimated.

The MIR to submm data we have cannot be used to determine the ISRF independently, so we did not include it in our grid model. From trial and error, we did find that N12 is slightly better fit when some protostellar luminosity is exchanged for a stronger ISRF, since the ISRF selectively boosts the part of the SED to the short-wavelength-side of the peak. Indeed, most PILS-Cygnus sources are expected to be exposed to a strong ISRF due to the proximity (≲300 pc) of various OB associations, especially Cygnus-OB2 (see cyan arrows in Fig. 2 pointing towards OB clusters). However, in these cases, the ISRF could not be boosted enough at the expense of intrinsic luminosity to fit data at λ ≲ 30 μm without convergence errors in Transphere, most likely because any gas that is far enough out and low enough in density to feel the effects of the ISRF is too rarefied to account for any significant fraction of the mass or emission. Moreover, the widths of the dust filaments and pillars in Fig. 2 suggest all the sources may be well-shielded out to radii well beyond r_out, even N12. This would imply that most or all of the MIR emission that could not be fit with Transphere is coming from a hot inner disk heated by the protostar itself, rather than any external emission from tens to hundreds of pc away.

3 Results and analysis

In this section we discuss the results of both SED-fitting and modelling with Transphere, starting with the fitting parameters in Sect. 3.1. We compare the results for each type of analysis, and discuss potential sources of discrepancies in Sect. 3.2. Then compare both sets of results to literature covering one or more of our eight targets, identifying which set of our results is more appropriate for comparison and justifying any significant inconsistencies (Sect. 3.3.1). Finally, we place our sources in the context of literature on similarly-young protostars across a wide range of masses and luminosities, and show how we determined the significance (or lack thereof) of apparent differences in trends between different mass regimes (Sect. 3.3.2).

3.1 SED-fitting and Transphere modelling parameters

Table 3 lists the single or multi-component greybody SED-fitting parameters for all sources. Half of the sources were fit with a single component, denoted as the “cold” component in sources where multiple components could be fit. The columns, from left to right, are: source name, the 3σ mean radius from MBS07 corrected for distance measurements from Rygl et al. (2012) (r_1.2mm), dust temperature of the cold (T_dust < 60 K) temperature component, H₂ column density in the cold dust component (N[H₂]), the dust emissivity index β, dust temperature and (H+)H₂ column density of any additional fitting components (T_{dust, i>1} and N_i>1[(H+)H₂], respectively) listed in order of increasing T_dust. The warmer temperature components rely on limited data, so they should be treated with caution, especially those peaking in the poorly-sampled region between 25 and 70 μm.

Table 4 gives the FIR luminosities and envelope masses derived from greybody fitting, the bolometric luminosity derived from spline-fitting the data in Table A.1, and the fitting parameters for Transphere modelling. The columns, from left to right, are: source name, bolometric luminosity from spline interpolation of data from 1.22 to 850 μm (L_bol; see also Fig. 3), luminosity of the FIR or cold dust component derived from greybody fitting (L_FIR), total gas mass derived from $N_{{\text{H}}_2}$ within the 3σ source area given the source dimensions in MBS07 at 1.2 mm (M_c), protostellar luminosity from Transphere modelling (L_TP), protostellar envelope mass from Transphere modelling (M_env), the fitted power law index of the Transphere model’s density profile (p), and the deconvolved fitted radius in the 450 μm SCUBA-2 data (r_out). Because of the wavelength restrictions on Transphere modelling, L_TP from greybody fitting is more directly comparable to L_FIR than L_bol. Note also that M_c is not directly comparable to M_env because it incorporates the warmer temperature components and makes different assumptions about the dimensions of each source. All temperature components warmer than the coldest component combined contribute < 1% of the total mass of the object, well within the uncertainties in mass, so M_env were not systematically smaller.

We used a range of possible r_out, defined from the minimum and maximum measured radii of flux contours at 10% of the background-subtracted maximum flux in the 450 μm images, before deconvolution with the 450 μm SCUBA-2 primary beam, to minimise χ² in p and r_out relative to the model radial flux profile. For more symmetric sources (flagged ‘A’ in Table 4) we radially averaged flux profiles from the 450 μm SCUBA-2 data, while for sources with bright filamentary protrusions (flagged ‘B’), we selected a radial cut free of such structures. However, r_out turned out to be at least somewhat degenerate with all of the other Transphere fitting parameters. The resulting r_out fits, rounded to the nearest 5000 AU due to precision limits and then deconvolved by subtracting the 1.64σ 450 μm SCUBA-2 beam radiusin quadrature², produced better matches to the dust SEDs than a simple average of the measured minimum and maximum half-widths at 10% of the maximum flux. However, the uncertainties on the resulting r_out values were ill-defined – for small sources like N48, refining r_out by 1000 AU produced noticeable improvements, but for most sources we could only tell that the uncertainties were at least 20% of r_out. The uncertainties reported in Table 4 are the greater of either ± 0.2r_out, or the deconvolved differences between the fitted r_out and the upper and lower measured limits used to constrain the fitting. For χ² minimisation in envelope mass and luminosity, it was enough to adopt the nearest multiple of 10⁴ AU and refine the r_out by eye using the SED (see Sect. 2.4). The uncertainties on the masses and FIR luminosities from greybody fitting are propagated from the uncertainties in the parameters from Table 3. Uncertainties for quantities derived from Transphere modelling are only approximations from the ~ 95% confidence levels on the mapsof χ².

N12 inparticular was almost as well-fit with PACS-70 μm data as without, but the physical conditions derived in each case are dramatically different. Including PACS-70 μm data increases T_dust to 31 ± 2 K, and reduces β to 0.6 ± 0.1 and $N_{{\text{H}}_2}$ to 4.3 ± 0.6 × 10²² cm⁻². This fit is slightly poorer, but still within the uncertainties of the data. For consistency with the adopted parameters from Transphere modelling, however, we adopt the lower-T_dust, higher $N_{{\text{H}}_2}$ (higher mass, lower luminosity) results.

Fig. 3 Full SEDs for all sources, which were fitted, interpolated, and integrated by trapezoidal sums to calculate the bolometric luminosities. The curves running through data at λ ≥ 70 μm are the FIR/submm component of the SED fits discussed in Sect. 2.3, while the rest of the data has been log-linearly interpolated in lieu of integrating over the less-physical warmer SED components. For N30 and N51, no SED components passed through the 70 μm data point, so the SED fit was only used for λ ≥ 160 μm while the 70 μm data were linearly interpolated like the shorter wavelength data. The light red filled curves show the margins of uncertainty as constrained by the flux errors, and were likewise integrated to estimate the total uncertainty in luminosity.

Fig. 4 SEDs of N30. The dashed red line is the modified greybody fit to the data, shown as black and grey markers with error bars. The solid blue line is the best fit Transphere SED to only the data shown in black. The grey points at λ < 70 μm were excluded from the fit because the SED at these wavelengths likely traces a warm inner disk where the assumption of spherical symmetry breaks down.

3.2 Comparison of methods

Figure 4 shows the SEDs, both from greybody fitting and Transphere, of N30, and Fig. 5 shows N30’s radial flux, temperature, and density profiles given the parameters in Table 4. Similar figures for all other objects can be found in Appendices B and C. It is unclear from either SED fit whether the 70 μm data should have been excluded: both fits tracked each other closely at that wavelength for 5 out of eight sources (the objects with discrepancies between the two SED models at 70 μm are N48, N51, and N63). It is noteworthy that for the five lowest-T_dust sources, the β index of a (fixed-T_dust) SED fit to just the 450 and 850 μm data would be steeper than the measured β based primarily on the fits to Herschel data. We have considered a variety of explanations, none of which appear to work by themselves. It is highly likely that the shallow β coming predominantly from the Herschel data is due to our increasing inability to separate emission from the cores and that from their surrounding filaments. If that were the only factor, however, the SCUBA data should be either shallower still or entirely below the rest of the dust SED.

Half of our sources have mutually consistent M_c and M_env, but all eight have mutually consistent L_FIR and L_TP. The discrepancies between the envelope masses from each method can be at least partly accounted for by noting that $M_{\text{env}}\propto r_{\text{out}}^{3-p}$: for six of the eight sources, r_out as measured at 450 μm is substantially larger than the same at 1.2 mm (denoted r_1.2mm in Table 3) despite the sources supposedly being resolved at both 450 μm and 1.2 mm. We suspect this is because the 1.2 mm data resolve out some of the more extended structure that would be included at 450 μm, and it was often less clear what to exclude at the shorter wavelength. The two sources that were measured to be smaller at 450 μm, N48 and N53, are in crowded high-background areas, and may be over-subtracted at 450 μm or blended at 1.2 mm.

Some of the discrepancy may also come from the power-law indices for each method and how they relate mathematically (if not physically). Physically, there is little reason to expect that β and p should be related, and Transphere does not fit or incorporate β explicitly, but mathematical effects of β and p on the SED are similar: the smaller their absolute values, the wider the SED. When comparing radial profiles from Transphere modelling to the radial profiles of the sources at 450 μm, neither smoothing nor Gaussian deconvolution of sources could be used without altering the derived density profile, so we derived p from the “best” (lowest-χ² relative to the model) background-subtracted radial cut in high-contamination areas, and used the radially-averaged p in nominally isolated sources. Even in relatively isolated sources, p for the average and “best” radial cuts could be different by up to 50%. For greybody fitting, because the β is an average over a model of the whole source, β often more closely resembles the radially-averaged |p|. In fact, under certain conditions (namely optically thin dust emission, for T > 10 K, in the Rayleigh-Jeans approximation), p and β are physically related because of how intensity can be defined as a function of projected distance from the centre of a centrally-heated source. Under these assumptions, if the density and temperature profiles can be described as power laws of indices p and q respectively, then the intensity as a function of projected radial distance is also a power law with an index of 1 − (p + q) (Scoville & Kwan 1976; Beltrán et al. 2002). Then in the optically thin Rayleigh-Jeans limit, it can be shown analytically that the temperature power-law index is related to the dust emissivity index via q = 2∕(4 + β) (Chandler et al. 1998; Kenyon et al. 1993).

We had convolved the SCUBA legacy data³ to the radius of the secondary or “error” beam, which has a FWHM 40′′, and then attempted the same Gaussian deconvolution routine as with the Herschel data. Because these data were much more aggressively background-subtracted than Herschel data, we allowed the fitted flat background to accommodate negative average background levels so we could add it to the peak amplitude of the Gaussian fit. We still integrated over a Gaussian with zero constant background level. Nevertheless, we know that negative average background levels are more of an issue at 450 μm than 850 μm, so data at the shorter wavelength may have still been boosted too much despite removing the flat background from integration. That makes it less clear whether β from the SCUBA data alone would be more reliable than the values of β we derived, which are dominated by Herschel data. For most objects that were isolated or separable, the results seemed reasonable to within the 50% uncertainty of the 450 μm data as compared to the 350 and 500 μm Herschel/SPIRE data. We also checked the results against simple aperture photometry with SAO DS9 for both SCUBA and SCUBA-2 data, and with Herschel images where the clumps could be isolated. For S8 alone, we redid the photometry manually, cutting out the bright spiral-arm-like feature identified in MBS07 as S7. Unfortunately, PSF fitting of the SCUBA data was never an option since, at the native SCUBA and SCUBA-2 resolutions, all of the objects are resolved along at least one axis.

Given that the Transphere models incorporate more physics, the masses from those models are probably more accurate than the masses from multi-component greybody fitting, of which some components are not necessarily physical. The luminosities could only be fully accounted for by spline-fitting the SED as shown in Fig. 3; until the degeneracy between intrinsic luminosity and ISRF strength is broken, we will have to evaluate which luminosity is more consistent with other results on a case-by-case basis. Thus, wherever we adopt values for chemical modelling purposes, we prefer the Transphere model values for everything except the luminosity, for which L_bol is preferred given the relative lack of underlying assumptions.

Fig. 5 A. Radial emission profile comparison between the data (red dashed line) and the Transphere model (blue solid line), both normalised by their respective maxima. B. Resulting model temperature (red solid line and red axis) and density (blue dash-dotted line and blue axis) profiles.

3.3 Comparison with literature

3.3.1 Other studies of Cygnus-X

Among literature where the same objects are analysed at comparable size scales (as opposed to separately evaluating components of each object identified at mm wavelengths), there is precious little consistency between different groups using the same methods, let alone between different methods. Most of the scatter stems from the lack of consensus on how to define the sizes of these objects, upon which the masses in particular sensitively depend. In some cases the discrepancies are further exacerbated by the use of different power-law indices β and p, either because they were fixed parameters (as with, for example Duarte-Cabral et al. 2013; Ellsworth-Bowers et al. 2015; Cao et al. 2019; Ortiz-León et al. 2021) or because they were fitted using a different technique (for example Palau et al. 2014). We find that in general, the results of greybody fitting agree slightly more often with literature than the results of Transphere modelling, in large part because most other groups also performed observational SED-fitting of some kind, while very few attempted to model density profiles.

After correcting for their older distance measurements, we found that 6 out of 8 sources had at least one mass measurement consistent with the respective mass estimates in MBS07. For N30 and N48, neither M_c nor M_env were as large as the envelope masses computed by MBS07, though the discrepancy with N30 was small. The rest of the sources seem to agree more closely with the MBS07 masses when they have lower T_dust, since the temperatures we found were typically higher than MBS07 assume. All of this makes sense given the methods they used to determine mass and luminosity – they had no resolved data between 21 and 350 μm, so they derived the mass ranges using only their 1.2 mm fluxes and assumed temperatures between 15 and 25 K. Their assumed temperatures were usually lower than we found, and the opacity coefficient they use folds in an assumption that β = 1.8 (Ossenkopf & Henning 1994), which we found to be too steep. For a given temperature, a slightly steeper β significantly increases any fitted column or surface density parameters.

The most direct and comprehensive comparison we can make to our greybody fitting results is with Cao et al. (2019), who included all of the PILS-Cygnus targets except for S26 in their sample. They decomposed the emission into different spatial scales, fit single-component SEDs pixel-by-pixel to only data with λ > 70 μm assuming β = 2 and optically thin emission, and excluded data that were too far from the resulting curve (which usually included the 70 μm data). We use the same angular size measurements, but our method of deconvolving and fitting the targets is much more basic. One would expect that they would derive lower T_dust and luminosities, and our masses would be similar to or larger than their estimates. Indeed, we systematically derive higher temperatures for five of the seven sources largely because of our inclusion of the 70 μm data. However, wefound that our SED fits yielded significantly lower masses for all objects except for N51 and N53, which were within the uncertainties. Instead, their luminosities matched our greybody-derived luminosities within the uncertainties for four objects, and for a furthertwo (N51 and N53), the luminosities we derived with Transphere matched instead. Only N48 had a dramatically discrepant luminosity, and none of its other derived parameters matched, either–not unexpected given N48’s large uncertainties and our aggressive decontamination efforts. We believe most of the mass discrepancies are related to our different values of β, as discussed above with respect to MBS07; increasing the temperature decreases the column density required to match the amplitude of the SED further.

Only one other group attempted density profile fitting to derive p for more than one object in our sample. Palau et al. (2014) simultaneously performed the density profile fitting and SED-fitting with a single temperature component fit to all data with λ > 50 μm, including free-free cm-wave radiation where it was available. Their single-component SED fits were usually offset towards higher temperatures relative to both their data and our coldest SED components for the same objects, and it is unclear why that would happen with their wavelength restrictions. Our Transphere models used the same lower wavelength limit, and while the Transphere models were often narrower than the greybody fits, the peak positions of the two fits were seldom misaligned. The fitted values of p from Palau et al. (2014) are also systematically greater in magnitude, as are their β values except in the case of S26. These results would make sense if the narrower SED fits of Palau et al. (2014) were because their fits systematically ignored the 70 μm Herschel data, but except in the case of S26, their model SEDs appear to fit the 70 μm data and run parallel to but offset from the longer-wavelength data.

Fig. 6 Comparison of PILS-Cygnus Sources (large black circles with error bars) to objects in the literature from van der Tak et al. (2000) (purple x’s), Jørgensen et al. (2002) (magenta upward-pointing triangles), Mueller et al. (2002) (thin blue diamonds), Shirley et al. (2002) (maroon stars), Hatchell & van der Tak (2003) (green squares), Williams et al. (2005) (small orange circles), Crimier et al. (2009, 2010a,b) (wide navy diamonds), Kristensen et al. (2012) (cyan downward-pointing triangles), van der Tak et al. (2013) (red plus signs), and Cao et al. (2019) (pale gold error bars). Solid and dashed black lines represent single and broken power law fits to the data, respectively.

3.3.2 Comparison to protostars of all sizes

Following the example of Crimier et al. (2010a), in Fig. 6 we have plotted the Transphere-derived envelope masses, χ² -fitted 450-μm outer envelope radii, density power-law index, and gas number density at r = 1000 AU against bolometric luminosity for all of the PILS-Cygnus sources (black circles with error bars in Fig. 6). For comparison, we include similar results for about 180 other resolved protostellar sources ranging from < 0.1 to > 10⁴ M_⊙ from van der Tak et al. (2000, 2013), Jørgensen et al. (2002), Mueller et al. (2002), Shirley et al. (2002), Hatchell & van der Tak (2003), Williams et al. (2005), Crimier et al. (2009, 2010a,b), and Kristensen et al. (2012). In an attempt to further fill out the sampling shortage between ~50 and 1000 L_⊙, we also tried to include possible ranges of values from the extensive but resolution-limited survey by Cao et al. (2019) (olive points with large error bars), excluding prestellar sources and sources whose deconvolved sizes were smaller than the Herschel/PACS 160 μm beam, and using the mean and sample standard deviation (1.4 and 0.4, respectively) of the density power law index of the rest of the data to estimate the particle densities at 1000 AU. The resulting distribution of data from Cao et al. (2019) appears to have a slope that is flatter than the aggregate distribution of the other sources, but that is likely an artefact of the sharp lower limit in radius and mass, so we treat these data with caution and exclude them from curve fitting.

Concentrations of data at specific power law indices are artefacts of limited precision in their respective surveys; for example, Williams et al. (2005) rounded their power law indices to the nearest multiple of 0.5. Uncertainties for data other than ours are not shown on the plot, but in all cases, we either adopted the values stated in the literature or assume uncertainties based on contextual information in the original literature (for example if no luminosity uncertainty was given, we adopted the largest percentage uncertainty in the flux) for curve fitting purposes. A machine-readable table of all values used for plotting will be included with the online publication.

For the plots of M_env versus L_bol, r_out versus L_bol, and n(r = 1000 AU) versus L_bol, we minimised the negative log likelihoods for single and broken power law fits (solid and dashed black lines in Fig. 6, respectively), modelled their respective posterior distributions using emcee (Foreman-Mackey et al. 2013), and compared the sum of squared residuals for each fit. For each single power law fit, we solved for:

• the power law exponent,

• the normalisation at L_bol = 1 L_⊙, and

• lnϵ, a nuisance parameter that assumes any residual sources of scatter that are not accounted for by the input uncertainties in the other parameters are constant in L_bol.

The broken power law fit is similar to the single power law fit except that the normalisation is solved for at the power-law break point in L_bol instead of at 1 L_⊙. We somewhat arbitrarily set the break point at L_bol = 60 L_⊙ because the fits failed to converge when that parameter was free, and that value is the approximate centre (in the log scale) of the gap in reliable data for intermediate luminosity protostars. As expected given the amplitude of scatter around both sets of fits, we found that in all cases, residual sums of squares for the broken power law and single power law differed by < 10% (< 2% for M_env and r_out vs. L_bol), smaller thanthe uncertainties for nearly all available data. Thus, we provide the following single power law results as our final fits: $$\text{log}\left(\frac{M_{\text{env}}}{M_{\odot }}\right)={0.79}_{-0.02}^{+0.01}\text{log}\left(\frac{L_{\text{bol}}}{L_{\odot }}\right)+\text{log}({0.30}_{-0.06}^{+0.07})$$(2) $$\text{log}\left(\frac{R_{\text{out}}}{\text{AU}}\right)=0.23\pm 0.02\text{log}\left(\frac{L_{\text{bol}}}{L_{\odot }}\right)+\text{log}({5.4}_{-0.6}^{+0.7}\times {10}^3)$$(3) $$\text{log}\left(\frac{n_{1000\text{AU}}}{{\text{cm}}^{-3}}\right)={0.2}_{-0.2}^{+0.3}\text{log}\left(\frac{L_{\text{bol}}}{L_{\odot }}\right)+5.9\pm \mathrm{0.8.}$$(4)

For Eqs. (2) and (3), the values of lnϵ were $-{0.08}_{-0.11}^{+0.04}$ and − 1.0 ± 0.2, respectively. Because the uncertainties for most values of n(r = 1000 AU) are large, the corresponding nuisance parameter became so small (~ 10⁻¹⁷) that removing it improved the fit. Corner plots of the posterior distributions of all parameters in these fits are available in Appendix D.

For the plot of |p| against L_bol, we followed a similar procedure as for the other three plots to fit a log-linear function in L_bol with a nuisance parameter lnϵ, and to model the posterior distributions. Appendix D also contains the corresponding corner plots for these fitting parameters. We compared this fit (dashed black line in the bottom left panel of Fig. 6) to the mean and standard deviation of |p|, 1.4 and 0.4 respectively (solid and dotted black lines in the bottom left panel of Fig. 6), by performing a Student’s t-Test with the null hypothesis that the slope of the line is indistinguishable from zero. We acknowledge that the uncertainties in |p| for our data and those from the literature are not the same for all data, but they do not appear to depend systematically on L_bol and do not vary much over the full range of L_bol – typical uncertainties are about 0.3. Thus, we believe the Student t-Test should at least be indicative enough to verify the insignificance of the slope suggested by the fit remaining within a standard deviation of the mean over the full L_bol domain. The two-tailed p-value (not to be confused with the density power law index |p|) we calculated from the t-statistic was 0.87. If our use of this test is valid, then it is statistically impossible to distinguish the distribution of our and others’ |p| indices froma Gaussian with a mean of 1.4 and a standard deviation of 0.4. We also note that all values of |p| > 2 have uncertainties large enough that they are consistent with |p| = 2, and that having |p| > 2 is probably not physical, so the apparent upward trend in |p| with L_bol is even less likely to be real than it looks.

4 Discussion

In this section, we start by discussing the results and implications for the eight sources in Cygnus-X, especially their density profiles in the context of their individual shapes and immediate environmental conditions. Then we zoom out to discuss our meta-analysis of envelope parameters in the literature, with our eight sources included. There, we focus on the statistical significance of the scaling relationships derived, that is, whether the relationships for low- and high-mass (luminosity) stars are statistically distinguishable or consistent with one power law per parameter over the full range of sampled protostellar luminosities.

4.1 Profiles of Cygnus-X sources

All of our sources have strong emission peaks at wavelengths of order 10² μm, and less emission towards shorter wavelengths, as shown in Fig. 3. In the protostellar spectral classification scheme of Lada (1987) and Andre et al. (1993), Class 0 and I sources have flux densities increasing with wavelength from NIR to submm wavelengths, whereas Class II and III sources have flux densities decreasing with wavelength over the same range. The decline in flux density towards shorter wavelengths for all of our sources identifies them as Class 0 or I. The original criteria from Andre et al. (1993) for Class 0 sources are not directly readable from an SED, and under the strictest interpretation of the revised criteria in Barsony (1994), all of our sources would be Class I. However, non-detection at λ < 10 μm with Spitzer is very different from a non-detection in the late 1980s to early 1990s, when the most deeply embedded sources detectable in the MIR were up to three orders of magnitude brighter at 100 μm than 10 μm (Lada 1987). Three of our sources are more than three orders of magnitude fainter at λ ≲ 10 μm than at λ ~ 100 μm: N12, N53, and N63. Thus, we argue that these should be at least considered possible Class 0 sources.

The average and standard deviation of p for our eightsources alone are −1.1 and 0.3, respectively. That makes our sub-sample statistically inconsistent with the theoretical slope of −1.5 for isolated sources in free fall, even as all the sources we could find in the literature put together have a mean that is statistically consistent with −1.5 within their dispersion. However, as the three-color Herschel images in Fig. 2 show, our sources are not isolated. There are atleast four potential ways that the modelled envelope profiles could be influenced physically or observationally by their surroundings, any or all of which could be in effect:

Geometric

The shallower slopes are an artefact of spherically averaging over asymmetrical features, like bipolar outflows, core fragments, and spiral arms.

Observational

The density profiles of most PILS-Cygnus sources, like their fluxes, are convolved with those of the filaments in which they reside.

Physical

The envelope may have filaments feeding material into the clump.

FUV

The clumps’ envelopes may be inflated by a moderately strong ISRF.

The geometric factors encapsulate a great number of circumstances that vary from source to source. The PILS-Cygnus sample, though small, is quite diverse in morphology, and asymmetries are readily visible in the continuum images at both submm (Fig. 1) and FIR (Fig. 2) wavelengths. Several follow-up studies have examined members of the MBS07 sample in search of fragmentation, with special attention paid to the clumps that were IR-quiet (weak to non-existent emission in the 1 ≲ λ ≲ 30 μm range, roughly). Studies by Bontemps et al. (2010) and Duarte-Cabral et al. (2013) looked at substructure in N12, N48, N53, and N63 with millimetre continuum and CO, and found all of them had two or more components. In fact, the only sources in the PILS-Cygnus survey not known to contain multiple cores are N38, N51, and N54. S26 looks monolithic, but given its apparent distance and luminosity, it is highly likely to be a protocluster rather than a single protostar; if it is a single protostar, it is likely to have an outflow pointed directly down our line of sight (Trinidad et al. 2003), enhancing its apparent luminosity.

All of the sources have at least one pair of outflows (Skretas et al., in prep.), and some appear to originate from a location offset from the continuum peak – strong indicators that both the geometric and physical factors are at work since outflows are indirect tracers of accretion (see for example Blandford & Payne 1982; Pelletier & Pudritz 1992; Stahler 1994; Hennebelle & Fromang 2008). N30, N51, S8, and S26 are known to be associated with compact HII regions: N30 with W75N(B) (Haschick et al. 1981), especially components VLA 1–3 (Torrelles et al. 1997); N51 with MSX6C G081.7522+00.5906 (Urquhart et al. 2011); S8 with LBN 77.61+01.59 (Lynds 1965); and S26 with RAFGL 2591 VLA 3 (Grasdalen et al. 1983). It is perhaps noteworthy that the only protostellar envelope with |p| > 1.5 is a part of this group. Conversely, given the possibility of association between N48 and 2MASS J20390285+4222001, and the identification of the latter as an HII region by Maud et al. (2015), it is unclear whether or not N48 (|p| ≲ 0.5), should also be in this group. If N48 is associated with an HII region, unless the effects of multiplicity can be disentangled and wholly account for the extremely flat density profile, our limited statistics would no longer support a correlation of p with later-stage protostellar evolution.

It is also worth noting that including a warm inner disk and treating it as part of the spherical envelope is expected to steepen the emission profile, ergo p. This is a likely explanation for the handful of sources with |p| > 2 from Mueller et al. (2002), Hatchell & van der Tak (2003), and Crimier et al. (2010a), although all have sufficiently large uncertainties in p that they are still consistent with |p| = 2. Since we mostly observe flattened profiles among our 8 sources, flattening effects on the density profile (for example core fragmentation) for these sources appear to totally overwhelm any countervailing steepening effects.

The observational and physical factors are especially relevant for N30, N48, N51, and N53, as they are embedded in the larger filamentary structure known as DR21. At the FIR wavelengths covered by Herschel, N53 is mostly likely at least somewhat blended with N54, and N48 is blended with at least two other objects in DR21 South. N30 has a large filament bisecting it from equatorial north-east to south-west, plus at least one additional spur extending towards the east, the latter of which was identified in MBS07 as N31. As discussed in Sect. 2, we tried to account for blending using Gaussian decomposition of the sources; however, because we could not effectively model the embedding filament, the filamentary contribution will affect the density profile, especially at large radii.

We are skeptical that FUV irradiation can have much of a flattening effect on the density profile; the observed impact of the ISRF on the surrounding filaments, where visible, looks more like compression than expansion. In Fig. 2, N12 sits in the head of pillar that appears to protrude into a bubble bounded on the equatorial west side by DR17. The head of this pillar shows a slight increase in opacity that seems indicative of dust and gas being compressed between N12 and radiation from the vicinity of DR17. N12 also has the smallest projected distance from Cygnus-OB2, at about 30 pc compared to ≳ 40 pc for most other sources, caveats listed in Sect. 2 notwithstanding. Yet the shapes of the SED and density profile of N12 are most similar to possibly the most deeply embedded source in our sample, N63. In any case, none of the SED-fitting and Transphere-modelling parameters for N12 or N63 are outliers for this sample, which tentatively suggests that the effects of external irradiation are either degenerate with or relatively insignificant compared to more intrinsic variations (for example mass and core multiplicity).

4.2 Connecting low- and high-mass envelopes

The fits to the aggregate of our data with literature discussed in Sect. 3.3.2 and shown in Fig. 6 suggest that, within the measurement uncertainties, high-mass protostellar envelopes only differ from their low-mass counterparts in scale. Our eight sources sit in the upper intermediate to high mass range and appear to follow established trends, where any are visible: envelope mass, outer envelope radius, and gas number density at 1000 AU all seem to increase as power laws in L_bol. Only the gas density power law index appears to be more or less independent of luminosity, with an uncertainty-weighted mean value of 1.4 and sample standard deviation of 0.4 for all sources put together. This mean and standard deviation is consistent with the theoretical |p| = 1.5 for sources in free-fall. Recall, however, that our eight sources alone had a mean |p| and σ_p of 1.1 ± 0.3, which are marginally inconsistent with the theoretical |p| = 1.5. This means that, while our 8 sources do help fill in the intermediate luminosity gap, they are not a very representative sample in p, at least given current observational limitations.

If the relative importance of, say, turbulence or magnetic fields as compared to gravity differed between protostellar mass regimes, one would expect that to show up as a systematic difference in typical values of |p|, which should propagate to the slopes of n(r = 1000 AU) and/or M_env versus L_bol. A cursory glance at Fig. 6, without the benefit of visible error bars on every point, does seem to vaguely suggest a transition or break at L_bol between a few tens and a few hundred L_⊙ for all relationships except the trend in r_out with L_bol. However, the uncertainties are so large that only a first order trend can be established with any confidence; any apparent change of trend with increasing L_bol is simply not statistically significant. More data are needed on new sources, and measurements for existing sources need to be refined further.

Fits to M_env, r_out, and n(r = 1000 AU) as functions of luminosity seem primarily driven by just the three most extensive data sets (after excluding all older measurements of the same objects): Mueller et al. (2002) and Williams et al. (2005) at the high-mass, high-luminosity end, and Kristensen et al. (2012) at the low-mass, low-luminosity end. At the high-luminosity end, much of the data from Williams et al. (2005) stands noticeably apart from most other data sets in the same luminosity regime. For r_out and n(r = 1000 AU), the uncertainties on these data are so large that the fitting routine does not give them much weight, but for the M_env–L_bol relationship, these data seem to skew the fit towards a steeper slope than the rest of the data appear to warrant. We repeated the fits of single and broken power laws to this relationship excluding the data from Williams et al. (2005) to see if the two functions then became statistically distinguishable. They did not, likely because some data from Hatchell & van der Tak (2003) also have large scatter towards higher M_env, even though they are not outliers in the trends of r_out and n(r = 1000 AU) with L_bol.

For the upper two panels of Fig. 6, we caution that envelope mass depends strongly on outer envelope radius, and methods for measuring the outer envelope radii varied from study to study. For most of the low-mass sources, particularly from Jørgensen et al. (2002) and Kristensen et al. (2012), r_out was defined as the radius where temperature dropped to 10 K. For high-mass sources, radii were typically measured in submm or short mm-wave images, either as an average radius of a contour above a background threshold or as some multiple of a Gaussian FWHM or σ; often the definition of r_out was unclear or not defined in terms of Gaussian parameters. For most high mass sources in general, it does not make much sense to try to calculate a 10 K radius because high-mass star-forming regions tend to be highly irradiated environments where only the heavily shielded centres of prestellar cores get that cold. The majority of our sources in particular are within 300 pc of Cygnus-OB2 or other OB clusters within Cygnus-X, so the ambient medium should be warmer, and indeed, the map of dust temperatures around DR21 by Hennemann et al. (2012) shows dust temperatures mostly between 15 and 22 K. These systematic differences in how r_out is measured indifferent mass and luminosity regimes may be at least partially responsible for how the mass-luminosity and radius-luminosity trends visually appear slightly flatter for high-mass sources than low-mass sources. As of now, this visual change of trend between luminosity regimes is not statistically significant, but if more data make it so, future astronomers will have to start trying both methods on the same sources wherever feasible to quantify the difference and see if there is a typical fractional discrepancy that can be applied as a correction to older data.

Additionally, the methods typically used for measuring the radii of high-mass sources yield somewhat wavelength dependent results, with less emission being recovered at mm wavelengths than closer to the wavelength of the peak of the dust SED. Furthermore, ground-based observations are limited by sky background, which will necessarily decrease the measured radius. For our sources, we again caution that the method we used to fit the envelope radii yields large and sometimes poorly-defined uncertainties due to the degeneracy with the density power law index and – in the case of N48 and S8 – complex, radially asymmetric morphology (see the discussion of Table 4).

5 Conclusions

In this study, we have deconvolved and fit greybody SEDs to archival IR data from 1.2 to 850 μm on eight intermediate- to high-mass Class 0-I protostellar sources. We used the derived luminosities and total mass estimates from SED fitting, and measurements of the sources’ sizes at 450 μm, to perform 1D radiative transfer modelling using Transphere, with which we derived our sources’ envelope masses, radial density structures, and temperature profiles. Finally, we combined these results with similar data in the literature for as many distinct sources as we could find, and fit the aggregated data to determine if low- and high-mass stars are statistically distinguishable. Our results are as follows:

1. Most of our eight sources have envelope masses between a few tens and a couple hundred M_⊙, and luminosities of several hundred to several thousand L_⊙, placing them in an intermediate mass and luminosity range that is otherwise poorly sampled in the literature.

2. The power law indices of our eight sources’ radial density profiles cluster around |p| = 1.1 with a dispersion of 0.3, slightly shallower than predicted for an isolated, gravitationally collapsing source (|p| = 1.5). We attribute the flattening to some combination of active accretion, spherical averaging of non-spherical features (for example disks and outflow cavities), and convolution with surrounding cloud structures. We briefly discuss a fourth possibility, inflation by a strong IRSF, but find it unconvincing.

3. When our sources are considered with other sources in the literature, we find the scatter is too large in envelope mass, outer radius, density power-law index, and gas density at 1000 AU to definitively establish if trends differ for low and high mass-sources. Low- and high-mass sources appear to differ significantly only in scale.

4. We establish to first-order that log(M_env)$\propto {0.79}_{-0.02}^{+0.01}$log(L_bol), log(r_out) ∝ 0.23 ± −0.02log(L_bol), and log(n(r = 1000 AU))$\propto {0.2}_{-0.2}^{+0.3}$log(L_bol) over seven decades of L_bol, with a scatter that is consistently about an order of magnitude above and below these fits.

5. For the aggregate of sources in the literature plus our sources, we determine that the distribution of density power law indices isconsistent with a constant |p| = 1.4 ± 0.4. With this much scatter, if low- and high-mass protostars get different relative amounts of turbulent or magnetic field support against gravity, we do not have the precision to detect the difference at a statistically significant level. It also implies that our eight sources were not a representative sample in this parameter.

Our inability to distinguish trends in the envelope parameters of low- and high-mass protostars (M_env, r_out, n(r = 1000 AU), p) clearly warrants further investigation, especially more and larger surveys that focus on protostars in the L_bol range of 50–500 L_⊙. Nonetheless,the results presented here suggest that the differences between low- and high-mass sources do not arise from variations inlarge-scale physical properties.

Acknowledgements

The research of R.L.P. and L.E.K. is supported by a research grant (19127) from VILLUM FONDEN. J.K.J. acknowledges support from the Independent Research Fund Denmark (grant number DFF0135-00123B). This research has made use of the NASA/IPAC Infrared Science Archive, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration (NASA). This research was made possible with the use of the Submillimeter Array. The Submillimeter Array is a joint project between the Smithsonian Astrophysical Observatory and the Academia Sinica Institute of Astronomy and Astrophysics and is funded by the Smithsonian Institution and the Academia Sinica. This research used archival data from Herschel, an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. This work is based in part on archival data obtained with the Spitzer Space Telescope, which was operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. Support for this work was provided by an award issued by JPL/Caltech. This research also made use of data products from the Midcourse Space Experiment. Processing of the data was funded by the Ballistic Missile Defense Organisation with additional support from NASA Office of Space Science. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the NASA and the National Science Foundation. This research also used the facilities of the Canadian Astronomy Data Centre (CADC) operated by the National Research Council of Canada with the support of the Canadian Space Agency. This research made use of Astropy (http://www.astropy.org), a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018).

Appendix A Measured and deconvolved continuum fluxes

Appendix B SEDs fits

Figures B.1-B.7show the data (black filled and hollow cir- cles with error bars), greybody SED fits (dashed red lines), and Transphere SED fits (solid blue lines) for all other objects in the PILS-Cygnus survey that could be modelled, besides N30. The parameters, processes, and caveats involved in each fitting routine are discussed in §2. Hollow circles indicate data that were excluded from the fit with Transphere, but were fit if pos- sible with our initial greybody fitting routine.

Fig. B.1 Same as Fig. 4, but forN12. It is uncertain whether the PACS-70 μm data are part of the same temperature component as the longer-wavelength data; the fits both with and without it are within the uncertainties, but not including the 70 μm data yields a more realistic β value.

Fig. B.2 Same as Fig. B.1, but for N48. N48 is located in an especially dense part of the DR21 ridge and is blended with two other clumps identified in MBS07 at λ > 250 μm. We advise caution in using the derived SED parameters for this object in any precision analysis.

Fig. B.3 Same as Fig. B.1, but for N51.

Fig. B.4 Same as Fig. B.1, but for N53.

Fig. B.5 Same as Fig. B.1, but for N63.

Fig. B.6 Same as Fig. B.1, but for S8. The distance to S8 is highly uncertain and it has a spiral arm feature thatcould only be accounted for to first order, so SED fitting parameters should be used with caution.

Fig. B.7 Same as Fig. B.1, but for S26. We note the change in y-axis scaling.

Appendix C Density power-law index fits, and temperature and density profiles

Figures C.1-C.7 are the equivalent of Fig. 5 for all other sources besides N30.

Fig. C.1 Same as Fig. 5, but for N12 (using the fit that excludes the 70 μm data).

Fig. C.2 Same as Fig. 5, but for N48.

Fig. C.3 Same as Fig. 5, but for N51.

Fig. C.4 Same as Fig. 5, but for N53.

Fig. C.5 Same as Fig. 5, but for N63.

Fig. C.6 Same as Fig. 5, but for S8.

Fig. C.7 Same as Fig. 5, but for S26. We note the change of x-axis scaling.

Appendix D Optimisation and posterior distributions of power-law fits to envelope parameters

In Fig. 6, we fitted power-law trends in M_env, r_out, |p|, and n(r = 1000 AU) with L_bol to see if trends for high- and low-mass sources were statistically distinguishable from a single trend for all sources. The details of the process are discussed in §3.3.2. The following figures show the posterior distributions of the pa- rameters of each of the pair of fits to each relationship, except for the fit to |p| as a function ofL_bol. As we discussed in the main text, we did not find a significant correlation between p and L_bol.

Fig. D.1 (a) Posterior distributions of the power-law exponent (ξ), normalisation (Mo), and residual variance parameter (lnϵM) for the single power-law fit to Menv vs. Lbol, shown asthe solid black line in the top left panel of Fig. 6. (b) Posterior distributions of the power-law exponents (ξ1 and ξ2), normalisation (Mbk), and residual variance parameter (lnϵM) for the broken power-law fit to Menv vs. Lbol, shown as the dashed black line in the top left panel of Fig. 6.

Fig. D.2 (a) Posterior distributions of the power-law exponent (α), normalisation (Ro), and residual variance parameter (lnϵR) for the single power-law fit to rout vs. Lbol, shown asthe solid black line in the top right panel of Fig. 6. (b) Posterior distributions of the power-law exponents (α1 and α2), normalisation (Rbk), and residual variance parameter (lnϵR) for the broken power-law fit to rout vs. Lbol, shown as the dashed black line in the top right panel of Fig. 6.

Fig. D.3 (a) Posterior distributions of the power-law exponent (η) and normalisation (no) for the single power-law fit to n(r = 1000 AU) vs. Lbol, shown as the solid black line the bottom right panel of in Fig. 6. (b) Posterior distributions of the power-law exponents (η1 and η2), normalisation (nbk), and residual variance parameter (lnϵn) for the broken power-law fit to n(r = 1000 AU) vs. Lbol, shown as the dashed black line in the bottom right panel of Fig. 6.

References

All Tables

All Figures

Fig. 1 SMA 868 μm (345 GHz) images of all 10 PILS-Cygnus sources (N53 and N54 are both featured in the middle right panel) contoured in red at 10, 40, and 70% of the maximum flux. The cyan ellipses in the upper right corners show the FWHMa of both axes and positionangle of the synthesised SMA beam, and the line segments alone the bottoms of each panel show how long 10 000 AU and its corresponding angular length appears at the adopted distance to each object. We note the change of scale for S26. All images are normalised to their peak flux densities, which are given in Table 1.

In the text

	Fig. 2 Herschel three-colour cut-outs showing the PILS-Cygnus sources in their natal environments at 70, 160, and 250 μm, with the DR objects labelled. Herschel beams are enlarged by two-fold for visibility.
In the text

Fig. 3 Full SEDs for all sources, which were fitted, interpolated, and integrated by trapezoidal sums to calculate the bolometric luminosities. The curves running through data at λ ≥ 70 μm are the FIR/submm component of the SED fits discussed in Sect. 2.3, while the rest of the data has been log-linearly interpolated in lieu of integrating over the less-physical warmer SED components. For N30 and N51, no SED components passed through the 70 μm data point, so the SED fit was only used for λ ≥ 160 μm while the 70 μm data were linearly interpolated like the shorter wavelength data. The light red filled curves show the margins of uncertainty as constrained by the flux errors, and were likewise integrated to estimate the total uncertainty in luminosity.

In the text

	Fig. 4 SEDs of N30. The dashed red line is the modified greybody fit to the data, shown as black and grey markers with error bars. The solid blue line is the best fit Transphere SED to only the data shown in black. The grey points at λ < 70 μm were excluded from the fit because the SED at these wavelengths likely traces a warm inner disk where the assumption of spherical symmetry breaks down.
In the text

	Fig. 5 A. Radial emission profile comparison between the data (red dashed line) and the Transphere model (blue solid line), both normalised by their respective maxima. B. Resulting model temperature (red solid line and red axis) and density (blue dash-dotted line and blue axis) profiles.
In the text

Fig. 6 Comparison of PILS-Cygnus Sources (large black circles with error bars) to objects in the literature from van der Tak et al. (2000) (purple x’s), Jørgensen et al. (2002) (magenta upward-pointing triangles), Mueller et al. (2002) (thin blue diamonds), Shirley et al. (2002) (maroon stars), Hatchell & van der Tak (2003) (green squares), Williams et al. (2005) (small orange circles), Crimier et al. (2009, 2010a,b) (wide navy diamonds), Kristensen et al. (2012) (cyan downward-pointing triangles), van der Tak et al. (2013) (red plus signs), and Cao et al. (2019) (pale gold error bars). Solid and dashed black lines represent single and broken power law fits to the data, respectively.

In the text

	Fig. B.1 Same as Fig. 4, but forN12. It is uncertain whether the PACS-70 μm data are part of the same temperature component as the longer-wavelength data; the fits both with and without it are within the uncertainties, but not including the 70 μm data yields a more realistic β value.
In the text

	Fig. B.2 Same as Fig. B.1, but for N48. N48 is located in an especially dense part of the DR21 ridge and is blended with two other clumps identified in MBS07 at λ > 250 μm. We advise caution in using the derived SED parameters for this object in any precision analysis.
In the text

	Fig. B.3 Same as Fig. B.1, but for N51.
In the text

	Fig. B.4 Same as Fig. B.1, but for N53.
In the text

	Fig. B.5 Same as Fig. B.1, but for N63.
In the text

	Fig. B.6 Same as Fig. B.1, but for S8. The distance to S8 is highly uncertain and it has a spiral arm feature thatcould only be accounted for to first order, so SED fitting parameters should be used with caution.
In the text

	Fig. B.7 Same as Fig. B.1, but for S26. We note the change in y-axis scaling.
In the text

	Fig. C.1 Same as Fig. 5, but for N12 (using the fit that excludes the 70 μm data).
In the text

	Fig. C.2 Same as Fig. 5, but for N48.
In the text

	Fig. C.3 Same as Fig. 5, but for N51.
In the text

	Fig. C.4 Same as Fig. 5, but for N53.
In the text

	Fig. C.5 Same as Fig. 5, but for N63.
In the text

	Fig. C.6 Same as Fig. 5, but for S8.
In the text

	Fig. C.7 Same as Fig. 5, but for S26. We note the change of x-axis scaling.
In the text

Fig. D.1 (a) Posterior distributions of the power-law exponent (ξ), normalisation (Mo), and residual variance parameter (lnϵM) for the single power-law fit to Menv vs. Lbol, shown asthe solid black line in the top left panel of Fig. 6. (b) Posterior distributions of the power-law exponents (ξ1 and ξ2), normalisation (Mbk), and residual variance parameter (lnϵM) for the broken power-law fit to Menv vs. Lbol, shown as the dashed black line in the top left panel of Fig. 6.

In the text

Fig. D.2 (a) Posterior distributions of the power-law exponent (α), normalisation (Ro), and residual variance parameter (lnϵR) for the single power-law fit to rout vs. Lbol, shown asthe solid black line in the top right panel of Fig. 6. (b) Posterior distributions of the power-law exponents (α1 and α2), normalisation (Rbk), and residual variance parameter (lnϵR) for the broken power-law fit to rout vs. Lbol, shown as the dashed black line in the top right panel of Fig. 6.

In the text

Fig. D.3 (a) Posterior distributions of the power-law exponent (η) and normalisation (no) for the single power-law fit to n(r = 1000 AU) vs. Lbol, shown as the solid black line the bottom right panel of in Fig. 6. (b) Posterior distributions of the power-law exponents (η1 and η2), normalisation (nbk), and residual variance parameter (lnϵn) for the broken power-law fit to n(r = 1000 AU) vs. Lbol, shown as the dashed black line in the bottom right panel of Fig. 6.

In the text