11institutetext: Niels Bohr Institute, Centre for Star & Planet Formation,
University of Copenhagen, Øster Voldgade 5-7, 1350 Copenhagen K, Denmark
11email<EMAIL_ADDRESS>
# Linking high- and low-mass star formation

Observation-based continuum modelling and physical conditions

R. L. Pitts, ORCiD: 0000-0002-7937-493111 L. E. Kristensen 11 J. K. Jørgensen
11 S. J. van der Walt

(Received September 28, 2021; accepted October 20, 2021)

###### Abstract

Context. Astronomers have yet to establish whether high-mass protostars form
from high-mass prestellar cores, similar to their lower-mass counterparts, or
from lower-mass fragments at the heart of a pre-protostellar cluster
undergoing large-scale collapse. Part of the uncertainty is due to a shortage
of envelope structure data on protostars of a few tens of solar masses, where
we expect to see a transition from intermediate-mass star formation to the
high-mass process.

Aims. We sought to derive the masses, luminosities, and envelope density
profiles for eight sources in Cygnus-X, whose mass estimates in the literature
placed them in the sampling gap. Combining these sources with similarly
evolved sources in the literature enabled us to perform a meta-analysis of
protostellar envelope parameters over six decades in source luminosity.

Methods. We performed spectral energy distribution (SED) fitting on archival
broadband photometric continuum data from 1.2 to 850 $\mu$mto derive
bolometric luminosities for our eight sources plus initial mass and radius
estimates for modelling density and temperature profiles with the radiative-
transfer package Transphere.

Results. The envelope masses, densities at 1000 AU, outer envelope radii, and
density power law indices as functions of bolometric luminosity all follow
established trends in the literature spanning six decades in luminosity. Most
of our sources occupy an intermediate to moderately high range of masses and
luminosities, which helps to more firmly establish the continuity between low-
and high-mass star formation mechanisms. Our density power law indices are
consistent with observed values in the literature, which show no discernible
trends with luminosity, and have a mean $p=-1.4\pm 0.4$. However, our sub-
sample, with a mean power law index of $-1.1\pm 0.3$, is slightly flatter than
would be expected for spherical envelopes in free fall ($p=-1.5$).

Conclusions. We attribute flattened density profiles for our eight sources to
one or more of the following: ongoing accretion from their natal filaments,
convolution of sources with neighbours or the larger filament, spherical
averaging of asymmetric features (for example fragments), or inflation of the
envelope by a moderate far-ultraviolet (FUV) field. Finally, we show that the
trends in all of the envelope parameters for high-mass protostars are
statistically indistinguishable from trends in the same variables for low-mass
protostars.

###### Key Words.:

Stars: formation – ISM: dust, extinction – (Stars:) circumstellar matter –
Submillimeter: stars – Submillimeter: ISM – Stars: protostars

##  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\gtrsim 8$ $M_{\odot}$—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_{\odot}<M<8$ $M_{\odot}$), and low-mass ($M\lesssim 1.5$
$M_{\odot}$) depending on whether or not the stars are small and cool enough
to sustain a convection in their outer layers (Hansen et al., 2004, §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 ($\sim 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 $\rho(r)\propto 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_{\odot}$ ($\sim 50$–1000 $L_{\odot}$). 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\sim 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 §2, we discuss how we obtained and pre-processed our continuum data to get
fluxes for SED fitting; how we used SED fitting to estimate parameters, like
envelope mass and luminosity, for bench-marking 1D envelope models; and how we
then used those SED fits and derived parameters to model the temperature and
density structures of our sources. In §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 §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 §5

##  2 Methods

As discussed in §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, H2, 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 H2 (as explained in §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 is much 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, and is 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.

###  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 H2O $2_{02}-1_{11}$
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
$\mu$m) in both compact and extended array configurations for a combined
spatial resolution of $0\aas@@fstack{\prime\prime}7$ and sensitivity of 0.15
Jy km s-1 (see van der Walt et al. 2021, to appear in A&A, for detailed
description of data reduction). Figure 1 shows the 868 $\mu$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 $\mu$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
$\mu$m to construct a believable SED.

Figure 1: SMA 868 $\mu$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 position angle of the synthesised SMA beam,
and the line segments alone the bottoms of each panel show how long 10000 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 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.

Table 1: Table of PILS-Cygnus source coordinates

Source(a) | RA(b) | Dec(b) | Associated Objects | D⊙ (pc)(c) | D${}_{CygOB2}^{proj}$ (pc)(d) | Fragments(e) |  $S^{max}_{868\mu m}$ (Jy/beam)(f)
---|---|---|---|---|---|---|---
N12 | 20:36:57.68 | 42:11:30.8 | J20365781+4211303 |  $1400\pm 100$(g) | $31\pm 2$ | 2+ | 0.50
N30 | 20:38:36.45 | 42:37:33.8 | W75N(B) | $1300\pm 70$ | $43^{+2}_{-3}$ | 4 | 1.29
N38(h) | 20:38:59.26 | 42:22:28.6 | DR 21(OH) W | $1500^{+80}_{-70}$ | $47\pm 2$ | 1 | 0.62
N48 | 20:39:01.46 | 42:22:05.9 | DR 21(OH) S | $1500^{+80}_{-70}$ | $47\pm 2$ | 5+ | 0.25
N51 | 20:39:01.97 | 42:24:59.2 | G081.7522+00.5906 | $1500^{+80}_{-70}$ | $48^{+2}_{-3}$ | 1+ | 0.42
N53 | 20:39:03.13 | 42:25:52.5 | W75S FIR3 | $1500^{+80}_{-70}$ | $48\pm 2$ | 2+ | 0.41
N54(h) | 20:39:04.03 | 42:25:41.1 |  | $1500^{+80}_{-70}$ | $48\pm 2$ | 1 | 0.27∗
N63 | 20:40:05.50 | 41:32:12.6 | BGPS G081.174-00.100 |  $1400\pm 100$(g) | $42^{+3}_{-2}$ | 2 | 1.66
S8 | 20:20:39.29 | 39:37:54.2 | Mol 121 |  $1400^{+2600}_{-1100}$(gi) |  $>80$(i) | 2+ | 0.30
S26 | 20:29:24.87 | 40:11:19.3 | RAFGL 2591 VLA 3 | $3330\pm 110$ |  $\cdots$(j) | 1+ | 0.60

  * a

Object names from MBS07.

  * b

Coordinates corrected using _Spitzer_ images overlaid with SMA contours and
ellipses defined by FWHMa from MBS07.

  * c

Heliocentric distances and uncertainties are from Rygl et al. 2012 and
references therein unless otherwise noted.

  * d

Projected from Cyg OB2 calculated using D⊙.

  * e

Number of core components visible in SMA images (868 $\mu$m,
$0\aas@@fstack{\prime\prime}7$ resolution). A plus (+) indicates that weak
and/or filamentary emission suggest there may be additional fragments.

  * f

Maximum 868 $\mu$m flux density among all core fragments. (∗Peak flux density
of N54 measured from image centred on N54.)

  * g

Rygl et al. 2012 did not explicitly measure the parallaxes at the locations of
these sources, so we adopted the average and uncertainty over Cygnus X.

  * h

SEDs could not be reconstructed for these sources. See text for details.

  * i

The distance to S8 is poorly constrained; the given uncertainties reflect the
full range of distance estimates in the literature (see review by Varricatt et
al. 2010, appendix A37), and the smallest possible projected distance can by
calculated by assuming both S8 and Cyg OB2 are at the latter’s minimum
distance from us ($\sim 1.3$ kpc).

  * j

The projected distance between S26 and Cyg OB2 is negligible compared to their
line-of-sight separation, $\gtrsim 1.57$ kpc.

|
---|---
|
Figure 2: _Herschel_ three-colour cut-outs showing the PILS-Cygnus sources in
their natal environments at 70, 160, and 250 $\mu$m, with the DR objects
labelled. _Herschel_ beams are enlarged by two-fold for visibility.

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 5.

Table 2: Instrument and filter specifications for all filters used.

Band |  $\lambda$ ($\mu$m) | Resolution | Pixel Size
---|---|---|---
2MASS-J | 1.235 | 4′′ | 2′′
2MASS-H | 1.662 | 4′′ | 2′′
2MASS-Ks | 2.159 | 4′′ | 2′′
IRAC-1 | 3.6 | 1$\aas@@fstack{\prime\prime}$66 | 0$\aas@@fstack{\prime\prime}$6(a)
IRAC-2 | 4.5 | 1$\aas@@fstack{\prime\prime}$72 | 0$\aas@@fstack{\prime\prime}$6(a)
IRAC-3 | 5.8 | 1$\aas@@fstack{\prime\prime}$88 | 0$\aas@@fstack{\prime\prime}$6(a)
IRAC-4 | 8.0 | 1$\aas@@fstack{\prime\prime}$98 | 0$\aas@@fstack{\prime\prime}$6(a)
MSX-A | 8.28 | 18$\aas@@fstack{\prime\prime}$3 | 6′′
MSX-C | 12.13 | 18$\aas@@fstack{\prime\prime}$3 | 6′′
MSX-D | 14.65 | 18$\aas@@fstack{\prime\prime}$3 | 6′′
MSX-E | 21.4 | 18$\aas@@fstack{\prime\prime}$3 | 6′′
MIPS-24 | 24 | 6′′ | 2$\aas@@fstack{\prime\prime}$45
IRAS-25 | 25 | 0$\aas@@fstack{\prime}$5 |  $\cdots$(b)
IRAS-60 | 60 | 1′ |  $\cdots$(b)
PACS-B | 70 | 5$\aas@@fstack{\prime\prime}$86 $\times$ 12$\aas@@fstack{\prime\prime}$16 | 3$\aas@@fstack{\prime\prime}$2
PACS-R | 160 | 11$\aas@@fstack{\prime\prime}$64 $\times$ 15$\aas@@fstack{\prime\prime}$65 | 3$\aas@@fstack{\prime\prime}$2
SPIRE-B | 250 | 17$\aas@@fstack{\prime\prime}$9 | 6′′
SPIRE-G | 350 | 24$\aas@@fstack{\prime\prime}$2 | 10′′
SCUBA-B(c) | 450 | 11′′$/$ 7$\aas@@fstack{\prime\prime}$9 | 3′′$/$ 2′′
SPIRE-R | 500 | 35$\aas@@fstack{\prime\prime}$4 | 14′′
SCUBA-R(c) | 850 | 19$\aas@@fstack{\prime\prime}$5 $/$ 13′′ | 6′′$/$ 3′′

  * a

Images of S8 had 0$\aas@@fstack{\prime\prime}$86 pixels.

  * b

Catalogued fluxes were used as upper limits.

  * c

Numbers to the left of the slash are for SCUBA Legacy data, which we used for
photometry. Numbers to the right are for SCUBA-2 data, which we used for
density profile fitting but not photometry because too much flux was lost to
secondary beam components.

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, so we 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$\sigma$ 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 $\mu$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 $\mu$m data, but our data do not allow confirmation.

###  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_{\textrm{dust}}$) and molecular hydrogen column
density ($N_{\textrm{H}_{2}}$), interpolates on grids of colour correction
factors111Colour correction functionality is only available for
$1.5\leq\beta\leq 2.0$ for _Herschel_ observations. Corrected fluxes are
typically different from the uncorrected values by less than the typical 10%
margin of uncertainty. Not colour correcting sources with temperatures ¿10 K
systematically depresses the SPIRE fluxes by $\lesssim 10\%$, but this is more
than counteracted by the tendency of our sources to be increasingly
contaminated at longer wavelengths., 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_{\mathrm{dust}})\left(1-\mathrm{exp}\left(-\left(\frac{\nu}{\nu_{0}}\right)^{\beta}\frac{\kappa_{0}N_{\mathrm{H}_{2}}}{\gamma}\mu
m_{\mathrm{H}}\right)\right)$ (1)

where $B_{\nu}$($T_{\textrm{dust}}$) is the Planck function, $\beta$ is the
dust emissivity index (typically with $1\lesssim\beta\lesssim 2$),
$\kappa_{0}$ is the dust opacity at reference frequency $\nu_{0}$, $\gamma$ is
the gas-to-dust ratio (assumed to be 100), $\mu$ is the mean molecular weight
per hydrogen molecule (2.8; see Appendix A of Kauffmann et al. 2008), and
$m_{\mathrm{H}}$ is the mass of a hydrogen atom. For each object, we took the
best-fitting value of $\kappa_{0}$ at 350 $\mu$m ($\nu_{0}=856$ GHz) for
grains with thin ice mantles from Ossenkopf & Henning 1994. Note that Eq. 1
does not assume optically thin dust emission, because for most PILS-Cygnus
sources, optical depth at 70 $\mu$m is not much less than 1.

Ordinarily $T_{\textrm{dust}}$ is somewhat degenerate with $\beta$ and
$N_{\textrm{H}_{2}}$ is degenerate with both $\gamma$ and $\kappa_{0}$, so
only $T_{\textrm{dust}}$ and $N_{\textrm{H}_{2}}$ were initially left as free
parameters. However, we found that we could get $\beta$ to converge by
disabling colour correction and fitting only data at wavelengths longer than
either 70 or 160 $\mu$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 $\lambda\gtrsim 70$ $\mu$m to derive
$\beta$, then fixed $\beta$ 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 $\beta$ 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 $\beta$ 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 $\mu$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 $\ll 1\%$ of the
total mass.

We also note that some sources—N12, N51, N53, and N63—have had emission at 70
$\mu$m excluded or assigned to a different temperature component of the SED
fit. In these cases, the 70 $\mu$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 $\mu$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 $\mu$m lie within the expected range of
values possible for a single cold greybody SED given the uncertainty in
$\beta$.

As usual, $\gamma$ and $\kappa_{0}$ 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_{\textrm{dust}}$ and
$N_{\textrm{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$\sigma$ 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
$\mu$m silicate absorption feature or the prominent emission lines attributed
to polycyclic aromatic hydrocarbons (PAHs) between 3 and 20 $\mu$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
$\lambda<50$ $\mu$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 $\lambda<50$ $\mu$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_{\textrm{TP}}$ controls the amplitude of the SED and the
position of the peak in wavelength-space via the usual $L_{\textrm{TP}}\propto
T^{4}$ and Wien’s Displacement Law. Increasing $L_{\textrm{TP}}$ also boosts
the temperature profile at all radii. Instead of $L_{\textrm{bol}}$, we used
only the luminosity calculated from _Herschel_ data and longer wavelengths as
initial values of $L_{\textrm{TP}}$ because of the aforementioned breakdown of
spherical symmetry at shorter wavelengths.

Stellar surface temperature $T_{\star}$ 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 10000 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_{\star}$ does slightly steepen the radial temperature profile at
small ($r\lesssim 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_{\star}=5000$ K for all PILS-Cygnus sources.

Envelope mass $M_{\textrm{env}}$ has a rather complex effect on the SED.
Increasing $M_{\textrm{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_{\textrm{env}}$, all else
equal, reaching higher temperatures at small radii and lower temperatures at
large radii (greater than a few $\times 10^{3}$ 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 $\times 10^{3}$ 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_{\textrm{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_{\textrm{in}}$ within about 20 AU of this radius, and contained one or more
discontinuities for $r_{\textrm{in}}$ between about 80 and 200 AU. We fixed
$r_{\textrm{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_{\textrm{out}}$ has a similar effect on the SED as
$M_{\textrm{env}}$, but not identical: increasing $r_{\textrm{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_{\textrm{out}}$ has
the effect of changing the concavity in log-log space. We initialised
$r_{\textrm{out}}$ to 20000 AU and fit it iteratively along with $p$.

Interstellar Radiation Field (ISRF) strength affects the height of the SED at
$\lambda\lesssim 50$ $\mu$m and the temperature profile in the outermost few
$\times 10^{3}$ AU (out of 1 to a few $\times 10^{4}$ AU). Increasing the ISRF
strength can boost the SED at $\lambda\lesssim 50$ $\mu$m, but never enough to
explain the observed fluxes at those wavelengths without convergence errors.
We fixed the ISRF strength at $1\,G_{0}$ with $G_{0}$ in Habing 1968 units.

To find the optimal envelope masses, luminosities, and density profiles of the
PILS-Cygnus sources, we took an iterative $\chi^{2}$-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_{\textrm{out}}$).
Second, we minimised $\chi^{2}$ between the model SEDs and data with
$\lambda>50$ $\mu$m by interpolating on a grid of $\chi^{2}$ as a function of
luminosity ($L_{\textrm{TP}}$) and envelope mass ($M_{\textrm{env}}$) only,
given a fixed $p$. Third, we created smaller grids of model intensity images
for each object, given $L_{\textrm{TP}}$ and $M_{\textrm{env}}$, for all
reasonable values of $p$ and $r_{\textrm{out}}$. Fourth, we interpolated on a
grid of $\chi^{2}$ as a function of $p$ and $r_{\textrm{out}}$ to minimise
$\chi^{2}$ between the model images and close-cropped cutouts from the
observational data. Finally, we took the resulting $p$ and $r_{\textrm{out}}$,
and repeated the second through fourth steps until $L_{\textrm{TP}}$,
$M_{\textrm{env}}$, $p$, and $r_{\textrm{out}}$ converged. Convergence
typically occurred by the third iteration.

Data at $\lambda\lesssim 30$ $\mu$m were excluded from the fit because these
data mostly trace warm ($T_{\textrm{dust}}$$\gtrsim 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 ($\lesssim$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 $\lambda\lesssim 30$ $\mu$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_{\textrm{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 §3.1. We compare the
results for each type of analysis, and discuss potential sources of
discrepancies in §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 (§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 (§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$\sigma$ mean radius from
MBS07 corrected for distance measurements from Rygl et al. 2012
($r_{1.2\textrm{mm}}$), dust temperature of the cold ($T_{\textrm{dust}}<60$
K) temperature component, H2 column density in the cold dust component
($N[\textrm{H}_{2}]$), the dust emissivity index $\beta$, dust temperature and
(H+)H2 column density of any additional fitting components
($T_{\textrm{dust},\;i>1}$ and $N_{i>1}[\textrm{(H+)H}_{2}]$, respectively)
listed in order of increasing $T_{\textrm{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
$\mu$m.

Table 3: Modified blackbody fitting parameters for PILS-Cygnus sources.

Source |  $r_{1.2\textrm{mm}}$ (AU) |  $T_{\textrm{dust}}$ (K) |  $N[\textrm{H}_{2}]$ (cm-2) | $\beta$ |  $T_{\textrm{dust},\;i>1}$ (K)(a) |  $N_{i>1}[\textrm{(H+)H}_{2}]$ (cm-2)(a)
---|---|---|---|---|---|---
N12(a) | $1.9\times 10^{4}$ | $22\pm 2$ | $9.5\pm 2.7\times 10^{22}$ | $1.2\pm 0.2$ | $\cdots$ | $\cdots$
N30 | $2.3\times 10^{4}$ | $37\pm 1$ | $7.8\pm 0.4\times 10^{22}$ | $1.3\pm 0.1$ | $396.4\pm 0.1$ | $5.0\pm 0.5\times 10^{15}$
N48 | $2.6\times 10^{4}$ | $40\pm 10$ | $7\pm 1\times 10^{22}$ | $0.6\pm 0.4$ | $\cdots$ | $\cdots$
N51 | $2.6\times 10^{4}$ | $22\pm 1$ | $9\pm 2\times 10^{22}$ | $1.7\pm 0.6$ |  $68\pm 3$, |  $<6.3\times 10^{20}$,
|  |  |  |  | $350\pm 20$ | $<10^{16}$
N53 | $2.2\times 10^{4}$ | $18\pm 5$ | $1.2\pm 1.1\times 10^{23}$ | $1.1\pm 0.7$ | $\cdots$ | $\cdots$
N63 | $1.2\times 10^{4}$ | $25\pm 4$ | $1.2\pm 0.4\times 10^{23}$ | $0.7\pm 0.3$ | $\cdots$ | $\cdots$
S8 | $2.9\times 10^{4}$ | $38\pm 1$ | $3.0\pm 0.2\times 10^{22}$ | $1.1\pm 0.2$ |  $98\pm 5$, |  $9\pm 2\times 10^{19}$,
|  |  |  |  |  $360\pm 20$, |  $1.0\pm 0.3\times 10^{16}$,
|  |  |  |  | $920\pm 40$ | $2.5\pm 0.5\times 10^{13}$
S26 | $5.3\times 10^{4}$ | $45\pm 2$ | $4.9\pm 0.3\times 10^{22}$ | $1.2\pm 0.1$ |  $140\pm 20$, |  $1.0\pm 0.1\times 10^{20}$,
|  |  |  |  | $433\pm 5$ | $4\pm 1\times 10^{17}$

  * a

$T_{\textrm{dust},\;i>1}$ and $N_{i>1}[\textrm{(H+)H}_{2}]$ are respectively
the temperatures and column densities of additional greybody components.

  * b

These are the adopted fitting parameters for N12, in which PACS-70 $\mu$m data
are excluded. When the 70 $\mu$mdata are included, $T_{\textrm{dust}}$,
$N_{\textrm{H}_{2}}$, and $\beta$ become $31\pm 2$ K, $4.3\pm 0.6\times
10^{22}$ cm-2, and $0.6\pm 0.1$, respectively.

Table 4 gives the FIR luminosities and envelope masses derived from greybody
fitting, the bolometric luminosity derived from spline-fitting the data in
Table 5, 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 $\mu$m ($L_{\textrm{bol}}$; see also
Fig. 3), luminosity of the FIR or cold dust component derived from greybody
fitting ($L_{\textrm{FIR}}$), total gas mass derived from $N_{\textrm{H}_{2}}$
within the 3-$\sigma$ source area given the source dimensions in MBS07 at 1.2
mm ($M_{\textrm{c}}$), protostellar luminosity from Transphere modelling
($L_{\textrm{TP}}$), protostellar envelope mass from Transphere modelling
($M_{\textrm{env}}$), the fitted power law index of the Transphere model’s
density profile ($p$), and the deconvolved fitted radius in the 450 $\mu$m
SCUBA-2 data ($r_{\textrm{out}}$). Because of the wavelength restrictions on
Transphere modelling, $L_{\textrm{TP}}$ from greybody fitting is more directly
comparable to $L_{\textrm{FIR}}$ than $L_{\textrm{bol}}$. Note also that
$M_{\textrm{c}}$ is not directly comparable to $M_{\textrm{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_{\textrm{env}}$ were not
systematically smaller.

We used a range of possible $r_{\textrm{out}}$, defined from the minimum and
maximum measured radii of flux contours at 10% of the background-subtracted
maximum flux in the 450 $\mu$m images, before deconvolution with the 450
$\mu$m SCUBA-2 primary beam, to minimise $\chi^{2}$ in $p$ and
$r_{\textrm{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 $\mu$m SCUBA-2 data, while for sources with bright filamentary
protrusions (flagged ’B’), we selected a radial cut free of such structures.
However, $r_{\textrm{out}}$ turned out to be at least somewhat degenerate with
all of the other Transphere fitting parameters. The resulting
$r_{\textrm{out}}$ fits, rounded to the nearest 5000 AU due to precision
limits and then deconvolved by subtracting the $1.64\sigma$ 450 $\mu$m SCUBA-2
beam radius in quadrature222The radius of a Gaussian at 10% of the maximum is
about $1.64\sigma$. For the SCUBA-2 beam at 450 $\mu$m, with a FWHM of 7.9′′,
the $1.64\sigma$ radius is about 5.5′′., 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_{\textrm{out}}$ values were ill-defined—for small sources like N48,
refining $r_{\textrm{out}}$ by 1000 AU produced noticeable improvements, but
for most sources we could only tell that the uncertainties were at least 20%
of $r_{\textrm{out}}$. The uncertainties reported in Table 4 are the greater
of either $\pm 0.2r_{\textrm{out}}$, or the deconvolved differences between
the fitted $r_{\textrm{out}}$ and the upper and lower measured limits used to
constrain the fitting. For $\chi^{2}$ minimisation in envelope mass and
luminosity, it was enough to adopt the nearest multiple of $10^{4}$ AU and
refine the $r_{\textrm{out}}$ by eye using the SED (see §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 $\sim 95\%$ confidence levels on the maps of
$\chi^{2}$.

Table 4: Parameters for PILS-Cygnus sources derived from SED fitting and
Transphere modelling.

Source |  $L_{\textrm{bol}}$ ($L_{\odot}$) |  $L_{\textrm{FIR}}$ ($L_{\odot}$) |  $M_{\textrm{c}}$ ($M_{\odot}$) |  $L_{\textrm{TP}}$ ($L_{\odot}$) |  $M_{\textrm{env}}$ ($M_{\odot}$) | $p$ | flag($p$)(a) |  $r_{\textrm{out}}$ (AU)
---|---|---|---|---|---|---|---|---
N12 | $740^{+80}_{-90}$ |  $330-490$(b) |  $22-54$(b) | $360^{+40}_{-50}$ | $130\pm 30$ | $-1.2^{+0.3}_{-0.1}$ | B | $2.4\pm 0.5\times 10^{4}$
N30 | $2.5^{+0.1}_{-0.3}\times 10^{4}$ | $1.5\pm 0.2\times 10^{4}$ | $220\pm 20$ | $1.4\pm 0.3\times 10^{4}$ | $190\pm 50$ | $-1.0\pm 0.1$ | A | $2.9^{+0.6}_{-0.7}\times 10^{4}$
N48 | $2000\pm 700$ | $1000\pm 300$ | $30\pm 20$ | $900\pm 200$ | $50\pm 10$ | $\gtrsim-0.5$ | B | $1.4^{+0.8}_{-0.3}\times 10^{4}$
N51 | $1300\pm 200$ | $400-700^{\dagger}$ | $130\pm 30$ | $750\pm 70$ | $90\pm 10$ | $-1.1\pm 0.1$ | B | $2.4^{+0.5}_{-0.8}\times 10^{4}$
N53 | $400\pm 200$ | $270\pm 80$ | $120\pm 90$ | $170^{+20}_{-30}$ | $240^{+40}_{-20}$ | $-1.0\pm 0.2$ | B | $1.3\pm 0.5\times 10^{4}$
N63 | $380^{+30}_{-40}$ | $220\pm 70$ | $39\pm 13$ | $310^{+40}_{-80}$ | $80^{+30}_{-10}$ | $-1.2^{+0.3}_{-0.1}$ | A | $2.4^{+0.5}_{-0.9}\times 10^{4}$
S8 | $6400^{+200}_{-1000}$ | $3500\pm 300$ | $63\pm 6$ | $3300\pm 500$ | $90\pm 10$ | $<-1.6$ | B | $3.9^{+0.9}_{-0.8}\times 10^{4}$
S26(c) | $2.10^{+0.06}_{-0.4}\times 10^{5}$ | $5.1\pm 0.7\times 10^{4}$ | $300\pm 30$ | $6.0\pm 0.2\times 10^{4}$ | $530^{+110}_{-150}$ | $-1.3\pm 0.1$ | A | $8\pm 2\times 10^{4}$

  * a

Flag($p$) indicates whether we estimated $p$ using the radially-averaged cut
(denoted A), or the “best” (lowest-$\chi^{2}$) cut (denoted B).

  * b

These ranges express the range of values plus or minus uncertainties that we
get depending on whether 70 $\mu$m data are included in the greybody SED fit.
See also the note about N12 under Table 3.

  * c

S26 is not only more than twice as far away as the other sources, but is also
thought to contain at least one outflow seen pole-on (Trinidad et al., 2003).

N12 in particular was almost as well-fit with PACS-70 $\mu$m data as without,
but the physical conditions derived in each case are dramatically different.
Including PACS-70 $\mu$m data increases $T_{\textrm{dust}}$ to $31\pm 2$ K,
and reduces $\beta$ to $0.6\pm 0.1$ and $N_{\textrm{H}_{2}}$ to $4.3\pm
0.6\times 10^{22}$ cm-2. 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_{\textrm{dust}}$, higher
$N_{\textrm{H}_{2}}$ (higher mass, lower luminosity) results.

Figure 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 $\lambda\geq 70$ $\mu$m are the FIR/submm
component of the SED fits discussed in §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 $\mu$m data point, so the SED fit was only used for $\lambda\geq 160$
$\mu$m while the 70 $\mu$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.

###  3.2 Comparison of Methods

Fig. 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
Appendix B and C . It is unclear from either SED fit whether the 70 $\mu$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 $\mu$m are N48, N51, and N63). It is noteworthy that
for the five lowest-$T_{\textrm{dust}}$ sources, the $\beta$ index of a
(fixed-$T_{\textrm{dust}}$) SED fit to just the 450 and 850 $\mu$m data would
be steeper than the measured $\beta$ 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 $\beta$ 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.

Figure 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 $\lambda<70\,\mu$m were excluded from the fit because the local ISRF is
poorly constrained, and we were generally unable to set it large enough to fit
the short-wavelength data without crashing Transphere.

(a) A. (b) B.

Figure 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.

Half of our sources have mutually consistent $M_{\textrm{c}}$ and
$M_{\textrm{env}}$, but all eight have mutually consistent $L_{\textrm{FIR}}$
and $L_{\textrm{TP}}$. The discrepancies between the envelope masses from each
method can be at least partly accounted for by noting that
$M_{\textrm{env}}\propto r_{\textrm{out}}^{3-p}$: for six of the eight
sources, $r_{\textrm{out}}$ as measured at 450 $\mu$m is substantially larger
than the same at 1.2 mm (denoted $r_{1.2\textrm{mm}}$ in Table 3) despite the
sources supposedly being resolved at both 450 $\mu$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 $\mu$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 $\mu$m, N48 and N53, are in crowded high-background
areas, and may be over-subtracted at 450 $\mu$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 $\beta$ and $p$ should be related, and
Transphere does not fit or incorporate $\beta$ explicitly, but mathematical
effects of $\beta$ 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 $\mu$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-$\chi^{2}$ 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 $\beta$ is an average over a model of the whole source, $\beta$
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 $\beta$ 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+\beta)$ (Chandler et al., 1998; Kenyon et al.,
1993).

We had convolved the SCUBA legacy data333The SCUBA-2 data have a better
resolution and smaller uncertainties, but the secondary beam component
contains a substantially larger fraction of the total flux (40% at 450 $\mu$m,
Holland et al. 2013) and is larger in radius. 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
$\mu$m than 850 $\mu$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 $\beta$ from the SCUBA data alone would be more
reliable than the values of $\beta$ 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 $\mu$m data as
compared to the 350 and 500 $\mu$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_{\textrm{bol}}$ is
preferred given the relative lack of underlying assumptions.

###  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 $\beta$ 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_{\textrm{c}}$ nor
$M_{\textrm{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_{\textrm{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
$\mu$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 $\beta=1.8$ (Ossenkopf & Henning, 1994), which we found to be
too steep. For a given temperature, a slightly steeper $\beta$ 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 $\lambda>70$ $\mu$m assuming $\beta=2$ and optically thin emission,
and excluded data that were too far from the resulting curve (which usually
included the 70 $\mu$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_{\textrm{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 $\mu$m data.
However, we found 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 further two (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 $\beta$, 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 $\lambda>50$ $\mu$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 $\beta$ 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 $\mu$m  _Herschel_ data, but
except in the case of S26, their model SEDs appear to fit the 70 $\mu$m data
and run parallel to but offset from the longer-wavelength data.

####  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, $\chi^{2}$-fitted 450-$\mu$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^{4}\;M_{\odot}$ from van der Tak et al. 2000, 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, 2010b, Kristensen et al.
2012, and van der Tak et al. 2013. In an attempt to further fill out the
sampling shortage between $\sim 50$ and 1000 $L_{\odot}$, 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 $\mu$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.

Figure 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, 2010b (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.

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_{\textrm{env}}$ versus $L_{\textrm{bol}}$,
$r_{\textrm{out}}$ versus $L_{\textrm{bol}}$, and $n$($r=1000$ AU) versus
$L_{\textrm{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_{\textrm{bol}}=1\,L_{\odot}$, and

  * •

ln$\epsilon$, 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_{\textrm{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_{\textrm{bol}}$ instead of at $1\,L_{\odot}$. We somewhat arbitrarily set
the break point at $L_{\textrm{bol}}=60\,L_{\odot}$ 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_{\textrm{env}}$ and $r_{\textrm{out}}$ vs. $L_{\textrm{bol}}$), smaller
than the uncertainties for nearly all available data. Thus, we provide the
following single power law results as our final fits:

$\mathrm{log}\left(\frac{M_{\mathrm{env}}}{M_{\odot}}\right)=0.79^{+0.01}_{-0.02}\mathrm{log}\left(\frac{L_{\mathrm{bol}}}{L_{\odot}}\right)+\mathrm{log}(0.30^{+0.07}_{-0.06})$
(2) $\mathrm{log}\left(\frac{R_{\mathrm{out}}}{\mathrm{AU}}\right)=0.23\pm
0.02\mathrm{log}\left(\frac{L_{\mathrm{bol}}}{L_{\odot}}\right)+\mathrm{log}(5.4^{+0.7}_{-0.6}\times
10^{3})$ (3)
$\mathrm{log}\left(\frac{n_{1000\mathrm{AU}}}{\mathrm{cm}^{-3}}\right)=0.2^{+0.3}_{-0.2}\mathrm{log}\left(\frac{L_{\mathrm{bol}}}{L_{\odot}}\right)+5.9\pm
0.8.$ (4)

For Eqs. 2 and 3, the values of ln$\epsilon$ were $-0.08^{+0.04}_{-0.11}$ and
$-1.0\pm 0.2$, respectively. Because the uncertainties for most values of
$n$($r=1000$ AU) are large, the corresponding nuisance parameter became so
small ($\sim 10^{-17}$) 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_{\textrm{bol}}$, we followed a similar
procedure as for the other three plots to fit a log-linear function in
$L_{\textrm{bol}}$ with a nuisance parameter ln$\epsilon$, 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_{\textrm{bol}}$ and do not vary much over the full range of
$L_{\textrm{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_{\textrm{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 from a 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_{\textrm{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^{2}$
$\mu$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 $\lambda<10$ $\mu$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 $\mu$m than 10 $\mu$m(Lada, 1987). Three of our sources are
more than three orders of magnitude fainter at $\lambda\lesssim 10$ $\mu$m
than at $\lambda\sim 100$ $\mu$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 eight sources 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 at least 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:

  1. 1.

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

  2. 2.

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

  3. 3.

Physical. The envelope may have filaments feeding material into the clump.

  4. 4.

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\lesssim\lambda\lesssim 30$ $\mu$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.
submitted), 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|\lesssim 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
§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
$\gtrsim 40$ pc for most other sources, caveats listed in §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 §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_{\textrm{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 $\sigma_{p}$ of $1.1\pm 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_{\textrm{env}}$ versus
$L_{\textrm{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_{\textrm{bol}}$ between a few tens and a few hundred $L_{\odot}$ for all
relationships except the trend in $r_{\textrm{out}}$ with $L_{\textrm{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_{\textrm{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_{\textrm{env}}$, $r_{\textrm{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_{\textrm{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_{\textrm{env}}$-$L_{\textrm{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_{\textrm{env}}$, even though they are not outliers in the trends of
$r_{\textrm{out}}$ and $n$($r=1000$ AU) with $L_{\textrm{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_{\textrm{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 $\sigma$; often the
definition of $r_{\textrm{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_{\textrm{out}}$ is measured in different 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 in 3).

##  5 Conclusions

In this study, we have deconvolved and fit greybody SEDs to archival IR data
from 1.2 to 850 $\mu$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
$\mu$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. 1.

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

  2. 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. 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. 4.

We establish to first-order that log($M_{\textrm{env}}$)$\;\propto
0.79^{+0.01}_{-0.02}$log($L_{\textrm{bol}}$),
log($r_{\textrm{out}}$)$\;\propto 0.23\pm-0.02$log($L_{\textrm{bol}}$), and
log($n$($r=1000$ AU))$\;\propto 0.2^{+0.3}_{-0.2}$log($L_{\textrm{bol}}$) over
seven decades of $L_{\textrm{bol}}$, with a scatter that is consistently about
an order of magnitude above and below these fits.

  5. 5.

For the aggregate of sources in the literature plus our sources, we determine
that the distribution of density power law indices is consistent with a
constant $1.4\pm 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_{\textrm{env}}$, $r_{\textrm{out}}$, $n$($r=1000$
AU), $p$) clearly warrants further investigation, especially more and larger
surveys that focus on protostars in the $L_{\textrm{bol}}$ range of 50–500
$L_{\odot}$. Nonetheless, the results presented here suggest that the
differences between low- and high-mass sources do not arise from variations in
large-scale physical properties.

## Acknowledgements

The research of RLP and LEK is supported by a research grant (19127) from
VILLUM FONDEN. JKJ 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,444http://www.astropy.org a
community-developed core Python package for Astronomy (Astropy Collaboration
et al., 2013, 2018).

## Data availability

The majority of the continuum data used for this project are publicly
available through IRSA, or in the case of the SCUBA-2 data, through the
Canadian Astronomy Data Centre.

## References

  * Alves et al. (2007) Alves, J., Lombardi, M., & Lada, C. J. 2007, A&A, 462, L17
  * Andre et al. (1993) Andre, P., Ward-Thompson, D., & Barsony, M. 1993, ApJ, 406, 122
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33
  * Barsony (1994) Barsony, M. 1994, in Astronomical Society of the Pacific Conference Series, Vol. 65, Clouds, Cores, and Low Mass Stars, ed. D. P. Clemens & R. Barvainis, 197
  * Beckwith et al. (1990) Beckwith, S. V. W., Sargent, A. I., Chini, R. S., & Guesten, R. 1990, AJ, 99, 924
  * Beerer et al. (2010) Beerer, I. M., Koenig, X. P., Hora, J. L., et al. 2010, ApJ, 720, 679
  * Beltrán et al. (2002) Beltrán, M. T., Estalella, R., Ho, P. T. P., et al. 2002, ApJ, 565, 1069
  * Benedettini et al. (2018) Benedettini, M., Pezzuto, S., Schisano, E., et al. 2018, A&A, 619, A52
  * Berlanas et al. (2019) Berlanas, S. R., Wright, N. J., Herrero, A., Drew, J. E., & Lennon, D. J. 2019, MNRAS, 484, 1838
  * Bernard et al. (2008) Bernard, J.-P., Reach, W. T., Paradis, D., et al. 2008, AJ, 136, 919
  * Blandford & Payne (1982) Blandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883
  * Bonnell et al. (2001) Bonnell, I. A., Bate, M. R., Clarke, C. J., & Pringle, J. E. 2001, MNRAS, 323, 785
  * Bonnell et al. (2004) Bonnell, I. A., Vine, S. G., & Bate, M. R. 2004, MNRAS, 349, 735
  * Bontemps et al. (2010) Bontemps, S., Motte, F., Csengeri, T., & Schneider, N. 2010, A&A, 524, A18
  * Cao et al. (2019) Cao, Y., Qiu, K., Zhang, Q., et al. 2019, ApJS, 241, 1
  * Chandler et al. (1998) Chandler, C. J., Barsony, M., & Moore, T. J. T. 1998, MNRAS, 299, 789
  * Crimier et al. (2010a) Crimier, N., Ceccarelli, C., Alonso-Albi, T., et al. 2010a, A&A, 516, A102
  * Crimier et al. (2009) Crimier, N., Ceccarelli, C., Lefloch, B., & Faure, A. 2009, A&A, 506, 1229
  * Crimier et al. (2010b) Crimier, N., Ceccarelli, C., Maret, S., et al. 2010b, A&A, 519, A65
  * Csengeri et al. (2017) Csengeri, T., Bontemps, S., Wyrowski, F., et al. 2017, A&A, 600, L10
  * Desert et al. (1990) Desert, F. X., Boulanger, F., & Puget, J. L. 1990, A&A, 500, 313
  * Di Francesco et al. (2008) Di Francesco, J., Johnstone, D., Kirk, H., MacKenzie, T., & Ledwosinska, E. 2008, ApJS, 175, 277
  * Downes & Rinehart (1966) Downes, D. & Rinehart, R. 1966, ApJ, 144, 937
  * Duarte-Cabral et al. (2013) Duarte-Cabral, A., Bontemps, S., Motte, F., et al. 2013, A&A, 558, A125
  * Dullemond et al. (2002) Dullemond, C. P., van Zadelhoff, G. J., & Natta, A. 2002, A&A, 389, 464
  * Egan et al. (2003) Egan, M. P., Price, S. D., Kraemer, K. E., et al. 2003, VizieR Online Data Catalog, V/114
  * Ellsworth-Bowers et al. (2015) Ellsworth-Bowers, T. P., Glenn, J., Riley, A., et al. 2015, ApJ, 805, 157
  * Fazio et al. (2004) Fazio, G. G., Hora, J. L., Allen, L. E., et al. 2004, ApJS, 154, 10
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gómez et al. (2021) Gómez, G. C., Vázquez-Semadeni, E., & Palau, A. 2021, MNRAS, 502, 4963
  * Grasdalen et al. (1983) Grasdalen, G. L., Gehrz, R. D., Hackwell, J. A., Castelaz, M., & Gullixson, C. 1983, ApJS, 53, 413
  * Griffin et al. (2010) Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3
  * Habing (1968) Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421
  * Hansen et al. (2004) Hansen, C. J., Kawaler, S. D., & Trimble, V. 2004, Stellar interiors : physical principles, structure, and evolution
  * Hanson (2003) Hanson, M. M. 2003, ApJ, 597, 957
  * Hartmann (1998) Hartmann, L. 1998, Accretion Processes in Star Formation
  * Haschick et al. (1981) Haschick, A. D., Reid, M. J., Burke, B. F., Moran, J. M., & Miller, G. 1981, ApJ, 244, 76
  * Hatchell & van der Tak (2003) Hatchell, J. & van der Tak, F. F. S. 2003, A&A, 409, 589
  * Heger et al. (2003) Heger, A., Fryer, C. L., Woosley, S. E., Langer, N., & Hartmann, D. H. 2003, ApJ, 591, 288
  * Hennebelle & Fromang (2008) Hennebelle, P. & Fromang, S. 2008, A&A, 477, 9
  * Hennemann et al. (2012) Hennemann, M., Motte, F., Schneider, N., et al. 2012, A&A, 543, L3
  * Hildebrand (1983) Hildebrand, R. H. 1983, QJRAS, 24, 267
  * Holland et al. (2013) Holland, W. S., Bintley, D., Chapin, E. L., et al. 2013, MNRAS, 430, 2513
  * Jørgensen et al. (2007a) Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2007a, ApJ, 659, 479
  * Jørgensen et al. (2007b) Jørgensen, J. K., Johnstone, D., Kirk, H., & Myers, P. C. 2007b, ApJ, 656, 293
  * Jørgensen et al. (2002) Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2002, A&A, 389, 908
  * Kauffmann et al. (2008) Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, N. J., I., & Lee, C. W. 2008, A&A, 487, 993
  * Kenyon et al. (1993) Kenyon, S. J., Calvet, N., & Hartmann, L. 1993, ApJ, 414, 676
  * Kong et al. (2021) Kong, S., Arce, H. G., Shirley, Y., & Glasgow, C. 2021, ApJ, 912, 156
  * Könyves et al. (2020) Könyves, V., André, P., Arzoumanian, D., et al. 2020, A&A, 635, A34
  * Kristensen et al. (2012) Kristensen, L. E., van Dishoeck, E. F., Bergin, E. A., et al. 2012, A&A, 542, A8
  * Lada (1987) Lada, C. J. 1987, in IAU Symposium, Vol. 115, Star Forming Regions, ed. M. Peimbert & J. Jugaku, 1
  * Lada & Lada (2003) Lada, C. J. & Lada, E. A. 2003, ARA&A, 41, 57
  * Lommen et al. (2008) Lommen, D., Jørgensen, J. K., van Dishoeck, E. F., & Crapsi, A. 2008, A&A, 481, 141
  * Lynds (1965) Lynds, B. T. 1965, ApJS, 12, 163
  * Markwardt (2009) Markwardt, C. B. 2009, in Astronomical Society of the Pacific Conference Series, Vol. 411, Astronomical Data Analysis Software and Systems XVIII, ed. D. A. Bohlender, D. Durand, & P. Dowler, 251
  * Massey & Thompson (1991) Massey, P. & Thompson, A. B. 1991, AJ, 101, 1408
  * Maud et al. (2015) Maud, L. T., Moore, T. J. T., Lumsden, S. L., et al. 2015, MNRAS, 453, 645
  * McKee & Tan (2003) McKee, C. F. & Tan, J. C. 2003, ApJ, 585, 850
  * McLaughlin & Pudritz (1997) McLaughlin, D. E. & Pudritz, R. E. 1997, ApJ, 476, 750
  * Molinari et al. (2016) Molinari, S., Schisano, E., Elia, D., et al. 2016, A&A, 591, A149
  * Motte et al. (2018) Motte, F., Bontemps, S., & Louvet, F. 2018, ARA&A, 56, 41
  * Motte et al. (2007) Motte, F., Bontemps, S., Schilke, P., et al. 2007, A&A, 476, 1243
  * Mottram et al. (2014) Mottram, J. C., Kristensen, L. E., van Dishoeck, E. F., et al. 2014, A&A, 572, A21
  * Mueller et al. (2002) Mueller, K. E., Shirley, Y. L., Evans, Neal J., I., & Jacobson, H. R. 2002, ApJS, 143, 469
  * Neugebauer et al. (1984) Neugebauer, G., Habing, H. J., van Duinen, R., et al. 1984, ApJ, 278, L1
  * Ortiz-León et al. (2021) Ortiz-León, G. N., Menten, K. M., Brunthaler, A., et al. 2021, arXiv e-prints, arXiv:2105.07471
  * Ossenkopf & Henning (1994) Ossenkopf, V. & Henning, T. 1994, A&A, 291, 943
  * Padoan & Nordlund (2002) Padoan, P. & Nordlund, Å. 2002, ApJ, 576, 870
  * Padoan et al. (2020) Padoan, P., Pan, L., Juvela, M., Haugbølle, T., & Nordlund, Å. 2020, ApJ, 900, 82
  * Palau et al. (2014) Palau, A., Estalella, R., Girart, J. M., et al. 2014, ApJ, 785, 42
  * Pelletier & Pudritz (1992) Pelletier, G. & Pudritz, R. E. 1992, ApJ, 394, 117
  * Pitts & Barnes (2021) Pitts, R. L. & Barnes, P. J. 2021, arXiv e-prints, arXiv:2106.02047
  * Pitts et al. (2019) Pitts, R. L., Barnes, P. J., & Varosi, F. 2019, MNRAS, 484, 305
  * Poglitsch et al. (2010) Poglitsch, A., Waelkens, C., Geis, N., et al. 2010, A&A, 518, L2
  * Reach et al. (2015) Reach, W. T., Heiles, C., & Bernard, J.-P. 2015, ApJ, 811, 118
  * Rieke et al. (2004) Rieke, G. H., Young, E. T., Engelbracht, C. W., et al. 2004, ApJS, 154, 25
  * Russeil et al. (2010) Russeil, D., Zavagno, A., Motte, F., et al. 2010, A&A, 515, A55
  * Rygl et al. (2012) Rygl, K. L. J., Brunthaler, A., Sanna, A., et al. 2012, A&A, 539, A79
  * Sanhueza et al. (2019) Sanhueza, P., Contreras, Y., Wu, B., et al. 2019, ApJ, 886, 102
  * Scoville & Kwan (1976) Scoville, N. Z. & Kwan, J. 1976, ApJ, 206, 718
  * Shirley et al. (2002) Shirley, Y. L., Evans, Neal J., I., & Rawlings, J. M. C. 2002, ApJ, 575, 337
  * Shu (1977) Shu, F. H. 1977, ApJ, 214, 488
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Smartt (2009) Smartt, S. J. 2009, ARA&A, 47, 63
  * Stahler (1994) Stahler, S. W. 1994, ApJ, 422, 616
  * Svoboda et al. (2016) Svoboda, B. E., Shirley, Y. L., Battersby, C., et al. 2016, ApJ, 822, 59
  * Tackenberg et al. (2012) Tackenberg, J., Beuther, H., Henning, T., et al. 2012, A&A, 540, A113
  * Tigé et al. (2017) Tigé, J., Motte, F., Russeil, D., et al. 2017, A&A, 602, A77
  * Torrelles et al. (1997) Torrelles, J. M., Gómez, J. F., Rodríguez, L. F., et al. 1997, ApJ, 489, 744
  * Trinidad et al. (2003) Trinidad, M. A., Curiel, S., Cantó, J., et al. 2003, ApJ, 589, 386
  * Urquhart et al. (2011) Urquhart, J. S., Morgan, L. K., Figura, C. C., et al. 2011, MNRAS, 418, 1689
  * van der Tak et al. (2013) van der Tak, F. F. S., Chavarría, L., Herpin, F., et al. 2013, A&A, 554, A83
  * van der Tak et al. (2000) van der Tak, F. F. S., van Dishoeck, E. F., Evans, Neal J., I., & Blake, G. A. 2000, ApJ, 537, 283
  * Varricatt et al. (2010) Varricatt, W. P., Davis, C. J., Ramsay, S., & Todd, S. P. 2010, MNRAS, 404, 661
  * Vázquez-Semadeni et al. (2009) Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A. K., Ballesteros-Paredes, J., & Klessen, R. S. 2009, ApJ, 707, 1023
  * Vázquez-Semadeni et al. (2019) Vázquez-Semadeni, E., Palau, A., Ballesteros-Paredes, J., Gómez, G. C., & Zamora-Avilés, M. 2019, MNRAS, 490, 3061
  * Williams et al. (2005) Williams, S. J., Fuller, G. A., & Sridharan, T. K. 2005, A&A, 434, 257

##  Appendix A Measured and deconvolved continuum fluxes

Table 5: Measured or deconvolved fluxes of PILS-Cygnus sources.

Filter |  $\lambda$ ($\mu$m) |  $F_{\lambda}$(N12) |  $F_{\lambda}$(N30) |  $F_{\lambda}$(N38)a |  $F_{\lambda}$(N48) |  $F_{\lambda}$(N51)
---|---|---|---|---|---|---
2MASS-J | 1.22 | $\cdots$ | $<1\times 10^{-3}$ | $\cdots$ | $\cdots$ | $<1\times 10^{-3}$
2MASS-H | 1.63 | $<3\times 10^{-3}$ | $8\pm 3\times 10^{-3}$ | $\cdots$ | $\cdots$ | $<3\times 10^{-3}$
2MASS-Ks | 2.19 | $4\pm 2\times 10^{-3}$ | $4\pm 2\times 10^{-3}$ | $\cdots$ | $\cdots$ | $3\pm 2\times 10^{-3}$
IRAC-1 | 3.6 | $1.243\pm.002\times 10^{-2}$ | $6.417\pm 0.004\times 10^{-1}$ | $2.13\pm 0.07\times 10^{-3}$ | $2.9\pm 0.5\times 10^{-4}$ | $1.352\pm 0.1\times 10^{-1}$
IRAC-2 | 4.5 | $3.344\pm 0.002\times 10^{-2}$ | $1.981\pm 0.001$ | $1.34\pm 0.02\times 10^{-2}$ | $4.87\pm 0.02\times 10^{-3}$ | $6.692\pm 0.004\times 10^{-1}$
IRAC-3 | 5.8 | $3.915\pm 0.007\times 10^{-2}$ | $3.902\pm 0.002$ | $1.97\pm 0.03\times 10^{-2}$ | $1.14\pm 0.03\times 10^{-2}$ | $1.617\pm 0.001$
IRAC-4 | 8.0 | $2.911\pm 0.006\times 10^{-2}$ | $5.361\pm 0.003$ | $1.22\pm 0.03\times 10^{-2}$ | $2.42\pm 0.06\times 10^{-2}$ | $2.495\pm 0.003$
MSX-A | 8.28 | $\cdots$ | $8\pm 4$ | $\cdots$ | $\cdots$ | $1.77\pm 0.07$
MSX-C | 12.13 | $\cdots$ | $14\pm 5$ | $\cdots$ | $\cdots$ | $3.6\pm 0.2$
MSX-D | 14.65 | $\cdots$ | $19\pm 6$ | $\cdots$ | $\cdots$ | $6.8\pm 0.4$
MSX-E | 21.4 | $\cdots$ | $54\pm 6$ | $\cdots$ | $\lesssim 28$ | $27\pm 2$
MIPS | 24 | $1.28\pm 0.03$ | $\cdots$ | $\cdots$ | $\cdots$ | $\gtrsim 17$
PACS-B | 70 | $130\pm 20$ | $4800\pm 300$ | $\cdots$ | $280\pm 130$ | $180\pm 50$
PACS-R | 160 | $230\pm 10$ | $4000\pm 200$ | $\cdots$ | $370\pm 110$ | $330\pm 30$
SPIRE-B | 250 | $163\pm 9$ | $2000\pm 100$ | $\cdots$ | $150\pm 80$ | $240\pm 20$
SPIRE-G | 350 | $99\pm 8$ | $660\pm 60$ | $\cdots$ | $\cdots$ | $130\pm 20$
SCUBA-B | 450 | $80\pm 30$ | $360\pm 60$ | $230\pm 110$ | $100\pm 50$ | $140\pm 70$
SPIRE-R | 500 | $53\pm 4$ | $320\pm 20$ | $\cdots$ | $\cdots$ | $60\pm 30$
SCUBA-R | 850 | $9\pm 2$ | $60\pm 20$ | $24\pm 6$ | $15\pm 5$ | $10\pm 5$
Filter |  $\lambda$ ($\mu$m) |  $F_{\lambda}$(N53) |  $F_{\lambda}$(N54)b |  $F_{\lambda}$(N63) |  $F_{\lambda}$(S8)c |  $F_{\lambda}$(S26)c
2MASS-J | 1.22 | $\cdots$ | $2\pm 1\times 10^{-3}$ | $\cdots$ | $2.0\pm 0.1\times 10^{-2}$ | $5\pm 1\times 10^{-3}$
2MASS-H | 1.63 | $\cdots$ | $2\pm 2\times 10^{-3}$ | $\cdots$ | $9.0\pm 0.2\times 10^{-2}$ | $7.1\pm 0.2\times 10^{-2}$
2MASS-Ks | 2.19 | $\cdots$ | $4\pm 2\times 10^{-3}$ | $\cdots$ | $2.58\pm 0.06\times 10^{-1}$ | $1.29\pm 0.03$
IRAC-1 | 3.6 | $8.56\pm 0.05\times 10^{-3}$ | $1.596\pm 0.004\times 10^{-2}$ | $1.77\pm 0.01\times 10^{-3}$ | $8.9\pm 0.2\times 10^{-1}$ | $73.0\pm 0.6$
IRAC-2 | 4.5 | $2.97\pm 0.01\times 10^{-2}$ | $6.913\pm 0.006\times 10^{-2}$ | $9.32\pm 0.02\times 10^{-3}$ | $1.2\pm 0.1$ | $176\pm 4$
IRAC-3 | 5.8 | $2.66\pm 0.03\times 10^{-2}$ | $1.992\pm 0.004\times 10^{-1}$ | $1.033\pm 0.006\times 10^{-2}$ | $3.2\pm 0.2$ | $324\pm 2$
IRAC-4 | 8.0 | $4\pm 2\times 10^{-3}$ | $3.18\pm 0.02\times 10^{-1}$ | $5.76\pm 0.07\times 10^{-3}$ | $5.4\pm 0.4$ | $\gtrsim 190$
MSX-A | 8.28 | $\cdots$ | $2.8\pm 0.1$ | $\cdots$ | $4.4\pm 0.2$ | $310\pm 30$
MSX-C | 12.13 | $\cdots$ | $3.3\pm 0.2$ | $\cdots$ | $14.6\pm 0.7$ | $570\pm 30$
MSX-D | 14.65 | $\cdots$ | $3.4\pm 0.2$ | $\cdots$ | $31\pm 2$ | $800\pm 50$
MSX-E | 21.4 | $\cdots$ | $11\pm 1$ | $\cdots$ | $145\pm 9$ | $1020\pm 60$
MIPS | 24 | $0.52\pm 0.04$ | $\cdots$ | $0.29\pm 0.01$ | $\lesssim 241$ | $\cdots$
IRAS-25 | 25 | $\cdots$ | $\cdots$ | $\cdots$ | $\lesssim 240$ | $\lesssim 1100$
IRAS-60 | 60 | $\cdots$ | $\cdots$ | $\cdots$ | $\lesssim 1050$ | $\lesssim 5300$
PACS-B | 70 | $50\pm 50$ | $\cdots$ | $57\pm 6$ | $1140\pm 60$ | $3600\pm 200$
PACS-R | 160 | $140\pm 30$ | $\cdots$ | $139\pm 8$ | $1040\pm 50$ | $1910\pm 90$
SPIRE-B | 250 | $140\pm 20$ | $\cdots$ | $106\pm 5$ | $390\pm 20$ | $580\pm 30$
SPIRE-G | 350 | $90\pm 20$ | $\cdots$ | $66\pm 4$ | $220\pm 20$ | $270\pm 10$
SCUBA-B | 450 | $110\pm 50$ | $30\pm 40$ | $45\pm 5$ | $160\pm 50$ | $220\pm 110$
SPIRE-R | 500 | $60\pm 10$ | $\cdots$ | $37\pm 3$ | $117\pm 7$ | $107\pm 7$
SCUBA-R | 850 | $12\pm 6$ | $12\pm 6$ | $8\pm 3$ | $30\pm 10$ | $24\pm 5$

  * a

In most filters, N38 could not be identified because the centroid of the
nearest source was almost as far away as N48. The few measurements given
should be treated with extreme caution.

  * b

In most filters, N54 could not be deblended from another object, N52, which
was included in the MBS07 survey but not later works.

  * c

Note that upper limits from IRAS are included for these two sources in
Cygnus-X South.

##  Appendix B SEDs fits

Figures 7-13 show the data (black filled and hollow circles 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 possible with our initial
greybody fitting routine.

Figure 7: Same as Fig. 4, but for N12. It is uncertain whether the PACS-70
$\mu$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 $\mu$m data yields a more realistic
$\beta$ value. Figure 8: Same as Fig. 7, 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 $\lambda>250\;\mu$m. We advise caution in using the
derived SED parameters for this object in any precision analysis. Figure 9:
Same as Fig. 7, but for N51. Figure 10: Same as Fig. 7, but for N53. Figure
11: Same as Fig. 7, but for N63. Figure 12: Same as Fig. 7, but for S8. The
distance to S8 is highly uncertain and it has a spiral arm feature that could
only be accounted for to first order, so SED fitting parameters should be used
with caution. Figure 13: Same as Fig. 7, but for S26. We note the change in
y-axis scaling.

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

Figures 14-20 are the equivalent of Fig. 5 for all other sources besides N30.

(a) A. (b) B.

Figure 14: Same as Fig. 5, but for N12 (using the fit that excludes the 70
$\mu$m data).

(a) A. (b) B.

Figure 15: Same as Fig. 5, but for N48.

(a) A. (b) B.

Figure 16: Same as Fig. 5, but for N51.

(a) A. (b) B.

Figure 17: Same as Fig. 5, but for N53.

(a) A. (b) B.

Figure 18: Same as Fig. 5, but for N63.

(a) A. (b) B.

Figure 19: Same as Fig. 5, but for S8.

(a) A. (b) B.

Figure 20: 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_{\textrm{env}}$,
$r_{\textrm{out}}$, $|p|$, and $n$($r=1000$ AU) with $L_{\textrm{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
parameters of each of the pair of fits to each relationship, except for the
fit to $|p|$ as a function of $L_{\textrm{bol}}$. As we discussed in the main
text, we did not find a significant correlation between $p$ and
$L_{\textrm{bol}}$.

(a) Single Power Law (b) Broken Power Law

Figure 21: (a) Posterior distributions of the power-law exponent ($\xi$),
normalisation ($M_{o}$), and residual variance parameter (ln$\epsilon_{M}$)
for the single power-law fit to $M_{\textrm{env}}$ vs. $L_{\textrm{bol}}$,
shown as the solid black line in the top left panel of Fig. 6. (b) Posterior
distributions of the power-law exponents ($\xi_{1}$ and $\xi_{2}$),
normalisation ($M_{bk}$), and residual variance parameter (ln$\epsilon_{M}$)
for the broken power-law fit to $M_{\textrm{env}}$ vs. $L_{\textrm{bol}}$,
shown as the dashed black line in the top left panel of Fig. 6.

(a) Single Power Law (b) Broken Power Law

Figure 22: (a) Posterior distributions of the power-law exponent ($\alpha$),
normalisation ($R_{o}$), and residual variance parameter (ln$\epsilon_{R}$)
for the single power-law fit to $r_{\textrm{out}}$ vs. $L_{\textrm{bol}}$,
shown as the solid black line in the top right panel of Fig. 6. (b) Posterior
distributions of the power-law exponents ($\alpha_{1}$ and $\alpha_{2}$),
normalisation ($R_{bk}$), and residual variance parameter (ln$\epsilon_{R}$)
for the broken power-law fit to $r_{\textrm{out}}$ vs. $L_{\textrm{bol}}$,
shown as the dashed black line in the top right panel of Fig. 6.

(a) Single Power Law (b) Broken Power Law

Figure 23: (a) Posterior distributions of the power-law exponent ($\eta$) and
normalisation ($n_{o}$) for the single power-law fit to $n$($r=1000$ AU) vs.
$L_{\textrm{bol}}$, shown as the solid black line the bottom right panel of in
Fig. 6. (b) Posterior distributions of the power-law exponents ($\eta_{1}$ and
$\eta_{2}$), normalisation ($n_{bk}$), and residual variance parameter
(ln$\epsilon_{n}$) for the broken power-law fit to $n$($r=1000$ AU) vs.
$L_{\textrm{bol}}$, shown as the dashed black line in the bottom right panel
of Fig. 6.