Young Stellar Objects in the Massive Star-Forming Regions W51 and W43

TL;DR

This study analyzes young stellar objects in W51 and W43, revealing their populations, spatial distributions, relative ages, and associations with massive star formation tracers using infrared data and clustering methods.

Contribution

It provides a detailed catalog of YSOs, identifies MYSO candidates, and examines their spatial and evolutionary relationships in two major star-forming regions.

Findings

Identified over 4000 YSOs across both regions.

Found similar cluster densities despite observational limitations.

Linked MYSO candidates with H II regions in W51.

Abstract

We present the results of our investigation of the star-forming complexes W51 and W43, two of the brightest in the first Galactic quadrant. In order to determine the young stellar object (YSO) populations in W51 and W43 we used color-magnitude relations based on Spitzer mid-infrared and 2MASS/UKIDSS near-infrared data. We identified 302 Class I YSOs and 1178 Class II/transition disk candidates in W51, and 917 Class I YSOs and 5187 Class II/transition disk candidates in W43. We also identified tens of groups of YSOs in both regions using the Minimal Spanning Tree (MST) method. We found similar cluster densities in both regions even though Spitzer was not able to probe the densest part of W43. By using the Class II/I ratios, we traced the relative ages within the regions and based on the morphology of the clusters we argue that several sites of star formation are independent of one another in terms of their ages and physical conditions. We used spectral energy distribution (SED)-fitting to identify the massive YSO (MYSO) candidates since they play a vital role in the star formation process and then examined them to see if they are related to any massive star formation tracers such as UCH II regions, masers and dense fragments. We identified 17 MYSO candidates in W51, and 14 in W43, respectively and found that groups of YSOs hosting MYSO candidates are positionally associated with H II regions in W51, though we do not see any MYSO candidates associated with previously identified massive dense fragments in W43.

Tables

Table 1. Table 1 Astronomical Observation Requests for W51

AORKEY	Date	R.A.	Decl.	Frame	IRAC
	(UT)	(J2000)	(J2000)	Time	Reduction
					Pipeline ver.
12241920	2004 Oct 10	19:21:33	14:20:08	2	S18.25.0
12242176	2004 Oct 10	19:20:49	14:00:17	2	S18.25.0
12243200	2004 Oct 10	19:19:22	13:20:24	2	S18.25.0
12243456	2004 Oct 10	19:20:05	13:40:20	2	S18.25.0
12244480	2004 Oct 10	19:22:17	14:40:03	2	S18.25.0
12244736	2004 Oct 11	19:25:15	15:59:30	2	S18.25.0
10127616	2004 Oct 29	19:21:40	14:04:36	30	S18.25.0
13828864	2005 Oct 21	19:21:40	14:40:36	30	S18.25.0
48640768	2013 Jun 26	19:21:40	14:40:36	30	S19.1.0
51112704	2014 Dec 04	19:21:45	13:49:21	30**HDR data	S19.1.0
51112448	2014 Dec 10	19:19:51	14:01:37	30**HDR data	S19.1.0

Table 2. Table 2 Astronomical Observation Requests for W43

AORKEY	Date	R.A.	Decl.	Frame	IRAC
	(UT)	(J2000)	(J2000)	Time	Reduction
					Pipeline ver.
9237248	2004 Apr 22	18:49:12	-01:06:39	2	S18.25.0
9241344	2004 Apr 22	18:43:55	-03:40:47	2	S18.25.0
9242880	2004 Apr 22	18:45:22	-02:59:36	2	S18.25.0
9244416	2004 Apr 22	18:44:40	-03:18:54	2	S18.25.0
9244672	2004 Apr 22	18:45:55	-02:43:12	2	S18.25.0
9240576	2004 Apr 22	18:47:43	-01:50:48	2	S18.25.0
9246208	2004 Apr 22	18:47:02	-02:10:14	2	S18.25.0
9246464	2004 Apr 22	18:48:28	-01:28:42	2	S18.25.0
21269248	2007 May 12	18:39:17	-01:32:16	2	S18.25.0
21272320	2007 May13	18:53:31	-02:28:02	2	S18.25.0
21326080	2007 May 13	18:54:17	-02:44:47	2	S18.25.0
21274624	2007/ May 13	18:54:13	-03:01:57	2	S18.25.0
45847296	2012 Nov 24	18:47:55	-01:35:37	30**HDR data	S19.1.0
45847552	2012 Nov 24	18:46:23	-02:31:00	30**HDR data	S19.1.0
45847808	2012 Nov 25	18:48:55	-02:23:00	30**HDR data	S19.1.0
45848064	2012 Nov 21	18:56:24	-01:48:00	30**HDR data	S19.1.0
45848320	2012 Nov 25	18:50:09	-01:27:54	30**HDR data	S19.1.0
51107840	2014 Nov 29	18:46:52	-01:51:23	30**HDR data	S19.1.0
51108352	2014 Jun 19	18:47:43	-01:48:24	30**HDR data	S19.1.0
51108608	2014 Nov 29	18:48:04	-02:16:50	30**HDR data	S19.1.0
51108864	2014 Nov 29	18:48:19	-01:46:17	30**HDR data	S19.1.0

Table 3. Table 3 Summary of Source Catalogs

		Number of Sources
Telescope/Instrument	Band	W51	W43
2MASS/UKIDSS	$J$	263,204	512,094
2MASS/UKIDSS	$H$	272,935	592,681
2MASS/UKIDSS	$K_s$	270,380	606,648
Spitzer/IRAC	3.6 $\mu$m	224,559	580,315
Spitzer/IRAC	4.5 $\mu$m	189,793	540,343
Spitzer/IRAC	5.8 $\mu$m	58,550	191,993
Spitzer/IRAC	8.0 $\mu$m	42,372	126,838
Spitzer/MIPS	24.0 $\mu$m	1505	3425

Table 4. Table 4 Source-matching results

			Number of matchesaaMatch percentages are shown in parentheses.
Catalog A	Catalog B	$r_{match}$	W51	W43
IRAC	2MASS	1$″$.0	214, 214 (75%)	370, 426 (58%)
IRAC+2MASS	UKIDSS	1$″$.0	276, 324 (97%)	622, 867 (97%)
IRAC+(2MASS/UKIDSS)	MIPS	1$″$.0	1505 (0.5%)	3425 (0.5%)

Table 5. Table 5 Source Classification Summary

	Number of Objects
Type	W51	W43
Class I	302	917
Class II	1101	4993
Transition Disks	77	194
Embedded protostars	56	144
Faint Class I**These are the YSO candidates with lower confidence.	100	303
Faint Class II**These are the YSO candidates with lower confidence.	1618	4549
Faint Transition Disk**These are the YSO candidates with lower confidence.	35	46
Class III/Photospheres	157, 165	429, 895
AGB candidates	180	866
Other ****Includes PAH emission-dominated sources, H2 shock emission dominated sources, broadline AGN candidates, and PAH galaxy candidates. Unclassified sources lack detection in four bands (either $H{K}_S$, IRAC 1 and 2, or IRAC 1, 2, 3, and 4) or a bright MIPS 24 $\mathrm{\mu m}$ detection.	524	1416
Unclassified	124, 094	198, 409
Total	285, 252	641, 732

Table 6. Table 6 Properties of the Clusters in W51

No.aaSubgroups identified with $d_c$ = 51$″$ are named under the main group with main number and a letter.	R.A. (J2000.0)bbCentral coordinates of the YSO groups or subgroups.	Decl. (J2000.0)bbCentral coordinates of the YSO groups or subgroups.	NIRccNumber of stars with infrared excess. Includes Class I, Class II, deeply embedded protostars, and transition disk candidates.	IddI; Class I, II; Class II candidates	IIddI; Class I, II; Class II candidates	II/IeeNumber in parentheses indicates uncertainty in ratio.	NembffNumber of deeply embedded protostar candidates.	NtdggNumber of transition disk candidates.	NMYSO	Diameter	Diameter
	($°$)	($°$)								($″$)	(pc)
1	290.628161	14.121185	217	53	158	2.98(0.48)	4	2	5	1178.69	30.94
1a	290.607449	14.137096	44	9	35	3.89(1.50)	0	0	1	162.35	4.26
1b	290.671961	14.218522	21	8	12	1.50(0.71)	0	1	1	238.02	6.25
1c	290.595251	14.069474	20	4	16	4.00(2.41)	0	0	0	159.93	4.20
1d	290.693620	14.168545	19	3	15	5.00(3.49)	1	0	1	55.81	1.46
1e	290.539286	14.046652	16	9	7	0.78(0.41)	0	0	1	102.86	2.70
1f	290.655596	13.961203	14	6	8	1.33(0.76)	0	0	0	95.26	2.50
1g	290.629889	14.185851	14	2	12	6.00(5.29)	0	0	0	39.95	1.05
2	290.895659	14.483122	175	55	107	1.94(0.32)	9	4	8	1254.22	32.92
2a	290.911695	14.532360	21	7	13	1.86(0.91)	1	0	2	255.73	6.71
2b	290.968399	14.413716	18	5	10	2.00(1.16)	0	3	0	87.22	2.29
2c	290.844332	14.459009	15	3	10	3.33(2.41)	1	1	0	88.78	2.33
2d	290.962397	14.474087	12	3	9	3.00(2.20)	0	0	0	59.63	1.56
2e	290.798509	14.440177	11	7	4	0.57(0.39)	0	0	0	85.10	2.23
2f	290.762695	14.480864	11	5	5	1.00(0.68)	1	0	4	55.42	1.45
3	291.075362	14.632441	37	8	24	3.00(1.27)	3	2	0	175.98	4.62
3a	291.071001	14.638581	11	2	9	4.50(4.05)	0	0	0	137.68	3.61
4	290.767228	14.327158	35	10	22	2.20(0.86)	3	0	1	111.64	2.93
4a	290.779470	14.342320	12	7	5	0.71(0.45)	0	0	0	214.01	5.62
5	290.863964	14.644746	26	9	12	1.33(0.61)	5	0	1	103.76	2.72
5a	290.868706	14.645758	20	7	9	1.29(0.68)	4	0	0	277.27	7.28
6	290.665391	14.411405	13	1	11	11.00(14.95)	0	1	0	114.24	3.00
7	290.776093	13.867202	12	0	8	$\cdots$	3	1	0	96.51	2.53
8	290.656547	14.319608	11	3	7	2.33(1.77)	1	0	0	216.18	5.67
9	290.743948	13.905521	10	3	6	2.00(1.55)	0	1	1	258.82	6.74

Table 7. Table 7 Properties of the Clusters in W43

No.aaSubgroups identified with $d_c$ = 59$″$ are named under the main group with main number and an asterisk.	R.A. (J2000.0)bbCentral coordinates of the YSO groups or subgroups.	Decl. (J2000.0)bbCentral coordinates of the YSO groups or subgroups.	NIRccNumber of stars with infrared excess. Includes Class I, Class II, deeply embedded protostars, and transition disk candidates.	IddI; Class I, II; Class II candidates	IIddI; Class I, II; Class II candidates	II/IeeNumber in parentheses indicates uncertainty in ratio.	NembffNumber of deeply embedded protostars.	NtdggNumber of transition disk candidates.	NMYSO	Diameter	Diameter
	($°$)	($°$)								($″$)	(pc)
1	281.799575	-2.183919	3047	561	2358	4.20(0.20)	73	55	11	412.90	142.07
1*	281.806380	-1.792551	217	59	143	2.42(0.38)	10	5	1	443.85	11.65
2*	281.912313	-2.081571	140	36	101	2.81(0.55)	2	1	3	554.08	14.54
…	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
2	281.726876	-1.506898	127	20	100	5.00(1.24)	3	4	$\cdots$	209.98	31.76
12*	281.695879	-1.525841	49	12	35	2.92(1.00)	2	0	$\cdots$	247.25	6.49
61*	281.769341	-1.654018	13	5	8	1.60(0.97)	0	0	$\cdots$	160.54	4.21
…	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$
3	281.379116	-3.130973	86	8	72	9.00(3.49)	2	4	$\cdots$	516.63	13.56
….	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$	$\cdots$

Table 8. Table 8 Clusters/Subgroups in W51

Parameter	SLF	$d_c$ = 51$″$
	Method
Number of clusters	9	16
Cutoff distance	82$″$.5(2.17 pc)	51$″$(1.34 pc)
Percent in clusters	60	32
Group diameter	96$″$.5-1178$″$ (2.5-30 pc)	40$″$-256$″$ (1-7 pc)
Class II/I Ratio	2.50(0.25)**Number in parentheses indicates uncertainty in ratio.	2.06(0.27)**Number in parentheses indicates uncertainty in ratio.

Table 9. Table 9 Clusters/Subgroups in W43

Parameter	SLF	$d_c$ = 59$″$
	Method
Number of clusters	51	95
Cutoff distance	81$″$(2.20 pc)	59$″$(1.61 pc)
Percent in clusters	73	43
Group diameter	52$″$-5412$″$ (1.4-142 pc)	29$″$-554$″$ (0.8-14 pc)
Class II/I Ratio	4.58(0.19)**Number in parentheses indicates uncertainty in ratio.	3.71(0.19)**Number in parentheses indicates uncertainty in ratio.

Table 10. Table 10 Physical Parameters of Massive YSO Candidates in W51

No.	Name	R.A.	Decl.	LuminosityaaWeighted mean results from the SED Fitting tool by Robitaille et al. (2007).	$\sigma$bbStandard deviation of the weighted mean values are reported here. ${\chi }^2$ values are used to weight the fits. Minimum 10 fits from the defined best fit criteria was required to calculate the weighted mean values.	MassaaWeighted mean results from the SED Fitting tool by Robitaille et al. (2007).	$\sigma$bbStandard deviation of the weighted mean values are reported here. ${\chi }^2$ values are used to weight the fits. Minimum 10 fits from the defined best fit criteria was required to calculate the weighted mean values.	IRAC	SED	SED
		J2000.0($°$)	J2000.0($°$)	(${10}^3$L⊙)	(${10}^3$L⊙)	(M⊙)	(M⊙)	Classcc0: deeply embedded, 1: Class I, 2: Class II, F: Flat spectrum, 3: transition disks.	Classcc0: deeply embedded, 1: Class I, 2: Class II, F: Flat spectrum, 3: transition disks.	Slopedd Values in parentheses signify error in last two digits of the SED slope value.
$1^e$	SSTOERC$\_$G049.4183+0.3211	290.252258	14.780608	1.14	0.43	6.84	0.94	3	1	1.74(0.06)
$2^e$	SSTOERC$\_$G049.4821-0.3507	290.895782	14.520227	4.12	1.58	8.65	1.27	1	1	2.55(0.2)
$3^e$	SSTOERC$\_$G049.3793-0.256	290.759369	14.474311	2.97	0.99	10.49	1.22	0	1	2.02(0.36)
$4^e$	SSTOERC$\_$G049.3884-0.2634	290.770569	14.478874	1.08	0.85	8.79	0.70	2	F	0.29(0.9)
5	SSTOERC$\_$G048.8743-0.5188	290.752411	13.905082	4.06	1.50	8.72	1.29	1	1	2.12(0.12)
6	SSTOERC$\_$G048.9022-0.2734	290.54303	14.045395	2.22	0.63	7.47	0.73	1	1	2.35(0.12)
7	SSTOERC$\_$G048.9916-0.3143	290.623596	14.104997	5.15	0.12	9.44	0.09	2	1	1.02(0.18)
8	SSTOERC$\_$G049.0693-0.3164	290.663391	14.172537	1.09	0.06	6.28	0.09	2	F	0.11(0.21)
9	SSTOERC$\_$G049.1138-0.2963	290.666656	14.221247	1.15	0.0	6.37	0.00	2	1	0.35(0.08)
10	SSTOERC$\_$G049.0327-0.3431	290.66983	14.127697	5.15	0.0	9.45	0.00	1	1	0.85(0.17)
11	SSTOERC$\_$G049.2213-0.3354	290.754608	14.297651	1.17	0.13	6.39	0.21	2	1	2.02(0.27)
12	SSTOERC$\_$G049.3820-0.2563	290.760986	14.476582	1.16	0.21	6.36	0.33	2	1	0.95(0.17)
13	SSTOERC$\_$G049.3800-0.2696	290.772064	14.468561	5.39	2.36	8.7	1.47	1	1	2.81(0.28)
14	SSTOERC$\_$G049.5688-0.2725	290.866852	14.633584	1.13	0.15	6.33	0.16	2	1	0.45(0.24)
15	SSTOERC$\_$G049.4252-0.354	290.870911	14.468496	1.61	0.59	6.85	0.66	1	1	2.12(0.18)
16	SSTOERC$\_$G049.4491-0.3611	290.889069	14.486297	1.20	0.31	6.43	0.30	2	1	0.63(0.24)
17	SSTOERC$\_$G049.4885-0.3628	290.909912	14.520178	9.47	2.45	11.38	1.1	1	1	2.39(0.13)

Table 11. Table 11 Physical Parameters of Massive YSO Candidates in W43

No.	Name	R.A.	Decl.	LuminosityaaWeighted mean results from the SED Fitting tool by Robitaille et al. (2007).	$\sigma$bbStandard deviation of the weighted mean values are reported here. ${\chi }^2$ values are used to weight the fits. A minimum of 10 fits from the defined best-fit criteria was required to calculate the weighted mean values.	MassaaWeighted mean results from the SED Fitting tool by Robitaille et al. (2007).	$\sigma$bbStandard deviation of the weighted mean values are reported here. ${\chi }^2$ values are used to weight the fits. A minimum of 10 fits from the defined best-fit criteria was required to calculate the weighted mean values.	IRAC	SED	SED
		J2000.0($°$)	J2000.0($°$)	(${10}^3$L⊙)	(${10}^3$L⊙)	(M⊙)	(M⊙)	Classcc0: deeply embedded, 1: Class I, 2: Class II, F: Flat spectrum, 3: transition disks.	Classcc0: deeply embedded, 1: Class I, 2: Class II, F: Flat spectrum, 3: transition disks.	Slopedd Values in parentheses signify error in last two digits of the SED slope value.
$1^e$	SSTOERC$\_$G029.8031-0.0563	281.481842	-2.810057	4.58	3.62	8.72	2.92	1	1	2.42(1.04)
$2^e$	SSTOERC$\_$G030.4203-0.2284	281.917053	-2.339325	1.52	0.51	8.24	0.46	1	1	2.92(1.14)
$3^e$	SSTOERC$\_$G030.9632+0.0828	281.887787	-1.714331	2.24	0.88	7.29	0.75	1	1	2.67(0.19)
4	SSTOERC$\_$G030.9238-0.1319	282.060944	-1.847242	1.43	1.45	6.48	1.24	1	1	1.25(0.45)
5	SSTOERC$\_$G030.5565-0.0894	281.855469	-2.154748	1.15	0.03	6.37	0.06	1	F	0.28(0.56)
6	SSTOERC$\_$G030.7019-0.0338	281.872314	-2.000047	1.44	0.94	6.58	0.99	1	1	0.64(0.22)
7	SSTOERC$\_$G030.6487-0.2003	281.996246	-2.123312	1.02	0.35	6.17	0.85	1	1	2.10(0.31)
8	SSTOERC$\_$G029.2101+0.1301	281.044617	-3.25246	1.46	1.12	6.52	1.24	1	1	1.48(0.34)
9	SSTOERC$\_$G030.0299+0.308	281.260956	-2.441933	1.04	0.02	6.21	0.05	1	1	0.37(0.13)
10	SSTOERC$\_$G031.4202-0.5835	282.689545	-1.61143	1.20	1.25	5.37	2.3	2	1	1.06(0.36)
11	SSTOERC$\_$G029.6445-0.1788	281.518433	-3.00707	1.09	0.25	6.24	0.48	2	F	0.23(0.28)
12	SSTOERC$\_$G030.6549-0.1449	281.949799	-2.092496	2.08	2.64	6.23	2.81	2	1	1.25(0.39)
13	SSTOERC$\_$G030.5984-0.1106	281.893402	-2.127159	1.61	0.05	6.95	0.06	2	1	1.71(0.3)
14	SSTOERC$\_$G029.8775-0.0379	281.49942	-2.735451	1.15	0.01	6.37	0.02	2	F	-0.11(0.34)

Table 12. Table 12 Clustering Analysis Summary

Region	Number of Clusters	Cutoff	Clustered	Class II/I	Method
		distance ($″$)	Fraction	Ratio
W51	9	82$″$.5(2.17 pc)	60	2.50(0.25)	aaStraight-line fit method
	16	51$″$(1.34 pc)	32	2.06(0.27)	bbThe method using the short branch length distance enables us to see subgroups in the region.
W43	51	81$″$(2.20 pc)	73	4.58(0.19)	aaStraight-line fit method
	95	59$″$(1.61 pc)	43	3.71(0.19)	bbThe method using the short branch length distance enables us to see subgroups in the region.
G305ccSaral et al. (2015)	15	48$″$(0.90 pc)	76	3.50(0.30)	aaStraight-line fit method
	32	30$″$(0.60 pc)	44	3.00(0.30)	bbThe method using the short branch length distance enables us to see subgroups in the region.
G333ccSaral et al. (2015)	11	73$″$(1.30 pc)	81	4.50(0.40)	aaStraight-line fit method
	45	52$″$(0.90 pc)	54	4.00(0.40)	bbThe method using the short branch length distance enables us to see subgroups in the region.
W49ccSaral et al. (2015)	3	96$″$(5.20 pc)	52	2.11(0.43)	aaStraight-line fit method
	7	40$″$(2.20 pc)	21	2.00(0.64)	bbThe method using the short branch length distance enables us to see subgroups in the region.
Rescaled W51	5	120$″$.94(6.47)	74	1.61(0.21)	aaStraight-line fit method
Rescaled W43	8	91$″$.82(4.91 pc)	73	2.79(0.21)	aaStraight-line fit method
Rescaled G305ccSaral et al. (2015)	4	52$″$(2.78 pc)	73	0.55(0.15)	aaStraight-line fit method
Rescaled G333ccSaral et al. (2015)	2	78$″$(4.22 pc)	70	0.03(0.04)	aaStraight-line fit method

Table 13. Table 13 Source List for W51

Source NameaaSources are named according to their Galactic longitude and latitude with the prefix SSTOERC, referring to Spitzer Space Telescope, Origin and Evolution of Rich Clusters project.	R.A.	Decl.	$J$	$H$	$K_s$	$\left[3.6\right]$	$\left[4.5\right]$	$\left[5.8\right]$	$\left[8.0\right]$	$\left[24\right]$			SED		R-
	J2000.0	J2000.0	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	Typebb0; deeply embedded, 1; Class I, 2; Class II, 3; transition disks, 4; faint Class I, 5; faint Class II, 6; faint transition disks, 99; Class III and photospheres, 9; shocked gas emission, 29; AGNs, 19; PAH galaxies, 20; PAH-dominated sources, 12; AGBs, -100; unclassified.	MethodccSources are flagged according to the method by which they were classified.1; Gutermuth et al. (2009) classification, 2; AGBs, faint YSO elimination method (see Section 2.3).	slope $\alpha$	Sigmadd Sigma values signify the error in last two digits of the SED slope value.	Squared
	($°$)	($°$)
SSTOERC$\_$G049.3793$-$0.2560	290.75937	14.47431	17.58(03)	15.17(02)	13.79(02)	11.76(20)	11.28(17)	$\cdots$	$\cdots$	0.62(17)	0	1	2.02	0.36	0.99
SSTOERC$\_$G048.8881$-$0.2638	290.52740	14.03759	17.28(03)	15.13(02)	13.90(02)	11.74(02)	11.09(02)	6.11.26(18)	$\cdots$	$\cdots$	1	1	0.81	0.29	0.84
SSTOERC$\_$G048.6453$-$0.1257	290.28412	13.88838	15.61(02)	14.09(02)	13.30(02)	12.80(03)	12.44(04)	$\cdots$	$\cdots$	$\cdots$	2	1	-1.75	0.06	1.00
SSTOERC$\_$G048.6287+0.2361	289.94690	14.04359	14.08(02)	12.27(03)	11.38(02)	10.79(01)	10.60(01)	10.47(07)	9.87(18)	4.50(16)	3	1	-1.71	0.30	0.90
SSTOERC$\_$G049.1638$-$0.0498	290.46670	14.38155	$\cdots$	18.25(10)	16.99(08)	13.70(12)	13.16(09)	$\cdots$	$\cdots$	$\cdots$	4	2	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G049.0780+0.5502	289.87799	14.58771	17.23(03)	16.03(03)	14.75(03)	14.01(09)	13.70(14)	$\cdots$	$\cdots$	$\cdots$	5	2	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G049.1458$-$0.0561	290.46362	14.36270	14.66(02)	13.97(02)	13.69(02)	13.29(01)	12.97(14)	$\cdots$	$\cdots$	2.11(10)	6	2	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G048.9561+0.4943	289.87006	14.45382	17.38(03)	15.56(02)	14.74(02)	13.86(10)	13.92(14)	$\cdots$	$\cdots$	$\cdots$	99	1	-1.48	0.33	0.99
SSTOERC$\_$G048.8942$-$0.4828	290.72949	13.93967	14.65(02)	13.68(02)	13.19(02)	12.40(03)	12.37(04)	11.70(09)	10.29(07)	$\cdots$	19	1	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G049.3857$-$0.0146	290.5426	14.5938	17.89(04)	16.65(03)	15.53(03)	13.30(04)	12.18(04)	11.77(09)	10.99(15)	$\cdots$	9	1	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G049.1445+0.2417	290.19168	14.5016	13.43(03)	10.24(03)	8.41(03)	6.53(00)	5.82(00)	4.68(00)	4.08(00)	$\cdots$	20	1	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G049.3055+0.6570	289.89069	14.83866	10.18(03)	8.48(05)	7.40(03)	6.68(01)	6.18(00)	5.59(00)	5.30(00)	$\cdots$	12	2	$\cdots$	$\cdots$	$\cdots$

Table 14. Table 14 Source List for W43

Source NameaaSources are named according to their Galactic longitude and latitude with the prefix SSTOERC, referring to Spitzer Space Telescope, Origin and Evolution of Rich Clusters project.	R.A.	Decl.	$J$	$H$	$K_s$	$\left[3.6\right]$	$\left[4.5\right]$	$\left[5.8\right]$	$\left[8.0\right]$	$\left[24\right]$			SED		R-
	J2000.0	J2000.0	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	(mag)	Typebb0; deeply embedded, 1; Class I, 2; Class II, 3; transition disks, 4; faint Class I, 5; faint Class II, 6; faint transition disks, 99; Class III and photospheres, 9; shocked gas emission, 29; AGNs, 19; PAH galaxies, 20; PAH-dominated sources, 12; AGBs, -100; unclassified.	MethodccSources are flagged according to the method by which they were classified.1; Gutermuth et al. (2009) classification, 2; AGBs, faint YSO elimination method (see Section 2.3).	slope $\alpha$	Sigmadd Sigma values signify the error in last two digits of the SED slope value.	Squared
	($°$)	($°$)
SSTOERC$\_$G031.2284$-$0.0191	282.09952	-1.52473	19.87(25)	18.72(21)	17.84(23)	$\cdots$	12.25(07)	10.63(11)	9.60(14)	4.58(08)	0	1	1.41	0.32	0.90
SSTOERC$\_$G030.4270+0.4739	281.29474	-2.01290	16.28(11)	18.09(09)	15.07(03)	12.87(08)	12.48(09)	$\cdots$	$\cdots$	5.89(24)	1	1	0.72	0.23	0.98
SSTOERC$\_$G029.5402+0.4199	280.93744	-2.82646	17.41(03)	15.09(02)	13.92(02)	12.60(04)	11.94(03)	11.29(07)	10.36(08)	7.58(12)	2	1	-0.35	0.04	0.98
SSTOERC$\_$G030.3208+0.9046	280.86285	-1.91069	$\cdots$	$\cdots$	$\cdots$	10.35(01)	10.29(01)	10.25(03)	10.01(05)	7.49(14)	3	1	-2.41	0.27	1.00
SSTOERC$\_$G030.3622$-$0.2990	281.95340	-2.42323	$\cdots$	18.28(08)	16.53(05)	13.81(09)	12.78(03)	$\cdots$	$\cdots$	6.28(28)	4	2	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G029.6605+0.1712	281.21396	-2.83302	17.38(03)	15.81(02)	14.66(02)	13.35(08)	12.81(07)	$\cdots$	$\cdots$	7.75(19)	5	2	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G029.4575$-$0.1408	281.39908	-3.15610	14.97(02)	14.00(02)	13.56(02)	13.20(09)	13.08(09)	$\cdots$	$\cdots$	7.59(12)	6	2	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G030.5678+0.8209	281.05017	-1.72920	12.72(03)	10.72(03)	9.89(02)	9.32(01)	9.45(01)	9.21(02)	9.17(03)	9.63(03)	99	1	-2.62	0.28	0.99
SSTOERC$\_$G030.6965+0.7979	281.12946	-1.62522	14.32(02)	13.39(02)	12.96(02)	12.39(03)	12.40(03)	11.68(10)	10.19(07)	7.56(08)	19	1	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G030.7902+0.2053	281.69977	-1.81238	$\cdots$	18.53(15)	15.41(03)	11.39(02)	10.05(02)	9.85(05)	9.61(09)	$\cdots$	9	1	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G030.2314+0.8732	280.84994	-2.00452	$\cdots$	$\cdots$	$\cdots$	8.89(01)	9.21(01)	8.84(01)	8.77(02)	9.33(02)	20	1	$\cdots$	$\cdots$	$\cdots$
SSTOERC$\_$G030.7677+0.4362	281.48395	-1.72702	$\cdots$	13.36(05)	10.94(04)	8.85(01)	8.58(01)	8.25(01)	7.82(01)	8.40(27)	12	1	$\cdots$	$\cdots$	$\cdots$

Equations

$$[4.5]-[5.8]>0.7\text{ and}[3.6]-[4.5]>0.7$$

$$[4.5]-[8.0]-{\sigma }_3>0.5$$

$$[3.6]-[5.8]-{\sigma }_4>0.35$$

$$[3.6]-[5.8]+{\sigma }_4\le \frac{0.14}{0.04}\times (([4.5]-[8.0]-{\sigma }_3)-0.5)+0.5$$

$$[3.6]-[4.5]-{\sigma }_4>0.15$$

$$[24]<7\text{ and}[X]-[24.0]>4.5$$

$$[5.8]-[24]>2.5\text{ or}[4.5]-[24]>2.5$$

$$3.5<[3.6]<9.5\text{ and}0.4<[3.6]-[8.0]<2.6$$

$$3<[3.6]<9.5\text{ and}0.2<[3.6]-[4.5]<1.25$$

$$\alpha \ge 0.3\text{ Class I}$$

$$0.3>\alpha \ge -0.3\text{ Flat Spectrum}$$

$$-0.3>\alpha \ge -1.6\text{ Class II}$$

$$-1.6>\alpha \ge -2.7\text{ Class III}$$

Full text

††slugcomment: 20 April 2017; Published in ApJ

1. Introduction

Near-infrared and mid-infrared observations have revealed a broad range of stellar densities of YSOs in star-forming regions (Bressert et al., 2010). While we have a detailed picture of the formation of low-mass stars with advances in observational techniques and improved theoretical knowledge in the last decade, our understanding of massive star formation is still lagging. Massive stars form and evolve quickly and they produce strong winds and outflows which heat, ionize, and disrupt their natal molecular cloud, which makes the observations of early evolutionary stages difficult. In addition, massive stars are forming in a clustered mode by interacting with each other as seen in Orion (Lada, 1992; Evans et al., 2009; Megeath et al., 2016) or Cygnus-X (Motte et al., 2007; Beerer et al., 2010; Kryukova et al., 2014; Schneider et al., 2016). In order to understand the star formation process, it is necessary to understand that the formation of massive clusters and the Galactic star-forming regions are the ideal places for this purpose, since they host many active sites of star formation in different environments at various stages of evolution.

We started a detailed study of a sample of mini-starburst complexes (Motte et al., 2003) to better understand how massive stars form and affect their environments and the effects that they have on the formation of lower-mass stars (Saral & Hora, 2013). Precursor clouds of young massive clusters (YMC) (Bressert et al., 2012; Ginsburg et al., 2012; Walker et al., 2016) are thought to be the potential descendant of mini-starburst clouds/clumps (Motte et al., 2003; Nguyen-Luong et al., 2011a, b; Louvet et al., 2014; Nguyen-Luong et al., 2016). YMCs (M $\gtrsim$ $10^{4}$ M⊙, ages $\lesssim$ 100 Myr) are stellar systems that host a very large number of YSOs with a range of masses including massive ones (Longmore et al., 2014). The first paper of this study (Saral et al., 2015) focused on W49, which is considered to be a mini-starburst region hosting a YMC precursor. We identified thousands of YSOs and tens of massive YSO (MYSO) candidates in W49, which shows a compact massive cluster in its center.

In this second paper, we investigate two YMC precursors, W51 (Walker et al., 2016) and W43, which are known to host two mini-starburst ridges (Nguyen-Luong et al., 2013). They are both among the brightest and most massive of the first Galactic quadrant. We present an imaging and photometric analysis of W51 and W43 with deep IR data from 1 to 24 $\mu$m in order to investigate the MYSOs and the embedded clusters forming around them.

W51 was first detected as an H II region by Westerhout (1958) and later identified as a molecular cloud by CO emission (Penzias et al., 1971). It hosts the star forming regions W51A and W51B and a supernova remnant W51C, with a total mass of $\sim 1.2\times 10^{6}$ (Carpenter & Sanders, 1998) and a stellar mass of $\sim 10^{4}$ (Kumar et al., 2004). Xu et al. (2009) used methanol masers to measure the trigonometric parallaxes and found a distance of $\displaystyle{5.1^{+2.9}_{-1.4}}$ kpc. Sato et al. (2010) reported the distance as $\displaystyle{5.4^{+0.31}_{-0.28}}$ kpc with higher accuracy by using $\textrm{H}_{2}\textrm{O}$ masers and the trigonometric parallax method. There are studies based on maser emissions and velocity measurements showing the interaction between the star-forming region W51B and the supernova remnant W51C (i.e., Green et al., 1997; Ginsburg et al., 2015). Kang et al. (2009) studied the YSO population of W51 by using Spitzer and 2MASS data and combining both the color selection criteria from Simon et al. (2007) and the spectral energy distribution (SED) slope classification (Lada, 1991; André et al., 1993; Greene et al., 1994), in addition to SED models from Robitaille et al. (2006). This resulted in the identification of 737 YSO candidates. Although there is a slight difference between the cataloged region in Kang et al. (2009) study and in this paper, they have 91% of the sources in common; and among them, 54% of the sources are YSO candidates in both studies. We compare our results with this earlier work in Section 2.3.1.

W43 is called a mini-starburst (Motte et al., 2003) with several tens of cores forming high-mass stars, among which are one UCH II region, a handful of maser sources, and at least three very massive Class 0s with $\sim 10^{4}$ ${\rm L}_{\sun}$ in the massive cloud. The starburst cluster was first detected as an unresolved source by Lester et al. (1985) in near-IR images and then confirmed as a massive cluster (W43-Main) by Blum et al. (1997). The complex hosts two bright clouds: W43-Main (G30.8+0.02), and W43-South (G29.96-0.02). The detailed analysis of molecular and atomic gas tracers by Nguyen-Luong et al. (2011a) showed that these two regions are connected and the complex extends over 140 pc at $l$ = (29*∘$-$32∘), $b$ = (-1∘* $-$ +1*∘*). The two clouds are surrounded by an atomic gas envelope with a diameter of 290 pc. The connection between W43-Main and W43-South has been confirmed by Carlhoff et al. (2013), and Motte et al. (2014) confirmed the presence of envelope, giving its diameter as 270 pc and a mass of $M_{\mathrm{HI,env}}\sim 3\times 10^{6}$ M⊙. On the other hand, recently, Bialy et al. (2017) argued that the large HI gas can be explained by several HI to H2 transition layers along the sightlines instead of one large HI envelope based on 21 cm observations and a theoretical model for the HI-to-H2 transition.Nguyen-Luong et al. (2011a) calculated its kinematic distance as 6 kpc and placed it at the meeting point of the Galactic bar and the Scutum-Centaurus Arm, while Zhang et al. (2014) redetermined the distance as $\displaystyle{5.49^{+0.39}_{-0.34}}$ kpc using the trigonometric parallax method by studying the masers G029.86-00.04, G029.95-00.01, G031.28+00.06 and G031.58+00.07.

In Section 2, we describe the observations, data reduction techniques, and our near- and mid-IR source catalog and YSO classification. In Section 3, we present the clustering analysis. Section 4 contains the SED fitting results for massive YSO candidates. In Section 5 we present the massive star formation tracers, we discuss the star formation history and compare W51 and W43 to other star-forming regions. Finally, in Section 6, we summarize our results.

2. Observations and Methods

2.1. IRAC Imaging

We used the mid-infrared Spitzer data from several programs obtained with the Spitzer IRAC instrument (Fazio et al., 2004) at $3.6$, $4.5$, $5.8$, and $8.0$ $\mu$m. The 2 s frame time data taken at all wavelengths in 2007 and prior years are from the GLIMPSE Legacy Survey (Churchwell et al., 2009, program ID 187). The data obtained in 2012 and afterward were performed during the warm mission, so only $3.6$ and $4.5$ $\mu$m data were obtained. We also used 30 s frame time data from other programs taken in October 2004, 2005, and June 2013, with the following program IDs: 2313 (PI: M. Kuchner), 20026 (PI: T. von Hippel), and 90095 (PI: K. Luhman). We also obtained 30 s data in High Dynamic Range (HDR) mode in December 2014 (project ID: 10012, PI: J. Hora), which performs consecutive individual observations with frame times of 1.2 s and 30 s. An RGB image of the W51 region constructed from the 2 s data is shown in Figure 1.

The IRAC data for the W43 region are from the GLIMPSE Legacy Survey (Benjamin et al., 2003; Churchwell et al., 2009, program IDs 186, 30570) and 30 s HDR data were taken with program IDs 10012 and 80058 (PI: J. Hora). An RGB image of the W43 is shown in Figure 2. The large image was constructed from the 2 s frames; the smaller images of W43-Main and W43-South were constructed from the 30 s HDR frames for the $3.6$ and $4.5$ $\mu$m wavelengths and the 2 s images at 8.0 $\mu$m. We list the dates and coordinates of each Astronomical Observation Request (AOR) in Tables 1 and  2.

We generated the mosaic images from the standard basic calibrated data (BCD) products (processed with the Spitzer IRAC pipeline version S18.25.0 and S19.1.0). Automated source detection and aperture photometry were carried out using PhotVis 1.10 (Gutermuth et al., 2004, 2008). PhotVis utilizes a modified DAOphot (Stetson, 1987) source-finding algorithm. Aperture photometry was performed with an aperture of 2.4$\arcsec$ radius and using a background annulus of inner and outer radii 2.4$\arcsec$, and 7.2$\arcsec$ respectively. The IRAC PSF is 1.66$\arcsec$, 1.72$\arcsec$, 1.88$\arcsec$, and 1.98$\arcsec$ at $3.6$, $4.5$, $5.8$ and $8.0$ $\mu$m, respectively (Fazio et al., 2004), corresponding to 0.044$-$0.052 pc for W51 (5.4 kpc) and 0.045$-$0.054 pc for W43 (5.6 kpc). In the first paper we had estimated the completeness magnitudes in the W49 region for $3.6$, $4.5$, $5.8$, and $8.0$ $\mu$m as 15.1, 14.7, 12.35, and 12.12, respectively, and since the exposure times are same and the distances are similar, we assume the completeness estimate will be similar for W51 and W43.

2.2. Source Catalog

In the W43 region, within an area of size $\triangle$$\alpha$ $\times$ $\triangle$$\delta$ = 2$\fdg$0 $\times$ 2$\fdg$3, centered at ($\alpha$, $\delta$) = (281$\fdg$8, -2$\fdg$2), 641,732 sources were detected in the IRAC images. Among these, 118,365 sources have photometry in all four IRAC filters. In W51, within an area of size $\triangle$$\alpha$ $\times$ $\triangle$$\delta$ = 1$\fdg$5 $\times$ 1$\fdg$5, centered at ($\alpha$, $\delta$) = (290$\fdg$6, 14$\fdg$2), 285,252 sources were detected. Among these, 37,932 sources have photometry in all four IRAC filters.

After performing the photometry in each of the four IRAC bands, we combined the catalogs of short (2 s) and long (30 s) frames using a band-merging process that takes all of the sources from the short frames that are brighter than a certain cutoff (magnitude 11, 10.4, 8.7, and 8.2 for $3.6$, $4.5$, $5.8$, and $8.0$ $\mu$m, respectively) based on the saturation limit of 30 s frames. All of the sources detected by Spitzer were matched to Two Micron All Sky Survey (2MASS) Point Source Catalog (Skrutskie et al., 2006) sources by the PhotVis photometry routine. We also cross-matched the IRAC Catalog with the UKIRT Infrared Deep Sky Survey (UKIDSS)555UKIDSS uses the UKIRT Wide Field Camera (Casali et al., 2007) on the United Kingdom Infrared Telescope, and the UKIDSS project is defined by Lawrence et al. (2007) Galactic Plane Survey (Lucas et al., 2008) data from Data Release 8 (DR8PLUS), which are deeper and have better spatial resolution than those of 2MASS. We used a maximum of 1$\arcsec$ for the matching radius to allow a reasonable error in source positions and also to prevent false matches. We took into account the saturation limit in the UKIDSS survey which is near 12.65, 12.5, and 12 mag in $J$, $H$, and $K_{s}$, respectively, and used 2MASS photometry when the source was saturated in UKIDSS (20% of the data in UKIDSS $JHK$ photometry was saturated in W51 and 25% in W43). The UKIDSS data have small but measurable zero-point photometric offsets from 2MASS, therefore we calculated the mean and standard deviations of the magnitude residuals between 2MASS and UKIDSS in both W43 and W51 catalogs and applied a mean offset to UKIDSS data of 0.06, $-$0.06, and $-$0.01 mag to the $J$, $H$, and $K_{s}$ bands, respectively. In addition, due to very small photometric errors in the UKIDSS catalog (e.g., $<$0.001 mag for 13 mag sources), we adjusted the values by adding 0.02 mag in quadrature, which imposes an error floor of 0.02 mag but does not affect the larger errors. Lastly, we matched our catalog with the MIPSGAL Archive 24 $\mu$m data (Gutermuth and Heyer, 2015). We used 1$\arcsec$ as a maximum radial tolerance based on Monte Carlo simulations that give a mean false match probability of 6 % and 5.5 % for W51, and W43, respectively.

In the case of multiple matches, we selected the closest IRAC object to the MIPS source. In W51, only 15 MIPS sources had a second IRAC object within the 1$\arcsec$ search radius that could possibly contribute to its 24 $\mu$m flux. However, these sources were either unclassified or classified as photospheres and only one source was classified as a YSO candidate, while the second IRAC object was an unclassified object. Similarly, in W43, 24 sources had a second IRAC object within the 1$\arcsec$ search radius. Within these matches, in only one pair were two IRAC sources were identified as YSO candidates that might be contributing to each other’s 24 $\mu$m flux. However, in the rest of the 11 pairs, we see that the second IRAC source is unclassified. As a summary, there are only two YSO candidates in W43 that might be contributing to each other’s 24 $\mu$m flux; we consider this as having a minimal effect on our analysis.

A summary of the source catalog and source-matching results can be found in Table 3 and in Table 4, respectively. In the end, we generated a final catalog that contains photometry of sources over a wavelength range from 1.2 to 24 $\mu$m and is a combination of 2MASS/UKIDSS, Spitzer/IRAC, and MIPS data. The final catalogs for W43 and W51 are presented in Table 13 and 14, respectively.

2.3. YSO Classification

2.3.1 IRAC-color Criteria

In this study we used the selection method based on color and magnitude criteria that we applied in Saral et al. (2015) for W49 and defined in Gutermuth et al. (2008, 2009) (hereafter IRAC-color method). The Gutermuth et al. (2009) IRAC-color method is updated to account for complications found in the more active and typically more distant regions. In Phase 1, by using all four IRAC bands we eliminated extragalactic sources such as star-forming galaxies and polycylic aromatic hydrocarbon (PAH) rich galaxies, since they can be classified from their very red $5.8$ and $8.0$ $\mu$m colors (Stern et al., 2005) due to strong PAH emission. In addition, we checked possible broadline active galactic nuclei (AGNs) candidates that have consistent mid-IR colors with YSOs (Stern et al., 2005) by using the [4.5]$-$[8.0] vs. [4.5] color-magnitude diagram as described in Gutermuth et al. (2009). In addition, we removed sources classified as knots of shock emission and PAH emission-dominated sources (Gutermuth et al., 2009). As a result, we eliminated one AGN candidate, 74 PAH-rich galaxies, five knots of shocked gas emission, and 444 PAH-contaminated apertures in the W51 region and 86 PAH-rich galaxies, 13 knots of shocked gas emission, and 1317 PAH-contaminated apertures in the W43 region (Figure 3). No sources were classified as AGN candidates toward the W43 region. After eliminating contaminant sources, by using a combination of different color-magnitude conditions of the four IRAC bands described by Gutermuth et al. (2009) we classified the YSOs into two categories: Class I sources (protostars with circumstellar disks and infalling envelopes), and Class II sources (pre-main-sequence stars with optically thick disks). Class I sources (shown with red color dots in Figure 4 and Figure 5) are classified by their red colors using the following criteria (Gutermuth et al., 2009):

[TABLE]

Class II candidates were classified using the following constraints as described in Gutermuth et al. (2009), while $\sigma$3 and $\sigma$4 are the errors of the colors [4.5]$-$[8.0] and [3.6]$-$[5.8], respectively:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

With this method we can differentiate YSOs from field stars/photospheres (Indebetouw et al., 2005; Flaherty et al., 2007). However, line-of-sight extinction and intrinsic variability can cause a reddened Class II source be misclassified as a Class I source. Therefore, the next step (Phase 2) in the Gutermuth et al. (2009) method is to estimate the extinction ($A_{K}$) simply by using the measured $JHK$ magnitude and, when $J$ is not available, $3.6$ and $4.5$  $\mu$m data. The extinction was measured by using the baseline colors of Classical T Tauri Stars (CTTS) locus from Meyer et al. (1997), standard dwarf colors from Bessell & Brett (1988), and baseline colors of YSO loci from Gutermuth (2005). Following that the measured 1-8  $\mu$m magnitudes were dereddened according to the reddening law by Flaherty et al. (2007). After separating the reddened Class II candidates from the Class I candidates, the third step (Phase 3) was to classify “deeply embedded sources,” which are Class I candidates with bright emission at 24 $\mu$m, and “transition disk candidates,” which are Class II sources with significant dust clearing within their disks which generates bright emission at 24 $\mu$m. Embedded protostar candidates are classified if they meet the following criteria:

[TABLE]

where [X] is the longest wavelength detection with IRAC. Finally, transition disk candidates are classified by the following criteria:

[TABLE]

Photospheres, transition disk candidatates, and embedded protostars can be seen in Figure 6.

Following this method we also took into account the greater distance, and therefore the possible higher level of contamination, from objects such as AGB stars, background sources, and other extragalactic contaminants that have probably been classified as YSOs. Since AGB stars are bright and, in general, slightly bluer than YSOs (Gutermuth et al., 2008; Koenig & Leisawitz, 2014), we reclassified the initially classified bright YSOs that also follow the following selection criteria as candidate AGB stars:

[TABLE]

or

[TABLE]

We classified 180 sources in the W51 region and 866 sources in the W43 region as candidate AGB stars, with their bright magnitudes shown in $[3.6]-[4.5]$ versus $[3.6]$ and $[3.6]-[8.0]$ versus $[3.6]$ color$-$magnitude diagrams (shown in the bottom panels of Figure 4 for W51 and Figure 5 for W43). The Gutermuth et al. (2009) criteria require detections at $5.8$ or $8.0$ $\mu$m to identify extragalactic and background contaminants. However, we have only 37,932 sources in W51 out of 285,252 and 118,365 sources in W43 out of 641,732 that are detected in all four IRAC bands. It is therefore possible to have contaminant sources in our list of YSO candidates that have been classified with only IRAC $3.6$ and $4.5$ $\mu$m photometry. Thus, we applied a selection cut of [$3.6$] ¿ 13 mag for the YSO candidates that are classified during the $JHK$[3.6][4.5] phase in order to separate faint YSO candidates from galaxies or background contaminants, and denoted the lower confidence in their identification by separating them from the rest of the YSO candidates. We identified 302 Class I YSOs and 1178 Class II/transition disk candidates in W51, and 917 Class I YSOs and 5187 Class II/transition disk candidates in W43. The difference in the YSO numbers is due to larger extent of W43. The color$-$color and color$-$magnitude diagrams of the identified YSOs, faint YSOs, and eliminated AGB star candidates are shown in Figure 4 and 5. The points are plotted without dereddening their photometry. In Table 5 we summarize the final source classification results.

Kang et al. (2009) used Spitzer and 2MASS data and classified 737 YSO candidates in W51 by using a combination of color criteria (Simon et al., 2007) and SED slope criteria (Lada, 1991; André et al., 1993; Greene et al., 1994). The class identification for the slope of log($\lambda$Fλ) vs. log($\lambda$) between $2$ and $\sim 20$ $\mu$m is:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

They used mainly the SED slope method and not the preferred method in Gutermuth et al. (2008). However, the Gutermuth et al. (2008) method has been slightly updated to better eliminate contaminant sources (Gutermuth et al., 2009), and we have also applied an additional color selection method as explained above to eliminate more background contaminants and AGB candidates. We cross-matched our catalog with that of Kang et al. (2009) using a matching radius of 1″and found 91% of the sources in common, the differences being due mainly to a slight difference of the cataloged region in both studies. Among the matched sources, 54% (366/680) of the sources are identified as YSO candidates in both catalogs, where 123 of them are in the YSO clusters we describe in Section 3. The rest of the sources are classified as contaminants or photospheres in our catalog. 17% (117/680) of them are not classified in our catalog, because they lack photometry in at least four bands, so they appear as unclassified. However, Kang et al. (2009) classified the sources based on their SED slopes of the four IRAC bands or the SED slope of log($\lambda$Fλ) vs. log($\lambda$) between $2$ and $\sim 24$ $\mu$m, without a constraint on the number of data points. 17% (119/680) of Kang et al. (2009) sources are classified as contaminants in our study, such as AGBs, PAH galaxies, PAH-dominated sources, or faint or possibly background sources (classified as YSOs in our first classification phase but then eliminated based on the Koenig & Leisawitz (2014) criteria). Finally, 12% (80/680) of the Kang et al. (2009) YSOs were classified as photospheres in our catalog since they do not have the infrared color excess that is described above for Class I or Class II candidates.

2.3.2 SED Slope Distribution

We calculated the SED slopes of the objects classified as YSOs using the IRAC-color criteria and determined their classifications using the SED slopes in order to compare the results from the two methods. We used all of the available photometry in the catalog except $J$ and $H$ and fit the slope of log($\lambda$Fλ) vs. log($\lambda$) (e.g., Figure 12) between $2$ and $\sim 20$ $\mu$m (Lada, 1987; André et al., 1993), including the $24$ $\mu$m data, in order to make a comparison between IRAC-identified embedded sources and transition disk candidates. We calculated the flux densities by using zero-point flux densities of 666.7 Jy for FKs (Cohen et al., 2003), and 280.9, 179.7, 115.0, 64.13 (IRAC Instrument Handbook, Version 2.0.3), and 7.17 Jy (MIPS Instrument Handbook, Version 3.0) for F3.6, F4.5, F5.8, F8.0, and F24, respectively. We used the isophotal wavelength values of 2.16 for [Ks], and 3.55, 4.44, 5.73, and 7.87 $\mu$m for [3.6], [4.5], [5.8], [8.0], respectively (Fazio et al., 2004), and 23.68 $\mu$m for [24] (MIPSGAL Data Delivery Version 2.0).

The distribution of the SED spectral index for both regions (for sources with R-squared values of the linear fit of $>$0.8) can be seen in Figure 7. In both the IRAC-color and SED slope methods, the YSO classifications refer to the observational parameters rather than the evolutionary phases and the choices of subdivisions are considered arbitrary (e.g., Dunham et al., 2014). When we compare the two, we find that the IRAC-color and SED slope methods give relatively consistent classifications for most objects. However, the comparison shows that using only the SED slope method does not help us to separate Class III or transition disk sources from Class II or embedded sources from Class I sources, while it does separate the Class III and Class II sources from the Class I or embedded sources.

Looking at the data for the IRAC-color Class I objects in W43, 98% of them have positive slopes as we expected. Similarly, we see that most of the embedded source candidates ($\sim$93%) have positive slopes, indicating that they are young sources. For the IRAC-identified Class II candidates, $\sim$70% of the SED slopes fall in the SED Class II range of -0.3 to -1.6. Most of them have negative slopes (89% in W43), indicating that they are evolved sources, and the 11% of them with positive slopes might be contaminant sources. Similarly, we see that 95% of transition disk candidates (Class II sources with disks with inner gaps) are consistent with being evolved sources with negative slopes. Approximately 11% of IRAC-identified transition disk candidates in W43 have positive slopes, which are most probably highly reddened transition disk candidates.

Overall, the histograms in Figure 7 show that, except for a small percentage of sources, the SED slope distribution is consistent with the YSOs identified using the IRAC-color method. In order to see the effect of extinction correction on classification, we used the $A_{K}$ values that were estimated during the classification process in Section 2.3.1, and dereddened the available photometry to calculate the SED slopes from $2$ to $24$ $\mu$m. While the extinction correction helps to reduce the small number of misclassifications of more evolved transition disk or Class II sources, it does not give a good result, especially for Class I or embedded sources; this is likely due to overcorrection. For IRAC-identified Class I candidates or embedded sources, extinction due to their envelopes and dusty disks plays an important role in the shape of their SEDs. Estimating an $A_{K}$ value and dereddening for these objects might be complicated, and it may cause misclassification due to their overcorrected colors and SED slopes. Indeed, the IRAC-color classification does not take into account the extinction for all sources except when double-checking Class II and Class I sources during Phase 2 if they are lacking detections at $5.8$ or $8.0$ $\mu$m); the comparison of IRAC-identified Class I sources according to the SED slopes without extinction correction gives more consistent results. The SED slope distribution in W51 yields similar results to W43, as can be seen in Figure 7. We provide the SED slope values that are calculated from the original photometric data with the errors and R-squared values in Table 13 and Table 14.

Both W51 and W43 are distant star-forming regions that contain significant amounts of high column density material, and young sources in these regions can be made extinct more easily by dust. In addition, there might be some misclassifications due to a mismatch of MIPS detections with their IRAC/near-IR counterparts. Due to the lower precision of MIPS astrometry, a transition disk candidate or a photosphere object can be classified as a deeply embedded source if it is mismatched with a MIPS source and lacks data in other bands. We have examined the catalog and images and removed obvious mismatched objects, but it is possible that chance alignments could lead to some that we cannot exclude with the available data.

3. Clustering Analysis

We used the Minimal Spanning Tree (MST) method (Cartwright & Whitworth, 2004) to examine the substructure in the W51 and W43 regions, such as groups/clusters of YSOs, without using any kinematic information, as we did for W49 (Saral et al., 2015). The MST method identifies clusters as a collection of stars that are connected to each other by short branches (with a maximum of cutoff distance of $d_{c}$) with a minimum number of members ($N$) in a group. It has been used to identify groups/clusters of YSOs (Koenig et al., 2008; Gutermuth et al., 2009; Beerer et al., 2010; Billot et al., 2011; Chavarría et al., 2014; Saral et al., 2015) and dense cores (Kirk et al., 2016) for nearby star-forming regions, and also to identify OB associations in other galaxies (Pietrzyński et al., 2001; Bastian et al., 2007).

There is no robust way to determine a branch length for the clustering determination. However, the most common way is to plot the distribution of branch lengths, fit straight lines through the long and short branch length domains, and choose the point of intersection (the Straight line Fit (SLF) method, see Figure 10 in Saral et al. (2015)). Choosing a branch length that falls between these two domains is also a good way to separate clusters from distributed sources (Gutermuth et al., 2008, 2009). Recently, Fischer et al. (2016) has identified the likely clusters in the Canis Major star-forming region by choosing the middle of the branch length range, which gave similar clusters. W51 and W43 are distant regions, and due to resolution effects and lower sensitivity, multiple stars may blend into single objects and therefore fewer stars might be detected. This effect causes the branch length value in the SLF method to increase, causing increased sensitivity to the largest-scale structures within the cloud. Therefore, it is likely that there are more clustered sources that are unresolved at this distance (e.g., Louvet et al., 2014).

To perform the clustering analysis we used a smaller field in both regions. This covers the main cloud structure but not the outer parts where there might be mainly line-of-sight contamination. The number of YSOs that we used for the clustering analysis is 857 and 5700 for W51, and W43, respectively. We used the MST method to perform a two-dimensional (2D) analysis by assuming that the projection effect would be small at these distances. However, in order to see the statistical significance of the identified MST groups in the 2D analysis, we performed Monte Carlo simulations of uniformly random distributed objects (see Section 3.3). Schmeja & Klessen (2006) also studied the effect of the 2D projections of the three-dimensional (3D) model clusters and concluded that although the clustering parameters can be different (and still giving similar mean values), the general behavior of the evolution of clusters is independent of the projection.

3.1. Clustering Results in W51

We first determined a branch length by using the SLF method, plotting the number of clustered objects versus branch lengths using a standard value of $N$ $\geq$ 10 for the minimum number of cluster members. Using a branch length of 82$\arcsec$.5 (which corresponds to 2.17 pc at the distance of W51), we found nine groups that contain 62% (536/857) of the total number of YSOs in the region. In order to investigate how clustering properties change with the requirement of minimum members in a group, we repeated the analysis with a mix of 5 to 15 members for a branch length of 82$\arcsec$.5. While we identified 26 MST groups with $N$ $\geq$ 5 and nine MST groups with $N$ $\geq$ 10, we identified only five MST groups with $N$ $\geq$ 15. When we decreased the number of minimum members from 15 to 5, we identified new small groups but the five largest MST groups remained unchanged.

To investigate how the numbers of identified groups change with different branch lengths, we plotted the number of groups containing 10 or more stars against cutoff distances from 1$\arcsec$ to 300$\arcsec$, with steps of 1$\arcsec$, in the left panel of Figure 8. A maximum of 16 MST groups can be identified using a branch length of 51$\arcsec$ (1.34 pc), while 33% (286/857) of the YSO candidates are associated with a cluster. The MST groups identified by both branch length values can be seen on the upper panel of Figure 9 and their properties are summarized in Table 6. The results show that using a shorter branch length allows us to identify subgroups and see the hierarchical structure. We also investigated the effects of changing the requirement of minimum members for 51$\arcsec$ and we identify 33 MST groups with $N$ $\geq$ 5, 16 groups with $N$ $\geq$ 10, and nine groups with $N$ $\geq$ 15. As expected with the increasing number of minimum members, we continue to identify the same largest groups, but lose the smaller ones. Our clustering analysis shows that the MST groups in W51 have 10$-$217 sources with diameters of 2.5$-$30 pc, and the subgroups within them have 10$-$44 sources with diameters of 1$-$7 pc (Table 8).

3.2. Clustering Results in W43

We identified the MST groups and subgroups within the W43 region by using the same method. We first determined a branch length with the SLF method by using standard value of $N$ $\geq$ 10 for minimum number of cluster members. 51 MST groups were identified with a branch length of 81$\arcsec$ (which corresponds to 2.20 pc at the distance of W43) and a total of 73% (4176/5700) of YSOs are associated with a cluster. We identified 134 MST subgroups with $N$ $\geq$ 5, 51 groups with $N$ $\geq$ 10, and 27 groups with $N$ $\geq$ 15. These 27 largest MST groups remained unchanged, while we identify new small groups with the decreasing number of minimum members.

As we did for W51, we investigated the change in the number of clusters with different branch lengths (see the right panel of Figure 8). A maximum of 95 MST groups can be identified with a branch length of 59$\arcsec$ (1.61 pc), while 43% (2448/5700) of the YSO candidates are associated with a cluster. The MST groups identified by both branch length values can be seen in Figure 11 and the properties of the five largest groups are summarized in Table 7. We identified 242 MST subgroups with $N$ $\geq$ 5, 95 MST groups with $N$ $\geq$ 10, and 53 MST groups with $N$ $\geq$ 15 for 59$\arcsec$; the largest 53 MST groups remained unchanged with any chosen minimum number of group members. Our clustering analysis shows that identified MST groups in W43 have 10$-$3047 sources with diameters of 1.4$-$142 pc and subgroups within them have 10$-$217 sources with diameters of 0.8$-$14 pc (Table 9).

3.3. Random Clustering Analysis

To determine the statistical significance of the identified MST groups, we performed Monte Carlo simulations of uniformly random distributed objects in order to examine their clustering properties and compare them with the observed ones. We created 1000 realizations of 857 objects spread randomly in a $0.5\times 1.3$ degree field to match the number of YSOs and the size of the W51 region. We used both the branch length determined from the MST method and the branch length giving the maximum number of groups (subgroups) in order to compare both results with the observed groups. The largest random group among all of the simulations has 58 members, with an average maximum group size of $24\pm 6$ members. In comparison, the observed largest group in W51 is significantly larger, with 215 members. The second largest random group in each of the simulations has on average $19\pm 4$ members, with a maximum of 36, while the observed second largest group is still larger, with 172 members. Similarly, the average sizes of third, fourth, and fifth largest groups in each simulation are smaller than the observed ones. However, the smaller observed groups (sixth to ninth largest) have similar sizes to the average size of random groups, which indicates that they might be random associations. For the subgroups identified with the branch length of 51$\arcsec$, we see that all of the observed subgroups are larger than the mean size of random subgroups ($8\pm 4$), which indicates that they are different than random distributions. Similarly, we generated 1000 simulations of 5700 randomly distributed sources in a $2.8\times 2.1$ degree field, which is the size of the W43 region. Almost all of the observed groups and subgroups are larger than the average random group size, ($8\pm 4$).

As a final step we used a two-sided Kolmogorov-Smirnov (K-S) test (Eppler et al., 1999) to determine if the two datasets (observed and random clusters) differ significantly. We calculated the two-sided K-S statistic between each of the 1000 realizations of randomly distributed sources and the mean cumulative distribution function (CDF) of random clusters from all 1000 realizations together. The distribution of the K-S statistic values, for the branch lengths of $d_{c}$ = 82.5$\arcsec$ and $d_{c}$ = 51$\arcsec$, are shown in Figure 10 as green and blue histograms, respectively. In the same figure, the K-S statistic value, between the CDF of observed MST groups and the mean CDF of random clusters, is shown with a vertical dashed line. We find that for both branch lengths, the K-S statistic value between the observed group size CDF and the mean random one is well outside of the distribution of K-S statistic values between the mean random CDF and those of each of the 1000 realizations, Thus, we infer that the probability of the distribution of a real MST group being drawn from a random distribution is less than $10^{-3}$.

4. SED Models for MYSOs

We used SED models to identify the MYSO candidates in the W51 and W43 regions since it is difficult to identify them using colors. We used the Python package of the SED Fitting tool by Robitaille et al. (2007) and SED models from Robitaille et al. (2006) to model the available photometric data (2MASS/UKIDSS $J$, $H$, and $K_{s}$, Spitzer/IRAC $3.6$, $4.5$, $5.8$, and $8.0$ $\mu$m, and Spitzer/MIPS $24$ $\mu$m). This best-fit method is described in detail by Robitaille et al. (2007). Only objects with at least six data points were fitted, and we used $e^{-\chi^{2}/2}$ values as weights. We chose the fits that satisfy the criterion of (${\chi}^{2}$ - ${\chi}^{2}_{best}$)/$n_{data}$¡ 3 as in Robitaille et al. (2007). Only objects with a minimum of 10 fits satisfying the defined best-fit criteria were used when calculating the weighted mean values for the physical parameters. We allowed the $A_{V}$ to vary between 0 and 40 mag, which gives enough flexibility and prevents unphysical solutions, and is also consistent with the $A_{V}$ values for MYSOs that are forming UCH IIs (e.g., Hanson et al., 2002). We also let the distance to vary within 10% ranges of 5.4 kpc for W51, and 5.6 kpc for W43. For sources without a $24$ $\mu$m detection (except for sources in the saturated or high-background regions), we used an upper limit based on the completeness estimate at $24$ $\mu$m by Gutermuth and Heyer (2015) in order to better constrain the SED. The luminosity and temperature are directly derived from the SED, and then mass and age values are determined using the evolutionary tracks in a Hertzprung-Russel diagram. The lack of photometric points at wavelengths longer than $24$ $\mu$m means that the disk and envelope parameters are poorly constrained. The total dust mass, for example, can only be estimated if the location of the far-IR bump is known.

We have applied the Bayesian method described in Azimlu et al. (2015) to construct posterior probability distribution functions for the model parameters. We assume that the measurement errors in the photometry are normally distributed, which results in a likelihood function proportional to the $e^{-\chi^{2}/2}$ distribution. We use uniform priors for all of the model parameters. The resulting posteriors are sampled using a Markov Chain Monte Carlo (MCMC) method. Within the set of YSOs we selected the ones with $L$ $\geq$ $10^{3}$L⊙ from both methods. We give the weighted mean values for the physical parameters such as masses and luminosities from the Robitaille et al. (2006) models, and to provide some indication of the model accuracy, we report the standard deviation values in Table 10 and Table 11. These model grids do not include binary or multiple system cases, while the multiplicity in massive stars and in star clusters is known to be higher than that of young low-mass stars (Duchêne et al., 2001; Zinnecker & Yorke, 2007). We should also note that the Robitaille models were developed mainly for lower-mass protostars. Therefore, the physical parameters of the MYSO candidates given here should be used cautiously. However, De Buizer et al. (2016) compared the fit results from the Robitaille method to a Turbulent Core model-based method (Zhang & Tan, 2011) that was developed for intermediate- and high-mass stars, and they found that, although the accretion rates are significantly underestimated, the stellar masses are consistent. Moreover, Mottram et al. (2011) found that for sources with $L$ $>$ $10^{3}$L⊙, luminosities are better determined with the SED fitter compared to graybody fitting. However, MYSO candidates determined with luminosities over $10^{3}$L⊙ should still be confirmed by far-IR observations and spectroscopic data.

As a result, we identified 17 and 14 MYSO candidates in W51 and W43, respectively. In W51, 16 of the MYSO candidates are located within the clusters, and one of them is a distributed source. Seven of these MYSO candidates were classified as Class I, eight as Class II candidates, one as a transition disk candidate, and one as an embedded source candidate with IRAC-color classification. In W43, 11 of the MYSO candidates are located within the clusters and the rest are distributed outside. Among the MYSO candidates, nine were classified as Class I candidates and five were classified as Class II candidates using the IRAC-color classification. The distribution of MYSO candidates is shown in Figure 14 and Figure 15 and the SED plots for those having $24$ $\mu$m data are shown in Figure 12 and Figure 13.

5. Discussion

5.1. Massive Star Formation Tracers

5.1.1 UCH II Regions and Masers in W51

W51 hosts several H II regions: G49.5-0.4 and G49.4-0.3 associated with W51A, and G49.2-0.3, G49.1-0.4, G49.0-0.3, and G48.9-0.3 associated with W51B (Kumar et al., 2004) (Figure 14). In addition, there are several compact and UCH II regions that are not detected at ¡20$\mu$m because of heavy extinction (Gaume et al., 1993; Zhang & Ho, 1995, 1997).

W51A, the most active part of the region, hosts several protoclusters and UCH II regions (Gaume et al., 1993; Goldader et al., 1994; Zhang & Ho, 1995; Barbosa et al., 2008). G49.5-0.4 is the most luminous part in W51A and hosts two embedded protoclusters strongly correlated with radio continuum emission: IRS 1 and IRS 2 (Wood & Churchwell, 1989; Wynn-Williams, Becklin and Neugebauer, 1974). In the IRS 2 UCH II region, a stellar object that is not resolved in our images was identified by Figuerêdo et al. (2008) and confirmed as an early O-type object with spectroscopy by Barbosa et al. (2008). We identified one subgroup (2a in Table 6) that covers both protoclusters including the H II region G49.5-0.4 as seen in Figure 14. Half of the YSO candidates of subgroup 2a are positionally associated with the G49.5-0.4, and the other half surrounds the H II region with a filamentary shape. We identified two MYSO candidates, SSTOERC$\_$G049.4885-0.3628 and SSTOERC$\_$G049.4821-0.3507 (detected also at $24\mu$m), that are shown in Figure 14. The other main H II region in W51A is G49.4-0.3, which seems to be associated with subgroup 2f, where we identified 11 YSO candidates, four of which are MYSO candidates: SSTOERC$\_$G049.3820-0.2563, SSTOERC$\_$G049.3800-0.2696, SSTOERC$\_$G049.3884-0.2634 and SSTOERC$\_$G049.3793-0.256.

In W51B, the largest and brightest H II region is G49.2-0.3, which is indicated as an interacting region with the W51C supernova remnant (Green et al., 1997; Koo & Moon, 1997). We identified a MYSO candidate in the vicinity of this region, SSTOERC$\_$GO49.2213-0.3354, and five more YSO candidates. The other three H II regions, G49.1-0.4, G49.0-0.3, and G48.9-0.3, seem positionally but partially associated with subgroups 1d, 1a, and 1e, respectively, and are located along a filamentary structure parallel to the Galactic plane (right panel of Figure 14). We identified one MYSO candidate, SSTOERC$\_$G048.9916-0.3143, in the field of G49.0-0.3 H II region, which is a weaker H II region in W51B. And finally, G48.9-0.3, a strong UCH II region in W51B, is partially corresponds to the subgroup 1e and is also located at the edge of the subgroup 1c, which hosts the MYSO candidate, SSTOERC$\_$GO48.9022-0.2734 in its vicinity.

W51 also hosts many masers (Zhang & Ho, 1995; Imai et al., 2002; Fish and Reid, 2007; Henkel et al., 2013). Among them are methanol masers, which are known as good massive star formation tracers (Pestalozzi et al., 2002; Breen et al., 2013). Two main methanol maser regions toward W51A (G049.49-00.37 and G49.49-0.39; Caswell et al., 1995; Minier et al., 2000; Etoka et al., 2012) are shown in Figure 14. There are tens of methanol masers in that region clustered within a few arcsec, and most of them do not have known counterparts with any IR source. This indicates that these sources are in very early evolutionary stages that we cannot yet trace with IR photometry.

5.1.2 Compact Fragments and Masers in W43

W43 hosts a giant H II region dominated by its centimeter and infrared emission, which makes the detection of compact sources with associated UCH IIs difficult (Motte et al., 2003). Using 1.3 mm and 350 $\mu$m continuum emission, Motte et al. (2003) identified 50 compact fragments in the W43-Main region with sizes ranging from 0.09 to 0.56 pc, which tens of fragments of which are protoclusters (see Louvet et al., 2014). Motte et al. (2003) suggested that the fragments MM1, MM2, and MM3, associated with masers and outflows, are forming high-mass stars; Bally et al. (2010) confirmed this with High-GAL far-infrared observations. These sources should be in very early stages (Class 0-like) and be still very embedded. Consequently, current observations have not yet found many associated near infrared sources within the fragments. Recently, 20 high-mass protostars were found in the MM1 ridge, encompassing MM1 and MM9 (Motte et al. 2017, in preparation). We did not find any MYSO candidates associated with these fragments. However, in subgroup 3 we see a YSO candidate coincident with each of the following fragments: MM4, MM11, MM14, MM15, and MM20. Motte et al. (2003) indicates that MM11 is also a methanol source with no infrared or centimeter detection; however, we see a Class I candidate (SSTOERC$\_$GO30.7598$-$0.0525) that coincides with the filament. We also see a YSO candidate positionally associated with the fragments MM9 and MM16 fragments that corresponds to subgroup 62 (Figure 15), which could become more massive in the future given the large amount of gas surrounding it.

5.2. Star Formation History in W51

In Section 3, we identified YSO groups and subgroups by using the MST method. In order to understand the star formation history in the region, we should consider the relative ages of these groups. Assuming Class I sources represent relatively early stages of formation compared to Class II sources, we can use the ratio of Class II/I sources to compare the relative ages of different groups within the region.

The Class II/I ratio is $2.50\pm 0.25$ when we consider all of the MST groups identified with $d_{c}$ = 82$\arcsec$.5(2.17 pc; see Table 8). If we only consider the subgroups identified with $d_{c}$ = 51$\arcsec$, we find a Class II/I ratio of $2.06\pm(0.27)$, which is consistent with the ratio from the SLF method within the uncertainties. If we look at individual groups within the region, the biggest group, G49.01-0.31, consists of seven subgroups (the red group in Figure 9, corresponds to W51B) and the ratio of Class II/I is $2.98\pm 0.48$; if we consider only the subgroups, it is $2.56\pm 0.54$. Meanwhile, for the second largest group, G49.45-0.37 (the blue group in Figure 9, corresponds to W51A), the ratio of Class II/I is $1.94\pm 0.32$, and it consists of six subgroups with a Class II/I ratio of $1.65\pm 0.44$. Therefore, either considering the large MST groups or subgroups, G49.45-0.37 (where W51A is located) is relatively younger than G49.01-0.31 (where W51B is located). This is consistent with the Ginsburg et al. (2015) study that showed a low star formation potential of W51B, with a density of $n\sim 1\times 10^{4}$ $\textrm{cm}^{\textrm{-3}}$. Furthermore, they associated the high-mass density, n $\gtrsim$ $5\times 10^{5}$ $\textrm{cm}^{\textrm{-3}}$ in W51A with high-mass star formation.

Kumar et al. (2004) suggested that the star formation in W51 is triggered by galactic spiral density waves. Also, the interaction between the H II region G49.2-0.3 and the W51C supernova remnant W51C ($\sim 3\times 10^{4}$ years old; Koo et al., 1995) might be triggering a future star formation episode (Green et al., 1997; Koo & Moon, 1997). Furthermore, the existence of the supernova remnant and old, diffuse H II regions in W51B can be considered as signposts of older star-forming regions. We did not identify a subgroup corresponding to the H II region G49.2-0.3, but we identified the larger-scale cluster 4 with a Class II/I ratio of 1.70(0.39), which indicates a relatively younger age than W51B that is consistent with a triggered star formation scenario.

5.3. Star Formation History in W43

The class II/I ratio is found to be $4.58\pm 0.19$ for the W43 region when we consider all of the MST groups identified with d_{c}$$=81$\arcsec$ (Table 9). If we only consider the subgroups, we find a smaller ratio, $3.71\pm 0.19$, but it is still closer to the value from the SLF method. For just MST group 1 (the big red group in Figure 11), corresponding to the whole molecular cloud, the II/I ratio is $4.20\pm 0.20$. We see that group 1, identified using SLF method, consists of 62 subgroups with a Class II/I ratio of $4.19\pm 0.25$. For both the SLF method and the maximum number of subgroups, W43 has significantly larger class II/I ratio than W51. The findings of many YSOs and subgroups are also consisted with the substructures found by Carlhoff et al. (2013) and supporting the hypothesis of several HI-to-H2 transition layers that had been suggested by Bialy et al. (2017).

5.4. Comparison of Other Star-Forming Regions

W51, W43, and W49 are among the brightest star-forming regions in the Galaxy. Although they have similar gas masses and stellar masses (see Section 1), they have different morphologies and clustering properties that suggest a different star-formation history. W49 shows a hierarchical network of filaments at scales from $\sim$10 to 100 pc, suggesting that it is formed from global gravitational contraction (Madrid et al., 2013); it also hosts a large central cluster with tens of massive YSO candidates which might be triggering the star formation. On the other hand, both W51 and W43 show a younger and an older part, with more distributed clusters and distributed MYSO candidates within the cloud.

W51 has an extended structure that shows independent sites of star formation (Clark et al., 2009); on the other hand, Kumar et al. (2004) suggested that the star formation in W51 is triggered by galactic spiral density waves. Clark et al. (2009) derived an age range of 3-6 Myr (assuming the distance as 6 kpc) for a very luminous supergiant, suggesting at least 5 Myr of massive star formation in W51. On the other hand, by studying the old Wolf-Rayet stars Bik et al. (2005) determined an age range of 5-10 Myr. We have identified one large cluster, representing both W51A and W51B regions in the GMC with different Class II/I ratios, that suggests a relative age difference. In addition, the interaction between the H II region G49.2-0.3 and the W51C supernova remnant (Green et al., 1997; Koo & Moon, 1997) might trigger future star formation. Similarly, W43 extends over 140 pc and is surrounded by an atomic gas envelope with a diameter of 290 pc (see Section 1). Considering the similar distance of both regions, the class II/I ratios show that W51 is relatively younger, with the Class II/I ratio of 2.50 (0.25) than W43 which has a Class II/I ratio of 4.58(0.19). However, we should note that although the number of Class IIs are over-numerous in W43, there is a larger number of Class 0s and even Class Is which we cannot probe at near-IR.

Among these three regions, W49 might be the youngest, with a suggested age range of 1-2 Myr (Homeier & Alves, 2005). W51 is also relatively young, with a Class II/I ratio of 2.50. Even though W49 has a smaller Class II/I ratio, since W49 is a much more distant region it does not allow us to make a direct comparison as we discussed in Saral et al. (2015). In addition, we should note that several SF events are developing in these clouds. For instance, W43-Main could also be at an earlier stage than W49 since it has a very compact massive cluster, formed by the present starburst cluster, in its center. This compact massive cluster should be observed by the James Webb Space Telescope in order to fully resolve and detect the individual sources in the cluster.

In addition, we see that the fraction of YSOs in clusters is slightly higher in W43 (43-73%) than in W51 (32-60%). Clustered fractions were reported to be 53-67% in the Cygnus-X (Beerer et al., 2010), 44-70% in W5 (Koenig et al., 2008), 53-67% and 44-76% in G305 and 54-81% in G333, respectively (Saral et al., 2015). Although we can say that all of these regions have similar clustered fractions for all masses of YSOs, the clustered fractions for MYSOs are higher (94% in W51, 79% in W43, 100% in W49) which supports the theory that massive stars are being formed in clusters.

When we compare the branch length distances shown in Table 12, we can see that both W51 (82$\arcsec$.5(2.17 pc) and W43 (81$\arcsec$(2.20 pc) have similar branch lengths, suggesting that the density of clusters is the same. These differ from G305 and also G333 (Saral et al., 2015), which have smaller branch length distances, indicating denser clusters than those of W51, W43, and W49. W49 has the largest branch length; however, this is likely due to its larger distance and incompleteness effects. For example, even though W49 might be intrinsically denser than the W51 and W43 regions, it has a larger branch length (96$\arcsec$ or 5.2 pc) than W51 and W43, which indicates that the clusters in W49 are less dense. However, this might be due to the larger distance and incompleteness effects. In order to make a direct comparison, one needs to rescale the source positions and magnitudes in order to determine the projected appearances of these regions at 11 kpc.

As we did in Saral et al. (2015) for G305 and G333, by following the same strategy we determined the projected appearance of the W51 and W43 YSOs by rescaling the relative source positions and magnitudes for a distance shift from 5.4 kpc (W51) and 5.6 kpc (W43) to the 11 kpc distance of W49. Following the rescaling, we performed the source classification and clustering as described in Section 2.3 and Section 3. At reprojected distances we identified 458 and 1524 YSO candidates in W51 and W43, respectively, and branch lengths of $d_{c}$ = 120$\arcsec$.94 (6.47 pc) with five clusters in W51, and $d_{c}$ = 91$\arcsec$.82 (4.91 pc) with eight clusters in W43. As seen from the clustering results for the reprojected regions, the clustering analysis yields higher branch lengths due to larger distances and incompleteness effects. The larger derived branch length is a result of our inability to resolve the smaller distances between the closest sources. The rescaled W51, which has the largest branch length of $d_{c}$ = 120$\arcsec$.94 (6.47 pc), seems less dense than other regions. The rescaled W43 region, with a branch length of $d_{c}$ = 91$\arcsec$.82 (4.91 pc), is denser than both W49 and rescaled W51 but less denser than G305 and G333. Among these regions, G305 has the densest clusters, with a branch length of $d_{c}$ = 52$\arcsec$ (0.90 pc). The branch length and clustered fraction values are summarized in Table 12.

6. Summary

In this paper, we generated photometric catalogs for W51 and W43 with a wavelength range of 1.2 to 24 $\mu$m by combining 2MASS/UKIDSS, Spitzer IRAC, and MIPS data in order to identify and classify the YSO candidates. We identified 302 Class I candidates, 1178 Class II/transition disk candidates, and 56 embedded sources in W51, and 917 Class I candidates, 5187 Class II/transition disk candidates, and 144 embedded sources in W43. In addition, we compared the color classification method to the SED slope method and concluded that except for a small percentage of sources, both methods are consistent. Following the YSO classification, we identified groups and clusters based on their spatial distributions using the MST method. We identified nine clusters in W51, where 60% of the YSO candidates (including transition disks) belong to the clusters, each with a minimum of 10 members. We identified 51 clusters in W43 (over six times more than in W51), where 73% of the YSO candidates belong to the clusters. Similar branch length distances ($\sim$2 pc) in W51 and W43 indicate similar densities within clusters in both regions. However, with a larger Class II/I ratio, W43 (4.58) is a relatively older region than W51 (2.50), and W43 has at least six times more clusters. Based on the morphology of the clusters, which are distributed within the cloud and show relative age differences according to their Class II/I ratios and also the existence of Class 0 clusters studied in the FIR and sub-millimeter, there might be several sites of star formation that are independent of one another in terms of their ages, physical conditions, and evolutionary progress. We performed 1000 Monte Carlo simulations of randomly distributed objects in order to compare the properties of random clusters to real ones, and also applied a K$-$S test in order to test the statistical significance of the identified clusters. We find that for both regions there is a high probability that the identified MST groups are different than expected from a random distribution.

We identified 17 and 14 MYSO candidates ($L$ $\geq$ $10^{3}$L⊙) in W51 and W43, respectively, according to both the SED Fitting tool by Robitaille et al. (2007) and the Bayesian method by Azimlu et al. (2015). The clustered fraction of MYSOs are $\sim$94% and $\sim$79% in W51 and W43, respectively. In both W51A and W51B, we see subgroups with MYSO candidates associated with H II regions, which makes these sources strong candidates for follow-up spectroscopic observations. We do not see any of our MYSO candidates associated with previously identified dense fragments in W43, which might be a sign of several independent star formation events and bursts, while we see YSO candidates in subgroup 3 that are associated with dense fragments. These sources are good candidates for follow-up studies.

We compared W51 and W43 to each other and also to other star-forming regions such as G305, G333, and W49. While W51 is a smaller star-forming region than W43, W51 and W43 have many similarities, such as distributed clusters in the younger and older parts of each cloud, suggesting independent star formation throughout the clouds. Both W51 and W43 have similar YSO densities in clusters, while G305 and G333 have smaller branch length distances, indicating denser clusters. Compared to the others discussed in this study and W49, G305 and G333 are small star-forming regions (Saral et al., 2015). W49 is smaller than both W51 and W43 and has more centrally concentrated cluster formation.

This work is based in part on observations made with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA. Support for this work was provided by NASA. G.S. acknowledges support from NASA grant NNX12AI60G, Istanbul University grant BAP50195, Marie Sklodowska-Curie RISE grant ”ASTROSTAT” (project number 691164), and Swiss National Science Foundation Sinergia project ”STARFORM” (project number CH: CRSII2_160759). H.A.S and R.M.G. acknowledge partial support from NASA grants NNX10AD68G and NNX12AI55G. The authors thank the referee and Kathleen Kraemer for careful reading of the manuscript and valuable comments that helped improve the paper.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 National Aeronautics and Space Administration and the National Science Foundation and NASAs Astrophysics Data System Bibliographic Services. Facilities: 2MASS ($JHK_{S}$), Spitzer (IRAC, MIPS), UKIRT ($JHK_{S}$)