The SOFIA Massive (SOMA) Star Formation Survey. II. High Luminosity Protostars

TL;DR

This study uses multi-wavelength observations and radiative transfer modeling to analyze high-luminosity massive protostars, revealing their properties, environments, and potential formation conditions within the Turbulent Core Accretion framework.

Contribution

It provides detailed SED fitting and property estimates for high-mass protostars, expanding understanding of their formation and environment, especially in clustered regions.

Findings

Protostars have masses 12-64 M_sun and luminosities 10^4-10^6 L_sun.

Models fit most SEDs well, but clustering affects long-wavelength fits.

No strong evidence that high surface density is necessary for massive star formation.

Abstract

We present multi-wavelength images observed with SOFIA-FORCAST from $\sim$10 to 40 $\mu$m of seven high luminosity massive protostars, as part of the SOFIA Massive (SOMA) Star Formation Survey. Source morphologies at these wavelengths appear to be influenced by outflow cavities and extinction from dense gas surrounding the protostars. Using these images, we build spectral energy distributions (SEDs) of the protostars, also including archival data from Spitzer, Herschel and other facilities. Radiative transfer (RT) models of Zhang & Tan (2018), based on Turbulent Core Accretion theory, are then fit to the SEDs to estimate key properties of the protostars. Considering the best five models fit to each source, the protostars have masses $m_{*} \sim 12-64 \: M_{\odot}$ accreting at rates of $\dot{m}_{*} \sim 10^{-4}-10^{-3} \: M_{\odot} \: \rm yr^{-1}$ inside cores of initial masses $M_{c} \sim 100-500 \: M_{\odot}$ embedded in clumps with mass surface densities $\Sigma_{\rm cl} \sim 0.1-3 \: \rm g \: cm^{-2}$ and span a luminosity range of $10^{4} -10^{6} \: L_{\odot}$. Compared with the first eight protostars in Paper I, the sources analyzed here are more luminous, and thus likely to be more massive protostars. They are often in a clustered environment or have a companion protostar relatively nearby. From the range of parameter space of the models, we do not see any evidence that $\Sigma_{\rm cl}$ needs to be high to form these massive stars. For most sources the RT models provide reasonable fits to the SEDs, though the cold clump material often influences the long wavelength fitting. However, for sources in very clustered environments, the model SEDs may not be such a good description of the data, indicating potential limitations of the models for these regions.

Tables

Table 1. Table 1: SOFIA FORCAST Observations: Observation Dates & Exposure Times (seconds)

Source	R.A.(J2000)	Decl.(J2000)	$d$ (kpc)	Obs. Date	7.7$\mu m$	19.7$\mu m$	31.5$\mu m$	37.1$\mu m$
G45.12+0.13	19h13m27$\stackrel{s}{.}$859	$+$10$°$53$\prime$36$\stackrel{\prime \prime }{.}$645	7.4	2016 Sep 17	2443	882	623	1387
G309.92+0.48	13h50m41$\stackrel{s}{.}$847	$-$61$°$35$\prime$10$\stackrel{\prime \prime }{.}$40	5.5	2016 Jul 14	291	828	532	1691
G35.58-0.03	18h56m22$\stackrel{s}{.}$563	$+$02$°$20$\prime$27$\stackrel{\prime \prime }{.}$660	10.2	2016 Sep 20	335	878	557	1484
IRAS 16562-3959	16h59m41$\stackrel{s}{.}$63	$-$40$°$03$\prime$43$\stackrel{\prime \prime }{.}$61	1.7	2016 Jul 17	1461	772	502	1243
G305.20+0.21	13h11m10$\stackrel{s}{.}$49	$-$62$°$34$\prime$38$\stackrel{\prime \prime }{.}$8	4.1	2016 Jul 18	1671	763	539	1028
G49.27-0.34	19h23m06$\stackrel{s}{.}$61	$+$14$°$20$\prime$12$\stackrel{\prime \prime }{.}$0	5.55	2016 Sep 20	290	716	664	1307
G339.88-1.26	16h52m04$\stackrel{s}{.}$67	$-$46$°$08$\prime$34$\stackrel{\prime \prime }{.}$16	2.1	2016 Jul 20	1668	830	527	1383

Table 2. Table 2: Integrated Flux Densities

Facility	$\lambda$	Fλ,fix aaFlux density derived with a fixed aperture size of the 70 $\mu$m data.	Fλ,var bbFlux density derived with various aperture sizes.	Rap ccAperture radius.	Fλ,fix	Fλ,var	Rap	Fλ,fix	Fλ,var	Rap	Fλ,fix	Fλ,var	Rap	Fλ,fix	Fλ,var	Rap	Fλ,fix	Fλ,var	Rap	Fλ,fix	Fλ,var	Rap
	($\mu$m)	(Jy)	(Jy)	($″$)	(Jy)	(Jy)	($″$)	(Jy)	(Jy)	($″$)	(Jy)	(Jy)	($″$)	(Jy)	(Jy)	($″$)	(Jy)	(Jy)	($″$)	(Jy)	(Jy)	($″$)
		G45.12+0.13			G309.92+0.48			G35.58-0.03			IRAS 16562			G305.20+0.21			G49.27-0.34			G339.88-1.26
WISE	3.4	…	…	…	…	…	…	…	…	…	1.62	0.62	12.0	…	…	…	…	…	…	…	…	…
											(2.53)	(0.89)
Spitzer/IRAC	3.6	5.04	2.83	12.0	2.68	1.56	6.0	0.51	0.05	6.0	…	…	…	1.15	0.62	4.6	0.10	0.05	7.7	1.10	0.39	12.0
		(5.64)	(3.04)		(2.93)	(1.65)		(0.76)	(0.09)					(1.44)	(0.69)		(0.32)	(0.07)		(1.65)	(0.56)
Gemini/T-ReCS	3.8	…	…	…	…	…	…	…	…	…	…	…	…	0.82	0.83	2.0	…	…	…	…	…	…
														(1.10)	(0.87)
Spitzer/IRAC	4.5	7.30	4.79	12.0	4.75	2.92	6.0	0.52	0.12	6.0	…	…	…	2.89	2.00	4.6	0.77	0.55	7.7	2.90	1.75	12.0
		(7.87)	(5.08)		(4.98)	(3.12)		(0.73)	(0.17)					(3.18)	(2.11)		(0.99)	(0.59)		(3.67)	(2.01)
WISE	4.6	…	…	…	…	…	…	…	…	…	4.20	0.73	12.0	…	…	…	…	…	…	…	…	…
											(5.42)	(1.56)
Gemini/T-ReCS	4.7	…	…	…	…	…	…	…	…	…	…	…	…	2.95	2.77	2.0	…	…	…	…	…	…
														(3.35)	(2.85)
Spitzer/IRAC	5.8	41.05	22.97	12.0	18.39	15.11	14.5	2.22	0.30	6.0	…	…	…	6.95	4.14	4.6	2.07	1.55	7.7	…	…	…
		(45.17)	(24.95)		(19.66)	(16.18)		(3.91)	(0.60)					(9.24)	(4.56)		(3.75)	(1.73)
SOFIA/FORCAST	7.7	103.95	57.19	12.0	48.04	35.46	14.5	8.95	0.92	6.0	77.53	58.03	12.0	24.90	13.70	4.6	4.10	2.77	7.7	11.11	1.09	6.0
		(92.89)	(60.84)		(47.17)	(37.55)		(5.83)	(1.62)		(81.72)	(60.63)		(26.31)	(14.95)		(4.20)	(3.02)		(17.69)	(1.67)
Gemini/T-ReCS	7.9	…	…	…	…	…	…	…	…	…	…	…	…	9.39	12.05	2.0	…	…	…	…	…	…
														(12.88)	(12.28)
Spitzer/IRAC	8.0	72.80	25.68	12.0	28.21	20.29	14.5	5.15	0.73	6.0	…	…	…	17.92	8.62	4.6	2.34	1.78	7.7	…	…	…
		(84.64)	(30.73)		(33.57)	(22.56)		(9.58)	(1.46)					(23.82)	(9.75)		(7.00)	(2.18)
Gemini/T-ReCS	8.8	…	…	…	…	…	…	…	…	…	…	…	…	12.81	14.85	2.0	…	…	…	…	…	…
														(15.58)	(15.01)
Gemini/T-ReCS	9.7	…	…	…	…	…	…	…	…	…	…	…	…	17.14	17.24	2.0	…	…	…	…	…	…
														(17.65)	(17.41)
Gemini/T-ReCS	10.4	…	…	…	…	…	…	…	…	…	…	…	…	25.66	25.19	2.0	…	…	…	…	…	…
														(26.33)	(25.41)
Gemini/T-ReCS	11.7	…	…	…	79.11	78.60	6.0	2.19	2.15	6.0	…	…	…	40.91	44.94	2.0	…	…	…	…	…	…
					(80.37)	(79.26)		(2.25)	(2.21)					(47.11)	(45.33)
Gemini/T-ReCS	12.3	…	…	…	…	…	…	…	…	…	…	…	…	54.23	56.18	2.0	…	…	…	…	…	…
														(59.14)	(56.68)
Gemini/T-ReCS	18.3	…	…	…	…	…	…	…	…	…	…	…	…	137	161	2.0	…	…	…	…	…	…
														(174)	(164)
SOFIA/FORCAST	19.7	1128	976	12.0	380	345	10.0	21.78	18.40	7.0	254	212	12.0	282	194	4.6	2.97	2.12	11.0	26.94	24.87	6.0
		(1087)	(988)		(376)	(350)		(25.61)	(19.32)		(241)	(214)		(280)	(201)		(2.84)	(2.25)		(19.63)	(25.66)
Gemini/T-ReCS	24.5	…	…	…	…	…	…	…	…	…	…	…	…	311	375	2.0	…	…	…	…	…	…
														(428)	(385)
SOFIA/FORCAST	31.5	3077	2345	12.0	1896	1700	12.0	276	210	7.7	2078	1758	16.0	687	521	7.7	63.37	41.77	11.0	720	541	7.7
		(3048)	(2423)		(1899)	(1735)		(275)	(221)		(2073)	(1797)		(710)	(546)		(69.24)	(45.96)		(714)	(566)
SOFIA/FORCAST	37.1	4126	2952	12.0	2601	2298	12.0	525	365	7.7	3015	2444	16.0	892	654	7.7	89.39	58.83	11.0	1202	815	7.7
		(4112)	(3082)		(2607)	(2352)		(531)	(394)		(3032)	(2531)		(925)	(694)		(86.18)	(62.34)		(1210)	(870)
Herschel/PACS	70.0	5848	5848	48.0	3403	3403	32.0	1538	1538	25.6	…	…	…	1250	1250	16.0	449	449	28.8	3610	3610	32.0
		(6205)	(6205)		(3536)	(3536)		(1647)	(1647)					(1617)	(1617)		(593)	(593)		(3846)	(3846)
Herschel/PACS	160.0	3517	3517	48.0	2088	2088	32.0	968	968	25.6	…	…	…	644	644	16.0	864	864	28.8	2723	2723	32.0
		(4045)	(4045)		(2454)	(2454)		(1193)	(1193)					(1032)	(1032)		(1198)	(1198)		(3046)	(3046)
Herschel/SPIRE	250.0	1506	1506	48.0	…	…	…	404	404	25.6	…	…	…	234	234	16.0	517	517	28.8	…	…	…
		(1796)	(1796)					(545)	(545)					(433)	(433)		(736)	(736)
Herschel/SPIRE	350.0	469	469	48.0	289	289	32.0	129	129	25.6	…	…	…	60	60	16.0	193	193	28.8	…	…	…
		(591)	(591)		(395)	(395)		(191)	(191)					(143)	(143)		(292)	(292)
Herschel/SPIRE	500.0	136	136	48.0	80	80	32.0	30.87	30.87	25.6	267	267	32.0	20.61	20.61	16.0	54.19	54.19	28.8	…	…	…
		(187)	(187)		(122)	(122)		(56.16)	(56.16)		(374)	(374)		(59.45)	(59.45)		(94.92)	(94.92)

Table 3. Table 3: Summary of Properties of the Protostellar Sources

Source	Radio emission?	Outflow?	Multiple (proto)stars within 20″?	What regulates the MIR morphology?
G45.12+0.13	UC HII	Two	ClusteraaBased on radio sources from Vig et al. 2006..	Ionized medium
G309.92+0.48	HC HII	…	MIR companion.	Outflow cavities
			Resolved.	or flared disk surface?
G35.58-0.03	UC HII	N	Nearby HII region with an outflow.	Outflows from nearby sources?
			Low-mass YSObbBased on multi-wavelength MIR and NIR data.?
IRAS 16562-3959	HC HII with jet	Two	ClusterccBased on X-ray sources from Montes et al. 2018..	Outflow cavities
G305.20+0.21	N	…	Nearby HII region. MIR companion.	Outflow cavities
			ResolvedddWe suspect here the resolved structures are more likely to be emission separated by optically thick disk rather than two distinct protostars.. Low-mass YSObbBased on multi-wavelength MIR and NIR data.?	or flared disk surface?
G49.27-0.34	Y	…	Radio companion. MIR companion.	…
G339.88-1.26	Jet	Two	MIR companion. ResolvedeeWe suspect here the resolved structures are emission tracing the outflow cavities rather than multiple distinct protostars. See also Zhang et al. 2019.. BinaryffBased on the fact of two outflows.?	Outflow cavities
			Low-mass YSObbBased on multi-wavelength MIR and NIR data.? Nearby HII region?	and extinction

Table 4. Table 4: Parameters of the Best Five Fitted Models

Source	${\chi }^2$/N	$M_c$	${\Sigma }_{\mathrm{cl}}$	$R_{\mathrm{core}}$	$m_{\ast }$	${\theta }_{\mathrm{view}}$	$A_V$	$M_{\mathrm{env}}$	${\theta }_{w,\mathrm{esc}}$	${\dot{M}}_{\mathrm{disk}}$	$L_{\mathrm{bol},\mathrm{iso}}$	$L_{\mathrm{bol}}$
		($M_{\odot }$)	(g ${\mathrm{cm}}^{-2}$)	(pc) ($″$)	($M_{\odot }$)	(°)	(mag)	($M_{\odot }$)	(deg)	($M_{\odot }$/yr)	($L_{\odot }$)	($L_{\odot }$)
G45.12+0.13	54.39	480	1.0	0.161 ( 4 )	64.0	34	0.0	325	32	1.2$\times {10}^{-3}$	6.5$\times {10}^5$	8.4$\times {10}^5$
$d$ = 7.4 kpc	63.23	480	1.0	0.161 ( 4 )	48.0	29	0.0	367	25	1.1$\times {10}^{-3}$	4.5$\times {10}^5$	5.4$\times {10}^5$
$R_{ap}$ = 48 ″	65.40	480	3.2	0.091 ( 3 )	24.0	13	0.0	441	12	2.0$\times {10}^{-3}$	1.1$\times {10}^6$	2.9$\times {10}^5$
= 1.72 pc	66.41	400	3.2	0.083 ( 2 )	24.0	13	0.0	362	13	1.9$\times {10}^{-3}$	1.3$\times {10}^6$	3.0$\times {10}^5$
	69.30	240	3.2	0.064 ( 2 )	32.0	29	0.0	175	23	1.9$\times {10}^{-3}$	4.5$\times {10}^5$	5.0$\times {10}^5$
G309.92+0.48	2.82	320	3.2	0.074 ( 3 )	24.0	22	12.1	277	15	1.8$\times {10}^{-3}$	3.3$\times {10}^5$	3.1$\times {10}^5$
$d$ = 5.5 kpc	3.90	480	1.0	0.161 ( 6 )	48.0	29	39.4	367	25	1.1$\times {10}^{-3}$	4.5$\times {10}^5$	5.4$\times {10}^5$
$R_{ap}$ = 32 ″	4.38	240	3.2	0.064 ( 2 )	32.0	34	17.2	175	23	1.9$\times {10}^{-3}$	3.2$\times {10}^5$	5.0$\times {10}^5$
= 0.85 pc	4.71	240	3.2	0.064 ( 2 )	24.0	29	0.0	194	18	1.6$\times {10}^{-3}$	2.6$\times {10}^5$	3.1$\times {10}^5$
	4.97	400	1.0	0.147 ( 6 )	48.0	34	4.0	289	29	1.0$\times {10}^{-3}$	3.0$\times {10}^5$	5.3$\times {10}^5$
G35.58-0.03	1.70	480	3.2	0.091 ( 2 )	24.0	22	16.2	441	12	2.0$\times {10}^{-3}$	2.9$\times {10}^5$	2.9$\times {10}^5$
$d$ = 10.2 kpc	2.14	400	3.2	0.083 ( 2 )	24.0	22	46.5	362	13	1.9$\times {10}^{-3}$	3.0$\times {10}^5$	3.0$\times {10}^5$
$R_{ap}$ = 26 ″	3.41	320	3.2	0.074 ( 1 )	24.0	29	35.4	277	15	1.8$\times {10}^{-3}$	2.7$\times {10}^5$	3.1$\times {10}^5$
= 1.27 pc	4.28	480	1.0	0.161 ( 3 )	48.0	34	39.4	367	25	1.1$\times {10}^{-3}$	3.0$\times {10}^5$	5.4$\times {10}^5$
	4.49	480	1.0	0.161 ( 3 )	64.0	39	72.7	325	32	1.2$\times {10}^{-3}$	3.6$\times {10}^5$	8.4$\times {10}^5$
IRAS16562	0.53	400	0.1	0.465 ( 56 )	32.0	44	100.0	304	29	1.5$\times {10}^{-4}$	9.2$\times {10}^4$	1.6$\times {10}^5$
$d$ = 1.7 kpc	0.64	480	0.1	0.510 ( 62 )	24.0	71	55.6	418	21	1.4$\times {10}^{-4}$	5.7$\times {10}^4$	8.7$\times {10}^4$
$R_{ap}$ = 32 ″	0.65	480	0.1	0.510 ( 62 )	32.0	48	100.0	391	26	1.6$\times {10}^{-4}$	9.8$\times {10}^4$	1.6$\times {10}^5$
= 0.26 pc	0.67	320	0.3	0.234 ( 28 )	16.0	22	17.2	283	16	2.5$\times {10}^{-4}$	5.3$\times {10}^4$	6.1$\times {10}^4$
	0.83	120	3.2	0.045 ( 6 )	16.0	29	100.0	90	21	1.1$\times {10}^{-3}$	1.0$\times {10}^5$	1.2$\times {10}^5$
G305.20+0.21	0.79	80	3.2	0.037 ( 2 )	24.0	48	14.1	35	37	1.1$\times {10}^{-3}$	7.5$\times {10}^4$	2.6$\times {10}^5$
$d$ = 4.1 kpc	0.92	100	3.2	0.041 ( 2 )	32.0	51	18.2	37	42	1.2$\times {10}^{-3}$	7.9$\times {10}^4$	3.5$\times {10}^5$
$R_{ap}$ = 16 ″	0.97	160	1.0	0.093 ( 5 )	32.0	44	13.1	88	39	5.9$\times {10}^{-4}$	8.2$\times {10}^4$	2.3$\times {10}^5$
= 0.32 pc	1.04	80	3.2	0.037 ( 2 )	16.0	34	8.1	50	27	9.5$\times {10}^{-4}$	7.2$\times {10}^4$	1.1$\times {10}^5$
	1.11	160	3.2	0.052 ( 3 )	48.0	58	16.2	59	45	1.6$\times {10}^{-3}$	9.0$\times {10}^4$	6.4$\times {10}^5$
G49.27-0.34	1.87	240	3.2	0.064 ( 2 )	12.0	22	54.5	219	12	1.2$\times {10}^{-3}$	4.5$\times {10}^4$	4.8$\times {10}^4$
$d$ = 5.55 kpc	1.96	200	3.2	0.059 ( 2 )	12.0	22	92.9	179	13	1.1$\times {10}^{-3}$	4.9$\times {10}^4$	5.2$\times {10}^4$
$R_{ap}$ = 29 ″	2.18	320	3.2	0.074 ( 3 )	12.0	22	0.0	302	10	1.3$\times {10}^{-3}$	4.7$\times {10}^4$	4.9$\times {10}^4$
= 0.77 pc	2.37	160	3.2	0.052 ( 2 )	12.0	29	77.8	139	15	1.0$\times {10}^{-3}$	4.4$\times {10}^4$	5.3$\times {10}^4$
	2.73	120	3.2	0.045 ( 2 )	12.0	34	73.7	99	18	9.6$\times {10}^{-4}$	3.6$\times {10}^4$	5.2$\times {10}^4$
G339.88-1.26	2.21	400	0.3	0.262 ( 26 )	12.0	22	17.2	373	11	2.3$\times {10}^{-4}$	3.7$\times {10}^4$	4.0$\times {10}^4$
$d$ = 2.1 kpc	2.30	320	0.3	0.234 ( 23 )	12.0	68	6.1	293	13	2.2$\times {10}^{-4}$	3.3$\times {10}^4$	4.0$\times {10}^4$
$R_{ap}$ = 32 ″	2.48	480	0.3	0.287 ( 28 )	12.0	22	7.1	459	10	2.5$\times {10}^{-4}$	3.8$\times {10}^4$	4.0$\times {10}^4$
= 0.33 pc	2.62	320	0.3	0.234 ( 23 )	16.0	22	90.9	283	16	2.5$\times {10}^{-4}$	5.3$\times {10}^4$	6.1$\times {10}^4$
	2.84	120	3.2	0.045 ( 4 )	12.0	44	0.0	99	18	9.6$\times {10}^{-4}$	3.3$\times {10}^4$	5.2$\times {10}^4$

Full text

ISM: jets and outflows — dust — stars: formation — stars: winds, outflows — stars: early-type — infrared radiation — ISM: individual (G45.12+0.13, G309.92+0.48, G35.58-0.03, IRAS 16562-3959, G305.20+0.21, G49.27-0.34, G339.88-1.26)

††slugcomment: DRAFT

1 Introduction

Massive stars play a key role in the regulation of galaxy environments and their overall evolution, yet there is no consensus on their formation mechanism. Theories range from Core Accretion (e.g., McLaughlin & Pudritz 1996; McKee & Tan 2003 [MT03]) in which massive stars form via a monolithic collapse of a massive core, to Competitive Accretion (e.g., Bonnell et al. 2001; Wang et al. 2010) in which massive stars have most of the mass reservoir joining later and form hand in hand with the formation of a cluster of mostly low-mass stars, to Protostellar Collisions (Bonnell et al. 1998). The confusion remains partly due to the difficulty of observations towards massive star formation given the typically large distances and high extinction of the regions.

Outflows appear to be a ubiquitous phenomenon in the formation of stars of all masses. They may limit the formation efficiency from a core since they expel material along polar directions. The resulting outflow cavities have been proposed to affect the appearance of massive sources in the mid-IR (MIR) up to $\sim$40 $\mu$m (De Buizer 2006; Zhang et al. 2013a) and this is seen in radiative transfer (RT) calculations of the Turbulent Core Model of MT03 (e.g., Zhang et al. 2013b, Zhang et al. 2014; Zhang & Tan 2018).

Motivated by the need of observations of a larger sample of massive protostars to test theoretical models of massive star formation, we are carrying out the SOFIA Massive (SOMA) Star Formation Survey (PI: Tan). The overall goal is to obtain $\sim$ 10 to 40 $\mu$m images with the SOFIA-FORCAST instrument of a sample of $\gtrsim$ 50 high- and intermediate-mass protostars over a range of evolutionary stages and environments, and then compare the observed spectral energy distributions (SEDs) and image intensity profiles with theoretical models. The results and SED analysis of the first 8 sources of the survey have been published by De Buizer et al. (2017) (hereafter Paper I).

In this paper, we now present the next seven most luminous protostars from the sample of completed observations, which are expected to be the highest-mass protostars. In this work we still focus on the SED analysis. Comparison with the image intensity profiles will be presented in a future paper. The observations and data used are described in §2. The analysis methods are described in §3. We present the MIR imaging and SED fitting results in §4 and discuss these results and their implications in §5. A summary is given in §6.

2 Observations

2.1 SOFIA Data

The following seven sources, listed in order of decreasing isotropic bolometric luminosity, were observed by SOFIA111SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institute (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. (Young et al. 2012) with the FORCAST instrument (Herter et al. 2013) (see Table 1): G45.12+0.13; G309.92+0.48; G35.58-0.03; IRAS 16562-3959; G305.20+0.21; G49.27-0.34; G339.88-1.26.

SOFIA data were calibrated by the SOFIA pipeline with a system of stellar calibrators taken across all flights in a flight series and applied to all targets within that flight series (see also the FORCAST calibration paper by Herter et al. 2013). Corrections were also made for the airmass of the sources. The main uncertainty in the SOFIA calibrations is caused by the apparent variability in the flux of the standard stars throughout the flight and from flight to flight due to changing atmospheric conditions. The calibration error is estimated to be in the range $\sim$ 3% - 7%.

2.2 Other IR Data

For all objects, data were retrieved from the Spitzer Heritage Archive from all four IRAC (Fazio et al. 2004) channels (3.6, 4.5, 5.8 and 8.0 $\mu$m). In some cases, the sources are so bright that they are saturated in the IRAC images and so these could not be used to derive accurate fluxes. For IRAS 16562, we used unsaturated WISE archival data (3.4 $\mu$m and 4.6 $\mu$m) as a substitute.

We also incorporated publicly-available imaging observations performed with the Herschel Space Observatory222Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. The Herschel data used in this paper are taken from the Level 2 (flux-calibrated) images provided by the Herschel Science Center via the NASA/IPAC Infrared Science Archive (IRSA), which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with NASA. (Pilbratt et al. 2010) and its PACS (Poglitsch et al. 2010) and SPIRE (Griffin et al. 2010) instruments at 70, 160, 250, 350 and 500 $\mu$m.

In addition to using these data for deriving multi-wavelength flux densities of our sources, the Spitzer 8 $\mu$m and Herschel 70 $\mu$m images are presented for comparison with our SOFIA images in §4.1. We note that the data being analyzed here were typically collected within a time frame of about 10 years (i.e., for the Spitzer, Herschel, and SOFIA observations).

We also present previously unpublished Gemini 8-m data taken with the instrument T-ReCS (De Buizer & Fisher 2004) for sources G309.92, G35.58, and G305.20. For both G309.92 and G35.58, only 11.7 $\mu$m data were taken, with on-source exposures times of 304s and 360s, respectively. For G305.20, we have images through ten T-ReCS filters from 3.8 $\mu$m (L-band) to 24.5 $\mu$m, all with an exposure time of 130s. Most T-ReCS filters have modest flux calibration errors (for MIR observations) with standard deviations between 2 and 10%. For instance, the 11.7 $\mu$m filter has a 1-sigma flux calibration error of 3%. Flux calibration through certain filters, however, is more difficult due to the presence of various atmospheric absorption lines contaminating the filter bandpass, some of which can be highly variable. Those filters most affected are the 7.7 $\mu$m (21%), 12.3 m (19%), 18.3 m (15%), and 24.6 m (23%) filters (De Buizer et al. 2005).

NIR images from the VISTA/VVV3 (Minniti et al. 2010) and the WFCAM/UKIDSS (Lawrence et al. 2007) surveys are also used to investigate the environments of the protostellar sources and look for association with the MIR counterparts.

2.3 Astrometry

The absolute astrometry of the SOFIA data comes from matching the centroids of point sources in the SOFIA 7 $\mu$m image with the Spitzer 8 $\mu$m image (or shorter IRAC wavelength, if saturated at 8 $\mu$m). The relative astrometry between the four SOFIA images is reduced to be better than 0.4″, which is around half a FORCAST pixel. Thus the astrometry precision is about 0.1″for the SOFIA 7 $\mu$m image and 0.4″for longer wavelength SOFIA images. The Herschel data can also be off in their absolute astrometry by up to 5$\arcsec$. For all targets in this survey, we were able to find point sources in common between the Herschel image and sources found in the SOFIA or Spitzer field of view that allowed us to correct the Herschel absolute astrometry. The astrometry is then assumed to have errors of less than 1$\arcsec$.

The Gemini images are calibrated using the Spitzer data and the astrometry precision is better than $\sim$ 0.2″. The archival WISE data and NIR data from the VVV survey and the UKIDSS survey were calibrated using 2MASS point source catalog and should have a positional accuracy $<$ 0.1″.

3 Methods

3.1 SED Construction

We follow the methods in Paper I and use PHOTUTILS, a Python package, to measure the flux photometry. When building the SEDs, we try two different methods. One is using fixed aperture size for all wavelengths, which is our fiducial case. The aperture size is mainly based on the Herschel 70 $\mu$m image, which is typically close to the peak of the SED, in order to capture the most flux from the source, while minimizing contamination from other sources. We assume this is the “core” scale from which the protostar forms as described in the Turbulent Core Model (MT03). If there is no Herschel data available, we use the SOFIA 37 $\mu$m image to determine the aperture size. Sometimes we see multiple IR peaks in the aperture at shorter wavelengths, but without corresponding resolved structures at longer wavelengths, as in G45.12, G309.92, G35.58, and G49.27. This is a combined effect of larger beam sizes at the longer wavelengths and the fact that the emission from the secondary sources appears to be weaker at longer wavelengths. Note that due to the limited size of the field of view of the Gemini images, even for the fixed aperture method, we adopt an aperture radius of 9″, 9″, and 10″ for the photometry of the Gemini images of G309.92, G35.58 and G305.20, respectively, which are the largest aperture sizes possible to allow for background subtraction in each image.

The alternate method is to use variable aperture sizes for each wavelength $<70\>\mu$m. In this case, we typically use smaller apertures at shorter wavelength to exclude secondary sources that appear resolved from the main massive protostar in the fiducial aperture in the Spitzer and SOFIA images and compare the effects on the SEDs. The aperture is always centered at the radio continuum source (or the location of the methanol maser if there is no radio emission as in G305.20), where we assume the protostar is located.

After measuring the flux inside the aperture, we carry out background subtraction using the median flux density in an annular region extending from 1 to 2 aperture radii, as in Paper I, to remove general background and foreground contamination and the effect of a cooler, more massive clump surrounding the core at long wavelengths. The aperture radii are typically several times larger than the beam sizes for wavelengths $\leq$ 70 $\mu$m (and by greater factors for the fixed aperture method that uses the 70 $\mu$m aperture radii across all bands). At wavelengths $>$ 70 $\mu$m, the fixed aperture radius set at 70 $\mu$m is always used, and the aperture diameter is still usually larger than the image resolution (except for G305.20 whose fixed aperture diameter becomes similar to the resolution at the longest wavelength 500 $\mu$m).

3.2 Zhang & Tan Radiative Transfer Models

We use Zhang & Tan (2018, [ZT18]) radiative transfer (RT) models (hereafter ZT models) to fit the SEDs and derive key physical parameters of the protostars. In a series of papers, Zhang & Tan (2011), Zhang et al. (2013b), Zhang et al. (2014) and ZT18 have developed models for the evolution of high- and intermediate-mass protostars based on the Turbulent Core Model (MT03). In this model, massive stars are formed from pre-assembled massive pre-stellar cores, supported by internal pressure that is provided by a combination of turbulence and magnetic fields. With various analytic or semi-analytic solutions, they calculate the properties of a protostellar core with different components, including the protostar, disk, infall envelope, outflow, and their evolutions, self-consistently from given initial conditions. The main free parameters in this model grid are: the initial mass of the core $M_{c}$; the mass surface density of the clump that the core is embedded in $\Sigma_{\rm cl}$; and the protostellar mass, $m_{*}$, which indicates the evolutionary stage. In addition, there are secondary parameters of inclination angle of line of sight to the outflow axis, $\theta_{\rm view}$, and the level of foreground extinction, $A_{V}$.

The evolutionary history of a protostar from a given set of initial conditions ($M_{c}$ and $\Sigma_{\rm cl}$) is referred to as an evolutionary track, and a particular moment on such a track is a specified $m_{*}$. Therefore the model grid is of three dimensions ($M_{c}$-$\Sigma_{\rm cl}$-$m_{*}$), including the entire set of tracks. Currently, $M_{c}$ is sampled at 10, 20, 30, 40, 50, 60, 80, 100, 120, 160, 200, 240, 320, 400, 480 $M_{\odot}$, $\Sigma_{\rm cl}$ is sampled at 0.1, 0.32, 1, 3.2 g $\rm cm^{-2}$, forming 60 evolutionary tracks. Then $m_{*}$ is sampled at 0.5, 1, 2, 4, 8, 12, 16, 24, 32, 48, 64, 96, 128, 160 $M_{\odot}$. Note that not all of these $m_{*}$ are sampled for each track. In particular, the maximum protostellar mass is limited by the final stellar mass achieved in a given evolutionary track. As a result, there are 432 different physical models defined by different sets of $M_{c}$, $\Sigma_{\rm cl}$ and $m_{*}$.

There are several things to note about the models. First, the models describe one protostar forming through monolithic collapse from the parent core. The formation of binary and multiple systems is not included in the models. Second, compared with the Robitaille et al. (2007) RT models that mostly focus on lower-mass protostars, the ZT18 model grid has broader parameter space relevant to high pressure, high density and thus high accretion rate conditions of massive star formation, while keeping the number of free parameters low. Third, the models do not explicitly include the clump component, which contributes to foreground extinction at short wavelengths and additional emission at long wavelengths. The former effect is compensated for by the free parameter $A_{V}$. The latter effect requires the model grid fitting to be done on clump-envelope-background-subtracted SEDs. Fourth, the aperture scale for the measured SED is not considered in the fitting process. The predicted SEDs in the model grid are total SEDs, which include modest contributions from parts of the outflow that extend beyond the core. We assume with the aperture adopted we also measure the total emission from the protostar and ideally the models that describe that observed SED best would predict a similar scale (this can be checked after the fitting results are returned). Fifth, PAH emission and thermal emission from transiently (single-photon) heated very small grains at $\lesssim$ 8 $\mu$m is not modeled, and so our method is to use the SEDs at these wavelengths as upper limits. Lastly, while the general trends of the features of the SEDs are determined by the initial/environmental conditions and evolution, some detailed features, such as the peak wavelength and long-wavelength spectral index, may be affected by the particular dust models used in the radiative transfer simulations.

3.3 SED Fitting

When fitting the SEDs to the models, we use our fiducial case, i.e., using fixed aperture size for all wavelengths, and set data points at wavelengths $\leq$ 8 $\mu$m as upper limits since the effects of PAH emission and thermal emission from very small grains are not included in the ZT RT models. For G309, the Spitzer 4.5 $\mu$m, 5.8 $\mu$m and 8 $\mu$m data have a ghosting problem. For G45.12 and IRAS 16562 all Spitzer data have ghosting problems. Thus we do not use these data for the SED fitting. The error bars are set to be the larger of either 10% of the clump background-subtracted flux density to account for calibration error, or the value of the estimated clump background flux density (see §3.1), which is used for background subtraction, given that order unity fluctuations in the surrounding background flux are often seen.

The fitting procedure involves convolving model SEDs with the filter response functions for the various telescope bands. Source distances are adopted from the literature. For each source, we present the five best-fitting models. Again we note that the SED model fitting performed here assumes that there is a single dominant source of luminosity, i.e., effects of multiple sources, including unresolved binaries, are not accounted for. This is a general limitation and caveat associated with this method as discussed in Paper I.

4 Results

The types of multi-wavelength data available for each source, the flux densities derived, and the aperture sizes adopted are listed in Table 2. $F_{\lambda,\rm fix}$ is the flux density derived with a fixed aperture size and $F_{\lambda,\rm var}$ is the flux density derived with a variable aperture size. The value of flux density listed in the upper row of each source is derived with background subtraction, while that derived without background subtraction is listed in brackets in the lower row. The SOFIA images for each source are presented in §4.1. General results of the SOFIA imaging are summarized in §4.2. The SEDs and fitting results are presented in §4.3.

[FIGURE:]

4.1 Description of Individual Sources

In this section we describe the MIR morphology of each source and also try to identify the nature of the structures revealed by our SOFIA or Gemini imaging, together with archival NIR data and other data from the literature.

4.1.1 G45.12+0.13

This UC HII region, also known as IRAS 19111+1048, has a measured far kinematic distance of 7.4 kpc (Ginsburg et al. 2011). The radio morphology of this region shows a highly inhomogeneous ionized medium (Vig et al. 2006), which is consistent with the extended MIR morphology revealed here in Figure 1. Vig et al. (2006) proposed the source is an embedded cluster of Zero Age Main Sequence (ZAMS) stars with twenty compact sources, including one non-thermal source, identified by their radio emission. The central UC HII source S14 is deduced to be of spectral type O6 from the integrated radio emission. They also found there are two NIR objects, IR4 and IR5, within the S14 region, while IR4 is at the peak of the radio emission and matches the OH maser position obtained by Argon et al. (2000). We see that most sources revealed at 8$\mu$m and 37$\mu$m in the central region have counterparts in NIR bands (see Figure 12), which also indicates that this site is probably a protocluster.

An extended bipolar outflow is revealed in CO(2–1), CO(3–2), CO(6–5), 13CO(2–1) and C18O(2–1) by Hunter et al. (1997). Higher resolution 13CO(1–0) observations resolve the system into at least two outflows. The highest velocity outflow appears centered on the UC HII region S14. The additional bipolar outflow was identified with a dynamical center lying offset (-8$\arcsec$, -3$\arcsec$) from S14, named “G45.12+0.13 west” by Hunter et al. (1997). Hunter et al. (1997) argued that G45.12+0.13 west most likely represents dust emission from a younger or lower-mass protostar that formed during the same epoch as the ionizing star of S14. They also argued the absence of H2O masers in the G45.12+0.13 cloud core suggests that both of the outflow sources have evolved beyond the H2O maser phase.

In our SOFIA images we see MIR to FIR emission peaking at the S14 position. We do not see a distinct source at the position of G45.12+0.13 west, though the MIR extension to the southwest of S14 could be due to the two blue-shifted outflows, which are also revealed in NIR (see Figure 12). There is a MIR peak $\sim$ 7.7″to the SE of S14, which is best revealed at 19 $\mu$m and further down $\sim$ 22″to the SW of S14 there is another MIR peak. The closer one is seen in all J, H, K bands while the further one is seen in H and K bands as shown in Figure 12. They could be more evolved low-mass protostars.

4.1.2 G309.92+0.48

This region is located at a distance of 5.5 kpc (Murphy et al. 2010). The MIR emission in this area was resolved into 3 sources with the CTIO 4-m at 10.8 $\mu$m and 18.2 $\mu$m, labeled 1 through 3 (see Figure 2 in De Buizer et al. 2000). In addition to these sources, our Gemini 11.7 $\mu$m data also shows three additional fainter point-sources, as shown in Figure 3, which we label 4 through 6. Note that all the sources that appear in the Gemini field in Figure 3 are located within the northern patch revealed by SOFIA 7.7 $\mu$m in Figure 2.

Source 1 is the brightest source in the MIR and is coincident with a cm radio continuum source believed to be a HC HII region (Phillips et al. 1998; Murphy et al. 2010). Our Gemini 11.7 $\mu$m image resolves Source 1 into two components as shown in Figure 3, which we name 1N and 1S. Since both sources are elongated at the same position angle, it may be that the dark lane between them is an area of higher obscuration. In fact, the radio continuum emission at 8.6 GHz (Walsh et al. 1998, Philips et al. 1998) and 19 GHz (Murphy et al. 2010) towards Source 1 shows a peak nearly in between mid-infrared Source 1N and 1S, possibly tracing the location of the highly embedded protostar. Both of the radio observations of Philips et al. (1998) and Murphy et al. (2010) show elongation in the same direction as the MIR-dark lane. However, in both cases the beam profile is also elongated in the same direction. The 8.6 GHz observations of Walsh et al. (1998) have similar resolution and a nearly circular beam, and do not show any elongation.

OH and Class II methanol masers are found to be distributed along an arc centered near the primary radio continuum peak (see inset in Figure 3) with increasingly negative line-of-sight velocities from north to south (Caswell 1997). Norris et al. (1993) considered this site to have a well-defined methanol maser velocity gradient and forwarded the idea that they are tracing a near-edge-on circumstellar disk. The MIR morphology seen in the Gemini data do not appear to support this idea. If the dark lane between elongated Sources 1N and 1S is indeed the location of the protostar as the radio peak suggests, then the morphology at 11.7 $\mu$m would be best explained as the emission from the walls of outflow cavities or flared disk surfaces, with the dark lane representing a nearly edge-on, optically-thick (in the IR), circumstellar disk. This disk plane would be perpendicular to the methanol maser distribution. Thus the Class II methanol masers may be coming from a region which experiences both strong shocks, but also a strong radiation field, which enables radiative pumping of the masers. To help infer the outflow orientation, De Buizer (2003) observed the field for signs of H2 emission, however, none was detected (note, however, that this H2 survey was relatively shallow). We could not find any additional outflow information about this region. Note that the extension of the central source in the Spitzer 8 $\mu$m image and the stripes in the SOFIA 31 $\mu$m and 37 $\mu$m images in Figure 2 are artifact features caused by very bright sources on the array.

With the NIR VVV data, we find there is little to no NIR emission from 1N, which suggests that it is the most obscured source seen in the MIR. In the J-band there is a compact emission source $\sim$2$\arcsec$ NE of the peak of Source 1N in the direction of Source 2, but no emission directly coming from Source 1N or 1S. The H-band image shows a source in this same location, but with the addition of an extended source with a peak coincident with 1S, and a “tail” to the SE. At Ks, there is only an extended source with a peak at 1S, and extended emission in the same direction as the tail seen in H-band, with emission also extending NE towards 1N. Source 2 lies to the northeast of Source 1 at a position angle of 53$\arcdeg$. Both Source 1 and 2 are seen at 8.6 GHz by Phillips et al. (1998) and in the NIR by Walsh et al. (1999). With the NIR VVV data, we find that Sources 2 and 3 are also seen at J, H and Ks. Source 6 is also seen at J, H, and Ks, Source 4 is seen at H and Ks, but Source 5 is not detected in the NIR. In our 7.7 $\mu$m SOFIA data we see fingers of emission reaching the area around Sources 3 and 5, as well as Source 6, though these are not detected at longer wavelengths in the SOFIA data.

In the larger field of view of the SOFIA data, we detect another extended (r$\sim$5$\arcsec$) emission region $\sim$18$\arcsec$ south of Source 1 at all SOFIA wavelengths. The nature of this region is unknown, however.

4.1.3 G35.58-0.03

The star-forming region G35.58-0.03 is located at the far kinematic distance of 10.2 kpc (Fish et al. 2003; Watson et al. 2003). Kurtz et al. 1994 resolved the 2 cm and 3.6 cm continuum emission here into two UC HII regions $\sim$2$\arcsec$ apart, with the western source named G35.578-0.030 and the eastern source named G35.578-0.031. G35.578-0.030 contains water and OH masers, but no methanol masers (Caswell et al. 1995). Zhang et al. (2014) found that there is an ammonia clump peaked co-spatially with their observed 1.3 cm radio continuum peak, which is $\sim$0.4$\arcsec$ north of the 2 cm peak of G35.578-0.030 (Kurtz et al. 1994; 1999). H30$\alpha$ shows evidence of an ionized outflow connecting to a molecular outflow seemingly centered on the radio continuum peak of G35.578-0.030. Only faint 1.3 cm continuum emission was found from the eastern source, G35.578-0.031, and no signs of outflow or ammonia emission.

De Buizer et al. (2005) presented $\sim$0.6$\arcsec$ resolution MIR images of this region at 10 and 20 $\mu$m, which showed a single source with some extension to the northwest. Due to poor astrometry of the data, it was unclear which UC HII region the mid-infrared emission was associated with. They argued that, due to the fact that the western source, G35.578-0.030, appears to have a similar extension to the northwest at 3.6 cm as seen by Kurtz et al. (1999), that the MIR emission is likely to be associated with that source.

Our data obtained at 11.7 $\mu$m from Gemini with $\sim$0.3$\arcsec$ resolution further resolve the MIR emission into a main bright peak with two fingers of extended diffuse emission to the north and northwest. Using Spitzer 8 $\mu$m images to confirm our astrometry, it is revealed that the MIR peak is not associated with the western UC HII region, but instead the eastern UC HII region, G35.578-0.031 (see Figure 5). The relative astrometric error between the Gemini 11.7 $\mu$m image and the radio data is better than 0.3$\arcsec$. No MIR emission is detected at the location of G35.578-0.030 out to 37 $\mu$m. The MIR peak is, however, close to the location of the redshifted outflow cavity of G35.578-0.030 seen in CO(2–1) by Zhang et al. (2014). However, if high extinction was causing the general lack of MIR emission from G35.578-0.030, it seems unlikely that the MIR emission we are seeing would come from the even more extinguished red-shifted outflow cavity of G35.578-0.030. It is more plausible that the MIR emission is coming solely from the eastern UC HII region, G35.578-0.031.

Our SOFIA images of this region (Figure 4) show a bright source peaked at the location of G35.578-0.031 and extended slightly to the northwest, as is seen in the higher spatial resolution Gemini 11.7 $\mu$m image Figure 5. The nature of this extension is unclear, since the outflow seen by Zhang et al. (2014) has an axis oriented east-west. A second compact source is detected in our SOFIA data (and in the Spitzer-IRAC data) located $\sim$10$\arcsec$ to the east of G35.578-0.031. There is also a hint of MIR extension to the west, which may be due to the outflow.

The eastern MIR source seen in the SOFIA data has a counterpart at K-band as can be seen from Figure 12. Thus it may be a more evolved protostar, closer to the end of its accretion. From the NIR image (see Figure 12) there are at least two K-band sources within the highest contour of the 37 $\mu$m emission. The southern K-band source is associated with the peak at 8 $\mu$m and the main bright peak at 11.7 $\mu$m, while the northern K-band source has some overlap with the northern finger in Gemini 11.7 $\mu$m image (not shown here). There could be one or two lower luminosity companion sources in that region together with the southern main massive protostar, but they are not well resolved in the MIR and FIR.

4.1.4 IRAS 16562-3959

This source (also known as G345.49+1.47) is located at a distance of 1.7 kpc (Guzmán et al. 2010). It is believed that the massive core hosts a high-mass star in an early stage of evolution, including ejection of a powerful collimated outflow (Guzmán et al. 2010). Guzmán et al. (2010) carried out ATCA observations to reveal five 6 cm radio sources: a compact bright central (C) component, two inner lobes that are separated by about 7$\arcsec$ and symmetrically offset from the central source, and two outer lobes that are separated by about 45$\arcsec$ (see Figure 4 in Guzman et al. 2010). The central radio source has a 3 mm counterpart, source 10 in Guzmán et al. (2014), and an X-ray counterpart, source 161 in Montes et al. (2018), and is associated with OH maser emission (Caswell 1998, 2004). It is interpreted as a HC HII region based on hydrogen recombination line (HRL) observations (Guzmán et al. 2014). The continuum at 218 GHz and CH3CN(12–11) (methylcyanide) observations by Cesaroni et al. (2017) revealed that the central source 10 actually consists of two peaks. The four other symmetrically displaced sources are interpreted as shock-ionized lobes (Guzmán et al. 2010) and are observed to move away from the central source at high speed (Guzmán et al. 2014).

On the other hand, the molecular observations of CO(6–5) and CO(7–6) show the presence of high-velocity gas exhibiting a quadrupolar morphology (Guzmán et al. 2011), most likely produced by the presence of two collimated outflows, one major outflow lying with a southeast-northwest (SE-NW) orientation, and the other with a N-S orientation, which may come from the unresolved mm source 13 in Guzmán et al. (2014) to the east of the central source. The SE-NW molecular outflow is aligned with the string of radio continuum sources. Extended Ks-band emission probably tracing excited $\rm H_{2}$-2.12 $\mu$m is also associated with the SE-NW flow.

In Guzmán et al. (2014), the molecular core in which the outflow is embedded presents evidence of being in gravitational contraction as shown by the blue asymmetric peak seen in HCO+(4–3). The emission in the SO2, 34SO, and SO lines exhibits velocity gradients interpreted as arising from a rotating compact ($\sim$ 3000 AU) molecular core with angular momentum aligned with the jet axis. López-Calderón et al. (2016) reported 13CO(3–2) APEX observations of this region and showed that the high-mass protostellar candidate is located at the column density maximum. Montes et al. (2018) decomposed the wider region into 11 subclusters with results from Chandra X-ray observations together with VISTA/VVV and Spitzer-GLIMPSE catalogs and the subcluster containing the high-mass protostar was found to be the densest and the youngest in the region with the high-mass protostar located near its center.

In our MIR images, the extended IR emission is likely tracing the illuminated inner outflow cavity containing the jet. There are two knots to the northeast of the central source revealed by SOFIA. The closer knot located $\sim$15$\arcsec$ NE of the central source is associated with the 92.3 GHz peak 18 in Guzmán et al. (2014), as well as a K-band source (see Figure 12). It may correspond to the X-ray source 178 in Montes et al. (2018). There is OH maser emission (Caswell 1998, 2004), but no radio continuum emission detected. Thus it may be a low-mass protostar. The farther knot, located $\sim$36$\arcsec$ northeast of the central source, has counterparts in all of the J, H, K bands. We did not find any associated X-ray source for this knot in the Montes et al. (2018) sample.

4.1.5 G305.20+0.21

G305.20+0.21 is a massive star-forming region located at a distance of $4.1^{+1.2}_{-0.7}$ kpc from parallax of 6.7 GHz methanol masers (Krishnan 2017). Class II methanol (CH3OH) masers were reported in two positions by Norris et al. (1993): G305.21+0.21 and G305.20+0.21 separated by approximately 22$\arcsec$. Walsh & Burton (2006) refer to these maser sites as G305A and G305B, respectively, and we will adopt that nomenclature here.

The brightest MIR source appears to be associated with the methanol masers of G305B, but does not possess detectable radio continuum emission (below a 4$\sigma$ detection limit of 0.9 mJy beam*-1* (beam $\sim$ 1.5$\arcsec$) at 8.6 GHz in Phillips et al. 1998, and a 3$\sigma$ detection limit of 0.09 mJy at 18 GHz in Walsh et al. 2007). Walsh et al. (2007) found no HC3N, NH3, OCS, or water at the position of G305B and proposed that it has evolved enough to the point that it has already had time to clear out its surrounding molecular material. By contrast, Boley et al. (2013) proposed G305B is a massive protostar in a pre-UCHII-region stage. Our SOFIA images show that G305B is the brightest MIR source out to 37 $\mu$m. Our high-spatial-resolution Gemini data (Figure 8) show G305B is resolved into two emission components, with the fainter secondary source (which we name G305B2) lying $\sim$1$\arcsec$ to the NE of the brighter source (G305B1). G305B2 is only visible at wavelengths greater than 8.8 $\mu$m. By contrast, G305B1 is seen to have emission in all Gemini images from 3.8 to 24.5 $\mu$m, and has a NIR counterpart as well (see Figure 12 and Walsh et al. 1999). Using four infrared sources seen in both the Gemini 3.8 $\mu$m image (but not shown in Figure 3) and the Spitzer 3.6$\mu$m image we were able to confirm the absolute astrometry of the Gemini data at all wavelengths to better than 0.2$\arcsec$. This places the Class II methanol maser reference feature (i.e., the brightest maser spot) from Phillips et al. (1998) $\sim$0.5″NE of the MIR peak (see the 9.7$\mu$m image in Figure 8). It is not clear what these masers are tracing.

What is the nature of the MIR double source associated with G305B? G305B2 could be a more embedded source, since it is not visible at shorter IR wavelengths. However, it appears to change shape considerably as a function of wavelength, flattening and becoming more diffuse at 18.3 and 24.5 $\mu$m. G305B1 also changes shape modestly with wavelength and its shape at 9.7 and 10.4 $\mu$m is peculiar. The northeast side of G305B1 is very flat, and almost completely straight at 9.7 and 10.4 $\mu$m (see white line in the 9.7$\mu$m panel of Figure 8 as reference). As these filters are sampling the wavelength of peak dust extinction (Gao et al. 2009), it may be that the morphologies of both sources could be explained if the dark lane between them is a “silhouette” of a circumstellar disk or toroid that is optically thick in the MIR. The brighter MIR source G305B1 would be the side of the disk or outflow cavity facing towards us, and G305B2 the side facing away which we only see at longer wavelengths due to extinction from the disk along the line of sight. We could corroborate the outflow cavity hypothesis if we had evidence of an outflow and knew its angle. Walsh et al. (2006) did image the area in commonly used outflow tracers 13CO and HCO+, and presented the data as integrated emission maps. However, the emission appears to peak on G305A and extends at larger scales in a direction parallel to the dark lane orientation, tracing the location of the extended 1.2 mm continuum emission (rather than an outflow). However, if the hypothesis of Walsh et al. (2007) is correct, i.e., that due to low chemical abundance this source is more evolved and has cleared much of its surrounding molecular material, then the source may have passed the stage where it would exhibit an active outflow. Conversely, a Class I methanol maser was detected by Walsh et al. (2007) 3$\arcsec$ due east of G305B, and they are generally only found in outflows.

Walsh et al. (2001) observed the 6.7 GHz methanol maser site G305A in the MIR (10.5 $\mu$m and 20 $\mu$m) and found that G305A is not associated with any MIR source. G305A is out of the field of our Gemini images. However, we see strong emission from G305A in our SOFIA images at 19 $\mu$m and longer, and it becomes the dominant source in the FIR starting at Herschel 70 $\mu$m. G305A is also not associated with any 8.6 GHz continuum emission with a flux density limit of 0.55 mJy beam*-1* (Phillips et al. 1998) or 18 GHz continuum emission with a detection limit of 0.15 mJy (Walsh et al. 2007), but is rich in molecular tracers (Walsh et al. 2007) indicating it is a source that is likely much younger and more embedded than G305B and in a hot core phase, prior to the onset of a UC HII region.

About 15$\arcsec$ to the southwest of G305B is an extended HII region, G305HII, with a flux of 195 mJy at 8.6 GHz (Phillips et al. 1998). We detect this source in all of our SOFIA images. We also detect an infrared source between G305A and G305B, which we call G305C, located $\sim$ 14″east of G305B. It is present at all wavelengths in the SOFIA images, but becomes less pronounced at longer wavelengths. It also has NIR counterparts, as shown in Figure 12, which seem to be resolved into three peaks. The nature of the source is uncertain, but it may be a low mass YSO. Besides the G305HII region there is no other radio emission in the field shown in Figure 7 revealed by the 18 GHz continuum in Walsh et al. (2007).

4.1.6 G49.27-0.34

This source, classed as an “extended green object” (EGO) is in an IRDC with near kinematic distance of 5.55 $\pm$ 1.66 kpc (Cyganowski et al. 2009). The MIR peak (see Figure 9) is associated with the 3.6 cm radio source CM2 in Cyganowski et al. (2011a). Towner et al. (2017) did not detect a 1.3 cm counterpart to CM2 at the a $4\sigma$ detection limit of 0.28 mJy beam*-1* (beam $\sim$ 1″). The MIR extension to the northeast is associated with a stronger radio source CM1 detected at 3.6 cm and 1.3 cm by Cyganowski et al. (2011a) and at 20 cm by Mehringer (1994).

We did not find any outflow information about this source. De Buizer & Vacca (2010) obtained Gemini L- and M-band spectra for this EGO, and detected only continuum emission (no H2 or CO). However, Cyganowski et al. (2011a) suspected that an outflow, perhaps driven by CM2 or by a massive protostar undetected at cm wavelengths, may exist, but is not detected, given the 44 GHz Class I CH3OH masers and 4.5 $\mu$m emission in the south. SiO(5–4), HCO+ and H13CO+ emission is detected toward this EGO with JCMT (Cyganowski et al. 2009). No 6.7 GHz CH3OH emission is detected towards this EGO (Cyganowski et al. 2009). Neither thermal nor maser 25 GHz CH3OH emission is detected (Towner et al. 2017).

There is a secondary component revealed by our SOFIA data to the south of the main MIR peak. It is also seen at 3.6 cm (Cyganowski et al. 2011a) but not at 1.3 cm (Towner et al. 2017). The nature of this source is unknown. We do not see obvious counterparts in the NIR image (see Figure 12).

4.1.7 G339.88-1.26

This source, also named IRAS 16484-4603 is located at $2.1^{+0.4}_{-0.3}$ kpc, determined from trigonometric parallax measurements of the 6.7 GHz methanol masers using the Australian Long Baseline Array (Krishnan et al. 2015).

De Buizer et al. (2002) resolved the central MIR emission of G339.88 into three peaks (1A, 1B, and 1C) at 10 and 18 $\mu$m that all lie within an extended MIR region elongated at a position angle of $\sim$120$\arcdeg$ (Figure 11a). Interferometric radio continuum observations have revealed an elongated, ionized jet/outflow at a position angle of $\sim$45$\arcdeg$ with a scale of 15″, approximately perpendicular to the elongation of the infrared emission (Ellingsen et al. 1996; Purser et al. 2016). Recent ALMA 12CO(2–1) observations by Zhang et al. (2019) also reveal a major molecular outflow with a E-W orientation and a tentative second outflow with a NE-SW orientation (at the same angle as the ionized outflow seen by Purser et al. 2016). Zhang et al. (2019) suggest that the 1.3 mm continuum peak, which is $\sim$ 0.5″ to the west of 1B, is the likely location of the origin of both outflows, which may indicate an unresolved proto-binary system. All of the 10 and 18 $\mu$m MIR emission seen by De Buizer et al. (2002) is therefore mainly tracing the outflow cavities of the molecular outflow seen at a position angle of $\sim$120$\arcdeg$ (Figure 11a). Our SOFIA data (see Figure 10) show an extension in this direction as well, seen best at 19.7 $\mu$m. At wavelengths longer than 20 $\mu$m, there is a faint pull of emission to the NE and another faint extension to the SW, both of which correspond to the radio lobes of the ionized outflow (Figure 11b). Therefore, both outflows are revealed in the IR, with the ionized outflow only showing up at longer wavelengths, which again may be due to extinction. Detection of red and blue-shifted emission on both sides suggests a near side-on view of the outflows.

There is a large half-moon feature to the east of the main MIR peak in our SOFIA data. It has radio continuum emission (see Ellingsen et al. 2005) and could be a cometary compact HII region. Closer to the main MIR peak, we also see a secondary source $\sim$ 10″to the south. There is no CO outflow associated with this source. We see a counterpart of this source in H and K band as seen in Figure 12. It could be a more evolved low-mass protostar. The source that is further SW, which is getting stronger at wavelengths longer than 31 $\mu$m, might be related to the ionized radio jet (Purser et al. 2016), though there is no hint of ionized emission from it in the study of Ellingsen et al. (2005).

4.2 General Results from the SOFIA Imaging

Overall in the sample of sources we have studied here, we often see that the MIR morphologies appear to be influenced by the presence of outflow cavities, which create regions of low dust extinction, and the presence of relatively cool, dense gas structures (potentially including disks and infall envelopes), which have high dust extinction, even at relatively long wavelengths. The presence of such structures is a general feature of Core Accretion models. A number of sources also appear to have companions, including from being in regions where a star cluster is likely forming, which can also complicate the appearance in the MIR.

In addition to the monochromatic images presented above, we also construct three-color images of all the sources, presented together in Figure 13. Note, however, that these RGB images have different beam sizes for the different colors (especially blue), with the effect being to tend to give small sources an extended red halo. In spite of this effect, in G45.12, G309.92, G35.58, IRAS 16562 and G339.88, short wavelength emission seems to dominate the extended structure. In IRAS 16562, we can see the near-facing outflow cavity appears bluer while the more extincted, far-facing outflow cavity appears redder. For the other sources we do not see obvious color gradients across the sources.

We summarize the properties of the protostellar sources in Table 3. The ordering of the sources is from high to low for the luminosity estimate (top to bottom). For two out of the three sources with detected outflows, the MIR morphology is significantly influenced by outflow cavities. For those lacking outflow information, we consider that it is still likely that the MIR emission is tracing outflows or flared disks. Especially in G309.92 and G305.20, high-resolution Gemini data reveals a flat dark lane, which could be the optically thick disk.

We see that at wavelengths $\gtrsim$ 19 $\mu$m, there is an offset between the radio emission, if that is where the protostar is located, and the MIR peaks in G309.92, G35.58, G49.27, G339.88. In Paper I we found that the MIR peaks appear displaced away from the protostar towards the blue-shifted, near-facing side of the outflow due to the higher extinction of the far-facing side at short wavelength. Here G339.88 may reveal a hint of this trend of the displacement. For the other sources, due to the lack of outflow information, the cause of the offset is not yet clear.

We have also found candidates of more evolved, probably lower-mass protostars in the company of the massive protostar in most regions, based on the common peaks seen in multi-wavelength MIR and NIR data and how their fluxes change with wavelength. With the caveat that our sample is likely incomplete, the projected separation between the massive protostar and the nearest lower-mass companion revealed by SOFIA is about 0.28 pc in G45.12, 0.49 pc in G35.58, 0.12 pc in IRAS16562, 0.28 pc in G305.20, and 0.10 pc in G339.88. Note that Core Accretion models, such as the Turbulent Core Model of McKee & Tan (2003), can be applied to conditions inside protoclusters, as well as to more isolated regions, while Competitive Accretion (Bonnell et al. 2001) and Protostellar Collision (Bonnell et al. 1998) models require the presence of a rich stellar cluster around the protostar. To the extent that some of the presented sources appear to be in relatively isolated environments is thus tentative evidence in support of Core Accretion models, but deeper observations to probe the low-mass stellar population are needed to confirm this.

4.3 Results of SED Model Fitting

4.3.1 The SEDs

Figure 14 shows the SEDs of the seven sources that have been discussed in this paper. The figure illustrates the effects of using fixed or variable apertures, as well as the effect of background subtraction. Our fiducial method is that with fixed aperture and with background subtraction carried out. This tends to have moderately larger fluxes at shorter wavelengths than the variable aperture SED especially for G35.58, IRAS16562 and G339.88 where emission from secondary sources can be significant at wavelengths $\leq$ 8 $\mu$m. However, as in Paper I, the $\leq$ 8 $\mu$m flux is in any case treated as an upper limit in the SED model fitting, given the difficulties of modeling emission from PAHs and transiently heated small grains. The flux density derived from the two methods between 10 $\mu$m and 70 $\mu$m is generally close. For flux densities longer than 70 $\mu$m, the influence of secondary sources is not illustrated by the variable aperture method. However, we tried measuring the SEDs up to 37 $\mu$m of the MIR companions alone, which are resolved from the emission of the main protostar, and found that their flux density at each wavelength is $\leq$ 5% of that of the main protostar (except that the 19 $\mu$m flux density of the southern patch in G49.27 is $\sim$ 20% of that of the massive protostar). Moreover, all of them have a SED peak $\leq$ 31 $\mu$m except that the southern patch in G49.27 has a nearly flat rising slope between 31 $\mu$m and 37 $\mu$m. Thus the influence of secondary sources is generally not severe at long wavelengths that control the SED fitting.

Again, as mentioned in §3.1, for the cases where there seem to be multiple sources in the fiducial aperture, the model assumes one source dominates the luminosity and the key is to measure the flux from the same region across all wavelengths. If a source is isolated, then the fixed aperture at shorter wavelengths, which tends to be larger than the source appears, may include more noise and make the photometry less accurate than the variable aperture method. However, since we set the clump background emission as the magnitude of the uncertainty, this effect should be very minor.

The peaks of the SEDs are generally between 37 $\mu$m and 70 $\mu$m. In particular, the SED peaks of G45.12, G309.92, G305.20 appear to be closer to 37 $\mu$m, while the peaks of G35.58, G49.27 and G339.88 appear to be closer to 70 $\mu$m. This may be related to the evolutionary stage and/or viewing angle of the sources (see §4.3.2).

4.3.2 ZT Model Fitting Results

Figure 15 shows the results of fitting the ZT protostellar radiative transfer models to the fixed aperture, background-subtracted SEDs. Note that the data at $\leq$ 8 $\mu$m are considered to be upper limits given that PAH emission and transiently heated small grain emission are not well treated in the models. The parameters of the best-fit ZT models are listed in Table 4. From left to right the parameters are reduced $\chi^{2}$ (i.e., normalized by the number of data points in the SED, $N$), the initial core mass ($M_{c}$), the mean mass surface density of the clump ($\Sigma_{\rm cl}$), the initial core radius ($R_{\rm core}$), the current protostellar mass ($m_{*}$), the viewing angle ($\theta_{\rm view}$), foreground extinction ($A_{V}$), current envelope mass ($M_{\rm env}$), half opening angle of the outflow cavity ($\theta_{w,\rm esc}$), accretion rate from the disk to the protostar ($\dot{m}_{*}$), the luminosity integrated from the unextincted model SEDs assuming isotropic radiation ($L_{\rm bol,iso}$), and the inclination corrected, true bolometric luminosity ($L_{\rm bol}$). For each source, the best five models are shown, ordered from best to worst as measured by $\chi^{2}$. Note that these are distinct physical models with differing values of $M_{c}$, $\Sigma_{\rm cl}$, and/or $m_{\rm*}$, i.e., we do not display simple variations of $\theta_{\rm view}$ or $A_{V}$ for each of these different physical models.

The best-fit models imply the sources have protostellar masses $m_{*}\sim 12-64\>M_{\odot}$ accreting at rates of $\dot{m}_{*}\sim 10^{-4}-10^{-3}\>M_{\odot}\>\rm yr^{-1}$ inside cores of initial masses $M_{c}\sim 100-500\>M_{\odot}$ embedded in clumps with mass surface densities $\Sigma_{\rm cl}\sim 0.1-3\>\rm g\ cm^{-2}$ and span a luminosity range of $10^{4}-10^{6}\>L_{\odot}$.

In most sources the best five models have similar values of $\chi^{2}$, but there is still significant variations in the model parameters even for G305.20 which has the most SED data points. As stated in Paper I, this illustrates the degeneracy in trying to constrain the protostellar properties from only MIR to FIR SEDs, which would be improved by extended SEDs fitting including centimeter continuum flux densities (Rosero et al. 2019) and image intensity profile comparison. From the SED shape the most variation between models appears at shorter wavelengths. Here more data points can help better constrain the models, as in G305.20. Again we note that although sometimes the $\chi^{2}$ may look high, as in G45.12+0.13, here we focus more on the relative comparison of $\chi^{2}$ between the models available in the model grid, which still give us constraints on the protostellar properties. At wavelengths $>\ 70\ \mu$m the models tend to be lower than the data points in many sources. Note the values of $R_{\rm core}$ returned by the models are usually much smaller than the aperture radii. This would indicate that, even after a first attempt at clump background subtraction, the measured flux still has significant contribution from the cool surrounding clump. Recall that this component is not included in the ZT radiative transfer models and can thus lead to the offset at long wavelengths, i.e., with models under-predicting the observed fluxes.

We also tried fitting the SEDs with variable apertures across wavelength. Most sources have $R_{\rm core}$ similar to that derived in the fiducial case and still the models appear lower than the data points at long wavelengths for G309.92, G35.58, G305.20 and G49.27.

We note that $m_{*}$ appears quite high for G45.12+0.13, G309.92+0.48, G35.58-0.03, IRAS 16562, G305.20+0.21. This is likely due to there being more than one protostar inside the aperture, even with the variable aperture case, like the source G35.20-0.74 in Paper I, where the stellar mass returned by the models is around the sum of the two binary protostars in the center (Beltrán et al. 2016, Zhang et al. 2018).

The location of the SED peak is thought to show a dependence on the evolutionary stage of the source. We compare the current age derived from the models and the corresponding total star formation time scale based on Eq. (44) in MT03 assuming a star formation efficiency of 0.5. G305.20 appears to be the most evolved followed by G309.92 and G45.12. G339.88 appears to be the least evolved followed by G49.27. G339.88 is still deep embedded with high dust extinction while G49.27 is an IRDC source. They are likely the youngest YSOs among the seven sources. The evolutionary stage revealed by the models is consistent with the picture that more evolved sources have a SED peak located at shorter wavelengths, as described in §4.3.1. However, orientation effects may also be playing a role, since the peak of the SED shifts to shorter wavelengths when viewing sources at angles closer to their outflow axis.

Next we describe the fitting results of each individual source and compare with previous literature results.

G45.12+0.13: This is our most luminous source (almost $10^{6}\>L_{\odot}$) and hits the boundary of the parameter space of the ZT model grid, which is partly why the models do not seem to fit the data points very well, as shown in Figure 15, since there are only a few models around $10^{6}\>L_{\odot}$ (Zhang & Tan 2018). As an experiment, we tried changing the distance from 7.4 kpc to 1 kpc and were able to obtain fitting results that have much smaller values of $\chi^{2}$. On the other hand, this region is likely to be a protocluster hosting many ZAMS stars. Thus the assumption of one source dominating the luminosity may not work well here. The current best models indicate high initial core mass $M_{c}\sim 500\>M_{\odot}$, high $\Sigma_{\rm cl}\gtrsim 1.0\rm\ g\>cm^{-2}$ clump environment and high protostellar mass $m_{*}\geq 24\>M_{\odot}$ for the dominant source. The accretion rate is $\sim 10^{-3}M_{\odot}\>{\rm yr}^{-1}$. The current envelope mass is also typically as high as $\sim$ 400 $M_{\odot}$. The foreground extinction $A_{V}$ is estimated to be very low, but this may be an artefact of other problems of the model fitting. The best five models all give a $\theta_{\rm view}$ close to $\theta_{w,\rm esc}$, which leads to high levels of short wavelength emission.

G309.92+0.48: The best models prefer a massive protostar of $\sim$ 24 to 48 $M_{\odot}$ accreting at $\sim 10^{-3}M_{\odot}\>{\rm yr}^{-1}$ in a massive core of $\sim$ 240 to 480 $M_{\odot}$ in high $\Sigma_{\rm cl}\gtrsim 1.0\ \rm g\ cm^{-2}$ clump environments. The protostar is slightly inclined $\sim$ 30$\arcdeg$. Walsh et al. (1997) concluded that if the region were powered by a single star, it would have to be an O5.5 star with a luminosity of $3.1\times 10^{5}L_{\odot}$, which agrees well with the isotropic luminosities returned by our models. The viewing angle is close to the outflow half opening angle, resulting in a relatively flat SED shape at shorter wavelengths.

G35.58-0.03: The best models prefer a massive protostar of $\sim$ 24 to 64 $M_{\odot}$ accreting at $\sim 10^{-3}M_{\odot}\>{\rm yr}^{-1}$ in a massive core of $\sim$ 320 to 480 $M_{\odot}$ in high $\Sigma_{\rm cl}\gtrsim 1.0\>\rm g\>cm^{-2}$ clump environments. We also tried fitting the SEDs with the flux measured in variable apertures without setting short wavelength data as upper limits, which exclude the flux from the secondary source to the east at short wavelengths. The best five models have almost the same range for $M_{c}$, $\Sigma_{\rm cl}$, $m_{*}$, $\dot{m}_{*}$ and $L_{\rm bol,iso}$ (there is one model having $m_{*}\sim 96\>M_{\odot}$) as our fiducial case. An early-type star equivalent to an O6.5 star is postulated to have formed within the HC HII region based on the derived Lyman continuum photon number in Zhang et al. (2014). The molecular envelope shows evidence of infall and outflow with an infall rate of $0.05\ M_{\odot}\>\rm yr^{-1}$ and a mass loss rate of $5.2\times 10^{-3}\>M_{\odot}\>\rm yr^{-1}$, which is somewhat higher than our derived disk accretion rate, but may reflect infall on larger scales.

IRAS 16562-3959: There are only 4 fully valid data points constraining the fitting. Since we have 5 free parameters and the $\chi^{2}$ is derived by dividing the number of total data points including those as upper limits, the small number of fully valid data points largely leads to the relatively small $\chi^{2}$. The first four best models tend to give high core masses $\sim$ 320 to 480 $M_{\odot}$ and low $\Sigma_{\rm cl}\lesssim 0.3\>\rm g\>cm^{-2}$ clump environments, while the fifth best model gives a less massive initial core of $120\>M_{\odot}$ and a much denser $\Sigma_{\rm cl}\sim 3.2\>\rm g\>cm^{-2}$ clump. Note in the first three models the core radii are larger than the aperture radius. The bolometric luminosity of the source is reported to be $5-7\times 10^{4}\>L_{\odot}$ by Lopez et al. (2011), which agrees well with most of the models. Guzman (2010) also fit this source with Robitaille et al. (2007) models. The stellar mass of their result 14.7 $M_{\odot}$ is close to our fourth and fifth best models. Their disk accretion rate $5.5\times 10^{-4}M_{\odot}\ \rm yr^{-1}$ is closest to our fourth best model. Their envelope mass 1700 $M_{\odot}$ is much larger than our results. Guzmán et al. (2011) estimated the inclination angle of the SE-NW outflow to be 80$\arcdeg$, which is similar to our second best model.

G305.20+0.21: We have the most data for this source to constrain the model fitting. The initial core mass returned is moderate, ranging from 80 to 160 $M_{\odot}$. Consistently, the envelope mass for this source is also much lower than previous sources. The stellar mass ranges from 16 to 48 $M_{\odot}$, accreting at a high rate $\sim 10^{-3}M_{\odot}\>{\rm yr}^{-1}$. Four models give $\Sigma_{\rm cl}$ as high as $3.2\>\rm g\>cm^{-2}$ and one gives $\Sigma_{\rm cl}\sim 1.0\>\rm g\>cm^{-2}$. The viewing angle is close to the outflow half opening angle, resulting in a flat SED shape at short wavelengths. The extrapolated IRAS luminosity is $\gtrsim 10^{5}\>L_{\odot}$ (Walsh et al. 2001), which is consistent with the $L_{\rm bol}$ derived here.

G49.27-0.34: The models at short wavelengths are much lower than the data points, perhaps indicating significant PAH emission or small dust grain emission from additional heating sources in the region. The best five models all return $m_{*}$ of 12 $M_{\odot}$ and $\Sigma_{\rm cl}$ of $3.2\>\rm g\>cm^{-2}$. The initial core mass ranges from 120 to 320 $M_{\odot}$. The accretion rates are $\sim 10^{-3}\>M_{\odot}\>{\rm yr}^{-1}$.

G339.88-1.26: The best four models prefer a protostar of $\sim 12\>M_{\odot}$ accreting at $\sim 2\times 10^{-4}\>M_{\odot}\>\rm yr^{-1}$ in massive cores of 320 to 480 $M_{\odot}$ in clumps with low $\Sigma_{\rm cl}\sim 0.3\>\rm g\>cm^{-2}$. Alternatively, the fifth best model gives a less massive initial core mass of 120 $M_{\odot}$, but a much denser clump environment with $\Sigma_{\rm cl}\sim 3.2\>\rm g\>cm^{-2}$ and a higher accretion rate of $\sim 10^{-3}\>M_{\odot}\>\rm yr^{-1}$. The bolometric luminosity has been estimated to be 6.4$\times 10^{4}L_{\odot}$ from the SED fitting to infrared fluxes with Robitaille et al. (2007) models in Mottram et al. (2010, 2011), which is similar to the luminosities in our five best models.

Recent ALMA observations (Zhang et al. 2019) reveal collimated CO outflows with a half opening angle of $\sim$ 10$\arcdeg$. In particular, they determine the outflow to be much edge-on so the second model here with $i\approx 20\arcdeg$ is favored. They also estimate the dynamical mass from the gas kinematics as $\sim$ 11 $M_{\odot}$, which is also consistent with our results.

In summary, the massive protostellar sources investigated in this paper tend to have very massive initial cores, high protostellar masses and high accretion rates. The mass surface densities of the clump environments show significant variation. The high envelope masses indicate the protostars are still in an active stage of accretion. Viewing angles tend to be more face-on than edge-on. This allows shorter wavelength photons to more easily escape through the outflow cavities towards the observer, though still regulated partially by extinction of core infall envelope and foreground clump material. Since SOMA survey sources have been selected based on their previously known MIR emission, it is not surprising that the sample may have such a bias towards having more face-on inclinations. Future studies examining inclinations constrained from MIR image intensity profiles and outflow kinematics will allow better measurement of source orientations and a more thorough examination of this effect.

5 Discussion

Compared with the first eight protostars in Paper I, we have extended the upper limit of the luminosity range by one order of magnitude as shown in Figure 16. The seven sources in this paper are more luminous, and thus likely to be more massive protostars embedded in higher mass cores. However, there is the caveat of there being multiple sources sometimes present.

Figure 17 shows the $\chi^{2}$ distribution in $\Sigma_{\rm cl}$ - $M_{c}$ space, $m_{*}$ - $M_{c}$ space and $m_{*}$

• $\Sigma_{\rm cl}$ space for 6 of the sources, i.e., all except G45.12 due to its large $\chi^{2}$. These diagrams illustrate the full constraints in the primary parameter space derived by fitting the SED data, and the possible degeneracies among these parameters. Thus these diagrams give a fuller picture of potential protostellar properties than just the best five models.

Similar to Paper I, $M_{c}$ and $m_{*}$ are relatively well constrained, while $\Sigma_{\rm cl}$ usually spans the full range (for G49.27 the best five models return a universal $\Sigma_{\rm cl}$ of 3.2 g cm*-2* though). The best models ($\chi^{2}-\chi^{2}_{\rm min}<5$, within the red contours) tend to occupy a region with lower $M_{c}$ at higher $\Sigma_{\rm cl}$ and higher $M_{c}$ at lower $\Sigma_{\rm cl}$, similar to the sources in Paper I as discussed in ZT18. The black dashed line denotes a constant $R_{c}$ with $R_{c}=R_{\rm ap}$ using $R_{c}=0.057(\Sigma_{\rm cl}/{\rm g\ cm}^{-2})^{-1/2}\ (M_{c}/60M_{\odot})^{1/2}$ pc (MT03). Parameter sets higher than this line mean they have a $R_{c}$ smaller than $R_{\rm ap}$, which is more physical since we assume the aperture we choose covers the whole envelope. This line only appears in IRAS 16562 because in other sources the $R_{\rm ap}$ is so large that they all appear to the right of the available $\Sigma_{\rm cl}$ - $M_{c}$ space. We can see for most sources at least the best models satisfy this criterion.

In Figure 18 we show the bolometric luminosity spectral energy distributions of the seven high luminosity protostars of this paper, together with the sample from Paper I. Here the $\nu F_{\nu}$ SEDs have been scaled by $4\pi d^{2}$ so that the height of the curves gives an indication of the luminosity of the sources assuming isotropic emission. The ordering of the vertical height of these distributions is largely consistent with the rank ordering of the predicted isotropic luminosity of the protostars from the best-fit ZT models (the legend in Figure 18 lists the sources in order of decreasing ZT best model isotropic luminosity). The curve of G305.20 appears higher than IRAS 16562. However, if we look at all the five best models the isotropic luminosity of G305.20 and IRAS 16562 are actually quite close. The foreground extinction of G305.20 is also generally lower than IRAS 16562, which leads to a higher $4\pi d^{2}\nu F_{\nu}$. Similarly, the foreground extinction of G339.88 is on average lower than G49.27, so that G339.88 has a larger height of the bolometric luminosity SED.

We find no obvious systematic variation in SED shape with varying luminosity. This was investigated by plotting the slope between 19 $\mu$m and 37 $\mu$m versus the isotropic luminosity of the sources (not shown here). We also investigated the relation between $\Sigma_{\rm cl}$, $M_{c}$ and $m_{*}$ in Figure 19. To form high-mass stars naturally requires relatively massive cores (this assumption is built in to the models). However, $\Sigma_{\rm cl}$ does not have to be very high. However, the models with $\Sigma_{\rm cl}\sim 0.1\>\rm g\>cm^{-2}$ have $R_{\rm core}>R_{\rm ap}$ most of the time, which is physically inconsistent with the analysis method. The models with $\Sigma_{\rm cl}\sim 0.3\>\rm g\ cm^{-2}$ only have $R_{\rm core}>R_{\rm ap}$ occasionally, while the other models with higher $\Sigma_{\rm cl}$ do not have such a problem. Thus it is massive protostellar core models with $\Sigma_{\rm cl}\gtrsim 0.3\>\rm g\>cm^{-2}$ surrounding clump environments that are currently consistent with the observed sources.

Overall the ZT models can fit the observed SEDs reasonably well assuming a single protostar forming through an axisymmetric monolithic collapse from a massive core. Only in G45.12, which has stronger evidence for their being multiple protostars that are part of a forming cluster, do the models fare badly and have relatively large values of $\chi^{2}$ (although this may also be due to its extreme luminosity causing it to be near the edge of the ZT model grid). There are reported examples of quite ordered protostellar cores, i.e., with collimated, symmetric outflows: e.g., the case of the early-stage protostar C1-Sa (Tan et al. 2016) and G339.88-1.26 (Zhang et al. 2019, presenting follow-up ALMA observations of one of these SOMA sources). On the other hand, there are also cases that appear much more disordered in both their accretion flows (W51e2e, W51e8, and W51 north, Goddi et al. 2018) and outflows (Orion KL, Bally et al. 2017). The combination of MIR to FIR SED and image fitting with high resolution studies of infall and outflow morphologies for larger samples will allow us to better determine the limitations of simple axisymmetric protostellar core models for Galactic massive star formation studies.

6 Conclusions

We have presented the results of MIR and FIR observations made towards the next seven highest luminosity protostars in the SOMA survey, built their SEDs and fit them with RT models of massive star formation via the Turbulent Core Accretion model. Our goal has been to expand the observational massive protostar sample size to test the star formation models over a wider range of properties and environments and investigate trends and conditions in their formation. Compared with the first eight protostars in Paper I, the seven YSOs in this paper are more luminous, and thus likely to be more massive protostars. Some of the new sources appear to be in more clustered environments and/or have lower-mass companions relatively nearby. In summary, our main results and conclusions are as follows.

1. The MIR emission of massive protostars is strongly influenced by outflow cavities, where extinction is relatively low. We see MIR extension along detected outflows in IRAS16562 and G339.88. Away from these cavities, extinction can be very high and block MIR emission. There is also a hint that the MIR emission may reveal the presence of the optically thick disk perpendicular to the outflow as in G309.92 and G305.20, though more evidence of the position of the protostar from mm or radio continuum observations will be needed to confirm the disk. The high extinction in the MIR tells us that large quantities of high column density material is present close to the protostar, as expected in the Turbulent Core model.

2. The sources span a luminosity range of $10^{4}-10^{6}\>L_{\odot}$. Fitting the SEDs with RT models yields protostellar masses $m_{*}\sim 12-64\>M_{\odot}$ accreting at rates of $\dot{m}_{*}\sim 10^{-4}-10^{-3}\>M_{\odot}\>\rm yr^{-1}$ inside cores of initial masses $M_{c}\sim 100-500\>M_{\odot}$ embedded in clumps with mass surface densities $\Sigma_{\rm cl}\sim 0.1-3\>\rm g\ cm^{-2}$. The relatively high protostellar mass in several sources is possibly due to there being more than one protostar in the region and the $m_{*}$ derived could be the sum of multiple sources.

3. The SED shape, especially the slope at short wavelengths, appears related to the viewing angle and the outflow opening angle. When the viewing angle is close to the outflow opening angle, a relatively flat slope at short wavelengths results. However, the SED shape, especially the location of the SED peak, is also likely to be related to the evolutionary stage of the protostar: more evolved protostars tend to peak at relatively shorter wavelengths. So far we do not see obvious relations between SED shape and bolometric luminosity.

4. To form high-mass stars naturally requires high values of $M_{c}$, but not seem to require especially high values of $\Sigma_{\rm cl}$. We see high-mass protostars are able to at least form from $\Sigma_{\rm cl}\ \gtrsim\ 0.3\ \rm g\ cm^{-2}$ environments.

5. Radiative transfer models based on the Turbulent Core Accretion scenario can reasonably well describe the observed SEDs of most relatively isolated massive protostars, but may not be valid for the most luminous regions in the sample, which may be better treated as protoclusters containing multiple sources. Whether or not core accretion models can apply on smaller physical scales within these regions requires higher angular resolution MIR to FIR observations.

We thank the anonymous referee for reading the pa- per carefully and providing useful comments. We thank Crystal Brogan and Adam Ginsburg for helpful discussion and suggestions. M.L. and J.C.T. acknowledge funding from NASA/USRA/SOFIA. J.C.T. acknowledges support from NSF grant AST1411527 and ERC project 788829 - MSTAR. K.E.I.T. acknowledges support by NAOJ ALMA Scientific Research Grant Numbers 2017-05A. This work is additionally based on observations obtained at the Gemini Observatory [programs GS-2004A-Q-7 and GS-2005A-DD-5], which is operated by the Association of Universities for Research in Astronomy, Inc., under a cooperative agreement with the NSF on behalf of the Gemini partnership: the National Science Foundation (United States), the National Research Council (Canada), CONICYT (Chile), Ministerio de Ciencia, Tecnologa e Innovacin Productiva (Argentina), and Ministrio da Cincia, Tecnologia e Inovao (Brazil).