Massive 70 micron quiet clumps I: evidence of embedded low/intermediate-mass star formation activity

TL;DR

This study investigates 18 massive, 70 micron quiet clumps to understand their potential for star formation, revealing that many already host faint protostars and exhibit infall signatures indicative of ongoing accretion processes.

Contribution

The paper provides new observational evidence that most starless candidates already contain faint protostars and show infall motions, challenging previous assumptions about the earliest stages of massive star formation.

Findings

Half of the clumps embed faint 24 micron sources.

Infall signatures are detected in seven clumps, with estimated accretion rates.

Most clumps may already host low-mass protostars, influencing massive star formation theories.

Abstract

Massive clumps, prior to the formation of any visible protostars, are the best candidates to search for the elusive massive starless cores. In this work we investigate the dust and gas properties of massive clumps selected to be 70 micron quiet, therefore good starless candidates. Our sample of 18 clumps has masses 300 < M < 3000 M_sun, radius 0.54 < R < 1.00 pc, surface densities Sigma > 0.05 g cm^-2 and luminosity/mass ratio L/M < 0.3. We show that half of these 70 micron quiet clumps embed faint 24 micron sources. Comparison with GLIMPSE counterparts shows that 5 clumps embed young stars of intermediate stellar mass up to ~5.5 M_sun. We study the clump dynamics with observations of N2H+ (1-0), HNC (1-0) and HCO+ (1-0) made with the IRAM 30m telescope. Seven clumps have blue-shifted spectra compatible with infall signatures, for which we estimate a mass accretion rate 0.04 < M_dot < 2.0 x 10^-3 M_sun yr^-1, comparable with values found in high-mass protostellar regions, and free-fall time of the order of t_ff = 3 x 10^5 yr. The only appreciable difference we find between objects with and without embedded 24 micron sources is that the infall rate appears to increase from 24 micron dark to 24 micron bright objects. We conclude that all 70 micron quiet objects have similar properties on clump scales, independently of the presence of an embedded protostar. Based on our data we speculate that the majority, if not all of these clumps may already embed faint, low-mass protostellar cores. If these clumps are to form massive stars, this must occur after the formation of these lower mass stars.

Tables

Table 1. Table 1: Photometry results of the 18 clumps studied in this work. Col.1: Clump name; Col.2: source id number as in Traficante et al . ( 2015a ) catalogue; Cols. 3 − - 4: Coordinates of the clump centroids obtained from the Gaussian fit done at 250  μ 𝜇 \mu m; Cols. 5 − - 10: Clump fluxes at 160, 250, 350, 500, 870 and 1100  μ 𝜇 \mu m respectively; Cols. 11 − - 13: minimum, maximum FWHMs and PA of the 2d-Gaussian fit. When the fit gives FWHM min =FWHM max the source is circular and the PA if fixed to 90 ∘ superscript 90 90^{\circ} ; Col. 14: deblend parameter. 1 means that one (or more) source companion has been deblended before measuring the clump flux.

Clump	Source	RA	Dec	F160μm	F250μm	F350μm	F500μm	F870μm	F1100μm	FWHMmin	FWHMmax	PA	Deblend
		(∘)	(∘)	(Jy)	(Jy)	(Jy)	(Jy)	(Jy)	(Jy)	(″)	(″)	(∘)
15.631-0.377	1	18:20:29:1	-15:31:26	1.21	3.98	4.66	2.68	0.61	0.29	28.37	36.00	126.59	0
18.787-0.286	1	18:26:15.3	-12:41:33	9.50	32.76	28.67	15.41	3.40	1.17	29.58	36.00	223.89	0
19.281-0.387	1	18:27:33.9	-12:18:17	11.38	25.53	21.38	8.38	2.23	0.88	36.00	36.00	90.0	1
22.53-0.1921	1	18:32:59.7	-09:20:03	14.61	33.21	22.02	8.90	2.43	…	22.81	36.00	263.65	0
22.756-0.284	1	18:33:49.1	-09:13:04	10.66	21.25	15.23	7.10	1.36	0.67	18.04	36.00	237.51	0
23.271-0.263	7	18:34:38.0	-08:40:45	14.77	34.53	21.09	7.10	1.93	1.16	27.36	30.00	232.83	0
24.013+0.488	1	18:33:18.5	-07:42:23	10.81	33.36	32.25	13.31	1.97	1.04	28.86	36.00	262.17	1
24.528-0.136	1	18:36:31.0	-07:32:24	11.87	20.93	21.12	13.80	3.24	1.38	27.88	36.00	193.65	0
25.609+0.2281	3	18:37:10.6	-06:23:32	10.13	26.27	29.51	17.11	2.34	…	36.00	36.00	90.0	1
25.982-0.0561	1	18:38:54.5	-06:12:31	11.84	24.95	16.92	7.05	1.68	…	30.51	36.00	251.49	1
28.178-0.091	2	18:43:02.7	-04:14:52	34.49	71.11	47.44	29.17	2.79	1.10	29.58	36.00	115.91	0
28.537-0.277	1	18:44:22.0	-04:01:40	3.01	19.65	15.67	5.78	1.81	0.46	21.48	36.00	143.86	1
28.792+0.141	2	18:43:08.8	-03:36:16	2.44	11.79	6.49	3.10	0.69	0.30	27.05	27.05	90.0	1
30.357-0.837	2	18:45:40.6	-02:39:45	9.70	15.13	11.89	6.48	1.00	0.31	28.46	36.00	217.28	0
30.454-0.135	1	18:47:24.0	-02:16:01	13.15	18.72	15.76	8.17	2.85	0.90	36.00	36.00	90.0	0
31.946+0.076	2	18:49:22.2	-00:50:32	7.37	13.90	13.81	8.24	1.70	0.91	26.41	36.00	262.13	0
32.006-0.51	1	18:51:34.1	-01:03:24	1.88	8.44	9.09	3.83	0.54	0.31	31.72	36.00	212.56	1
34.131+0.075	2	18:53:21.5	+01:06:14	13.59	30.01	19.09	6.70	1.51	0.49	27.91	36.00	158.51	0

Table 2. Table 2: Sensitivity and FWHMs comparison between the surveys used in this work: MIPSGAL, Hi-GAL, ATLASGAL and the BGPS.

Survey	Wavelength	Sensitivity	FWHM
	( $\mu$m)	(MJy/sr)	(″)
MIPSGAL	24	2.71	5.9
Hi-GAL	70	272	10.0
Hi-GAL	160	702	13.5
Hi-GAL	250	162	18.0
Hi-GAL	350	62	24.0
Hi-GAL	500	62	34.5
ATLASGAL	870	73	19.2
BGPS	1100	2.54	33.0

Table 3. Table 3: Clumps properties derived from the SEDs. Col.1: Clump name; Col.2: Clumps mass with errors estimated from the SED fit and assuming a flux uncertainties of 20% for Hi-GAL fluxes (Traficante et al . , 2015a ) and of 40% on ATLASGAL and BGPS fluxes, plus an uncertainties on distance estimation of 10%; Col.3: Bolometric luminosity estimated integrating the flux in the range 160-1100  μ 𝜇 \mu m; Col. 4: Source radius derived from the geometric mean of the FWHMs in Table 1 ; Col. 5: Temperatures and associated uncertainties; Col. 6: χ 2 superscript 𝜒 2 \chi^{2} of the SED fits done using Hi-GAL, ATLASGAL and, where available, BGPS datasets; Col.7: Source distance taken from Traficante et al . ( 2015a ) .

Clump	Mass	Luminosity	Radius	$\Sigma$	Temperature	${\chi }^2$	Distance
	(M⊙)	(L⊙)	(pc)	(g cm-2)	(K)		(kpc)
15.631-0.377	269(79)	18	0.54(0.05)	0.06(0.01)	9.6(0.4)	0.60	3.47(0.35)
18.787-0.286	1915(550)	206	0.69(0.07)	0.27(0.06)	10.3(0.4)	2.96	4.36(0.44)
19.281-0.387	701(206)	123	0.67(0.07)	0.10(0.02)	11.4(0.5)	1.61	3.82(0.38)
22.53-0.192	1579(488)	349	0.80(0.08)	0.16(0.04)	11.8(0.6)	2.01	5.77(0.58)
22.756-0.284	655(194)	136	0.55(0.06)	0.14(0.03)	11.8(0.5)	0.67	4.43(0.44)
23.271-0.263	997(297)	285	0.72(0.07)	0.13(0.03)	12.3(0.6)	4.76	5.21(0.52)
24.013+0.488	1957(559)	294	0.81(0.08)	0.20(0.04)	10.9(0.4)	5.83	5.18(0.05)
24.528-0.136	2074(650)	215	0.80(0.08)	0.22(0.05)	10.6(0.5)	1.34	5.19(0.52)
25.609+0.228	3098(953)	299	0.97(0.10)	0.22(0.05)	10.3(0.4)	2.81	5.57(0.56)
25.982-0.056	889(275)	200	0.80(0.08)	0.09(0.02)	11.9(0.6)	1.33	5.00(0.50)
28.178-0.091	2092(610)	669	0.85(0.09)	0.19(0.04)	12.6(0.6)	8.72	5.35(0.54)
28.537-0.277	1153(327)	145	0.67(0.07)	0.17(0.03)	9.9(0.4)	12.63	4.96(0.50)
28.792+0.141	449(128)	72	0.61(0.06)	0.08(0.02)	10.6(0.4)	6.46	4.62(0.46)
30.357-0.837	372(111)	98	0.67(0.07)	0.06(0.01)	12.6(0.6)	3.46	4.30(0.43)
30.454-0.135	1326(419)	233	1.0(0.10)	0.08(0.02)	11.7(0.6)	1.40	5.84(0.58)
31.946+0.076	1431(440)	156	0.82(0.08)	0.14(0.03)	10.7(0.5)	0.59	5.51(0.55)
32.006-0.51	449(127)	50	0.70(0.07)	0.06(0.01)	10.0(0.3)	9.39	4.24(0.42)
34.131+0.075	480(142)	108	0.55(0.06)	0.11(0.02)	12.0(0.6)	0.86	3.56(0.36)

Table 4. Table 4: Emission of the 24  μ 𝜇 \mu m source for the clumps with a faint 24  μ 𝜇 \mu m counterpart identified. Col.1: Clump name; Col. 2-3: Coordinates of the 24  μ 𝜇 \mu m source counterpart; Col. 6: Integrated flux of the 24  μ 𝜇 \mu m source. The flux has been estimated in a radius of 3.5″ and corrected for an aperture correction of 2.78.

Clump	24 $\mu$m RA	24 $\mu$m Dec	24 $\mu$m flux
	(∘)	(∘)	(mJy)
18.787-0.286_11	18:26:15.2	-12:41:30	28.2
18.787-0.286_21	18:26:15.4	-12:41:38	4.45
19.281-0.387	18:27:33.9	-12:18:23	18.68
22.53-0.192	18:32:59.7	-09:19:59	29.31
24.013+0488_1	18:33:18.6	-07:42:24	9.46
24.013+0488_2	18:33:18.1	-07:42:33	18.62
24.013+0488_3	18:33:17.6	-07:42:44	13.11
28.178-0091	18:43:02.2	-04:14:32	24.08
30.357-0.837_1	18:49:40.9	-02:39:47	15.61
30.357-0.837_2	18:49:40.4	-02:39:50	5.49
30.454-0.1352	18:47:23.9	-02:15:55	9.27
31.946+0.076	18:49:22.3	-00:50:35	19.49
34.131+0.075_1	18:53:21.1	01:06:13	8.95
34.131+0.075_2	18:53:21.1	01:06:23	6.13

Table 5. Table 5: Physical parameters and classification of the central star in each clump with NIR/MIR counterparts. Col 1: Clump name; Col. 2: Source classification obtained as described in the text; Col. 3: Average magnitude; Col. 4: Average column density, obtained from A V using the conversion described in Bohlin, Savage & Drake ( 1978 ) ; Col. 5: Average mass of the central star; Col.6: χ 2 superscript 𝜒 2 \chi^{2} range of the models used to estimate the weighted parameters as described in the text; Col. 6: Number of models used to estimate the physical parameters. The average values and uncertainties associated with the estimation of A V , column density and mass of the central star have been evaluated as the weighted mean and variance of the distributions with weights equal to the inverse of the χ 2 superscript 𝜒 2 \chi^{2} value of each fit.

Clump	Class	AV	$\Sigma$	M∗	${\chi }^2$	Models
		(mag)	(g cm-2)	(M⊙)
18.787-0.286_1	I	56.4$\pm$4.6	12.7$\pm$1.0	4.65$\pm$0.62	7.69$-$16.42	50
18.787-0.286_2	I	75.9$\pm$16.9	17.0$\pm$3.8	2.87$\pm$0.98	0.38$-$5.00	1178
19.281-0.3871	I-II	66.5$\pm$8.3	14.9$\pm$1.9	4.60$\pm$0.55	0.36$-$5.00	73
22.53-0.192	I	93.4$\pm$9.3	21.0$\pm$2.1	5.52$\pm$0.73	8.68$-$11.00	50
24.013+0.488_1	I	63.0$\pm$2.2	14.1$\pm$0.5	4.15$\pm$0.35	0.02$-$8.10	50
24.013+0.488_2	I	69.6$\pm$14.0	15.6$\pm$3.1	4.94$\pm$0.96	12.90$-$27.70	21
30.357-0.837_1	II	36.8$\pm$4.1	8.3$\pm$0.9	3.95$\pm$0.45	0.69$-$5.00	50
30.357-0.837_2	II	25.7$\pm$5.9	5.8$\pm$1.3	2.73$\pm$0.82	0.12$-$5.00	206

Table 6. Table 6: N 2 H +  parameters. Col. 1: Clump name; Col. 2: N 2 H +  central velocity; Col. 3: Velocity dispersion measured as 1 / ( 8 l n 2 ) 1 / 2 × 1/(8\mathrm{ln}2)^{1/2}\times FWHM of the hyperfine fit; Col. 4: Excitation temperature; Col. 5: Total optical depth.

Clump	vLSR	$\sigma$	Tex	${\tau }_{tot}$
	(km s-1)	(km s-1)	(K)
15.631-0.377	40.2(0.1)	0.30(0.01)	3.41(0.70)	3.75(0.73)
18.787-0.286	65.7(0.1)	1.07(0.01)	5.20(0.39)	1.86(0.13)
19.281-0.387	53.6(0.1)	0.47(0.03)	3.15(0.80)	5.33(1.17)
22.53-0.192	76.3(0.1)	1.25(0.01)	5.01(0.13)	1.22(0.03)
22.756-0.284	105.1(0.1)	0.95(0.01)	4.25(0.52)	1.62(0.19)
23.271-0.263	82.3(0.1)	0.94(0.01)	8.32(2.03)	0.62(0.15)
24.013+0488	94.5(0.1)	0.91(0.02)	3.76(0.36)	3.67(0.31)
25.609+0.288	113.6(0.1)	1.05(0.01)	6.64(1.43)	0.61(0.13)
25.982-0.056	89.8(0.1)	0.69(0.01)	4.19(0.80)	1.58(0.29)
28.178-0.091	98.2(0.1)	1.07(0.02)	22.71(20.44)	0.10(0.09)
28.537-0.877	88.3(0.1)	0.78(0.01)	4.37(0.35)	2.28(0.17)
28.792+0.141	107.2(0.1)	0.99(0.03)	13.91(12.52)	0.10(0.09)
30.357-0.837	78.8(0.1)	0.57(0.13)	4.06(4.25)	0.90(0.81)
31.946+0.076	96.4(0.1)	1.19(0.01)	11.52(0.36)	0.32(0.01)
32.006-0.51	71.6(0.1)	0.32(0.02)	3.02(0.70)	12.19(2.19)
34.131+0.075	56.8(0.1)	0.74(0.06)	3.85(3.36)	1.34(0.90)

Table 7. Table 7: N 2 H +  and HCO +  column density and abundances of our clumps. Col. 1: Clump name; Col. 2: H 2 column density derived from the parameters in Table 3 ; Col. 3: N 2 H +  column density; Col. 4-5: HCO +  and HCN column densities estimated iterating RADEX as described in Section 5.2 ; Col. 6-8: N 2 H + , HCO +  and HNC abundances with respect to H 2 ; Col. 9-10: Optical depth of HCO +  and HNC spectra as obtained from the RADEX run.

Clump	N(H2)	N(N2H+)	N(HCO+)	N(HNC)	$X($N2H+)	$X($HCO+)	$X($HNC)	$\tau$(HCO+)	$\tau$(HNC)
	(1022 cm-2)	(1012 cm-2)	(1012 cm-2)	(1012 cm-2)	(10-10)	(10-10)	(10-10 )
15.631-0.377	1.58(0.26)	4.97(0.55)	2.00${}_{-0.90}{}^{}{}_{}{}^{+2.40}$	12.70${}_{-5.40}{}^{}{}_{}{}^{+13.70}$	3.14(0.62)	1.26${}_{-0.60}{}^{}{}_{}{}^{+1.53}$	8.02${}_{-3.66}{}^{}{}_{}{}^{+8.76}$	4.77${}_{-2.63}{}^{}{}_{}{}^{+10.23}$	17.29${}_{-10.17}{}^{}{}_{}{}^{+34.89}$
18.787-0.286	6.91(1.55)	15.94(0.94)	6.20${}_{-2.70}{}^{}{}_{}{}^{+7.40}$	29.10${}_{-12.70}{}^{}{}_{}{}^{+32.70}$	2.31(0.53)	0.90${}_{-0.44}{}^{}{}_{}{}^{+1.09}$	4.21${}_{-2.07}{}^{}{}_{}{}^{+4.82}$	3.03${}_{-1.70}{}^{}{}_{}{}^{+6.57}$	8.79${}_{-5.14}{}^{}{}_{}{}^{+17.90}$
19.281-0.387	2.68(0.52)	10.00(1.21)	2.30${}_{-1.00}{}^{}{}_{}{}^{+2.70}$	9.60${}_{-4.10}{}^{}{}_{}{}^{+10.40}$	3.73(0.85)	0.86${}_{-0.41}{}^{}{}_{}{}^{+1.02}$	3.58${}_{-1.68}{}^{}{}_{}{}^{+3.94}$	3.15${}_{-1.75}{}^{}{}_{}{}^{+6.86}$	8.33${}_{-4.82}{}^{}{}_{}{}^{+16.90}$
22.53-0.192	4.24(1.03)	11.55(0.23)	22.70${}_{-10.00}{}^{}{}_{}{}^{+28.20}$	84.50${}_{-37.20}{}^{}{}_{}{}^{+97.30}$	2.72(0.67)	5.35${}_{-2.70}{}^{}{}_{}{}^{+6.78}$	19.93${}_{-10.03}{}^{}{}_{}{}^{+23.46}$	8.14${}_{-4.73}{}^{}{}_{}{}^{+17.25}$	19.56${}_{-11.51}{}^{}{}_{}{}^{+40.15}$
22.756-0.1284	3.72(0.78)	9.18(0.74)	0.90${}_{-0.40}{}^{}{}_{}{}^{+0.70}$	9.40${}_{-4.10}{}^{}{}_{}{}^{+10.60}$	2.47(0.55)	0.24${}_{-0.12}{}^{}{}_{}{}^{+0.19}$	2.53${}_{-1.22}{}^{}{}_{}{}^{+2.90}$	0.62${}_{-0.37}{}^{}{}_{}{}^{+1.26}$	4.14${}_{-2.34}{}^{}{}_{}{}^{+8.72}$
23.271-0.263	3.31(0.78)	9.88(2.62)	10.00${}_{-4.50}{}^{}{}_{}{}^{+11.80}$	37.30${}_{-16.40}{}^{}{}_{}{}^{+40.00}$	2.99(1.06)	3.02${}_{-1.53}{}^{}{}_{}{}^{+3.64}$	11.28${}_{-5.62}{}^{}{}_{}{}^{+12.38}$	5.74${}_{-3.26}{}^{}{}_{}{}^{+12.04}$	13.70${}_{-8.02}{}^{}{}_{}{}^{+27.41}$
24.013+0.488	5.13(1.03)	16.81(0.93)	6.20${}_{-2.70}{}^{}{}_{}{}^{+7.40}$	26.40${}_{-11.90}{}^{}{}_{}{}^{+27.20}$	3.28(0.69)	1.21${}_{-0.58}{}^{}{}_{}{}^{+1.46}$	5.15${}_{-2.54}{}^{}{}_{}{}^{+5.41}$	3.86${}_{-3.16}{}^{}{}_{}{}^{+8.27}$	10.24${}_{-5.89}{}^{}{}_{}{}^{+20.21}$
25.609+0.228	5.66(1.29)	7.50(1.53)	5.70${}_{-2.50}{}^{}{}_{}{}^{+7.00}$	17.30${}_{-7.60}{}^{}{}_{}{}^{+20.00}$	1.33(0.41)	1.01${}_{-0.50}{}^{}{}_{}{}^{+1.26}$	3.06${}_{-1.51}{}^{}{}_{}{}^{+3.60}$	3.33${}_{-1.83}{}^{}{}_{}{}^{+7.31}$	6.65${}_{-3.81}{}^{}{}_{}{}^{+14.01}$
25.982-0.056	2.39(0.52)	6.37(0.80)	19.10${}_{-8.20}{}^{}{}_{}{}^{+23.60}$	46.40${}_{-19.10}{}^{}{}_{}{}^{+52.70}$	2.67(0.67)	8.00${}_{-3.85}{}^{}{}_{}{}^{+10.04}$	19.44${}_{-9.04}{}^{}{}_{}{}^{+22.48}$	13.08${}_{-7.78}{}^{}{}_{}{}^{+27.68}$	22.88${}_{-13.89}{}^{}{}_{}{}^{+47.06}$
28.178-0.091	4.98(1.03)	2.24(2.24)	6.30${}_{-2.80}{}^{}{}_{}{}^{+7.30}$	24.50${}_{-10.90}{}^{}{}_{}{}^{+27.30}$	0.45(0.45)	1.27${}_{-0.62}{}^{}{}_{}{}^{+1.49}$	4.92${}_{-2.42}{}^{}{}_{}{}^{+5.58}$	3.36${}_{-1.88}{}^{}{}_{}{}^{+7.21}$	8.46${}_{-4.86}{}^{}{}_{}{}^{+17.29}$
28.537-0.277	4.41(0.78)	11.04(0.60)	2.50${}_{-1.00}{}^{}{}_{}{}^{+3.00}$	20.90${}_{-9.10}{}^{}{}_{}{}^{+21.80}$	2.50(0.46)	0.57${}_{-0.25}{}^{}{}_{}{}^{+0.69}$	4.74${}_{-2.22}{}^{}{}_{}{}^{+5.01}$	2.12${}_{-1.18}{}^{}{}_{}{}^{+4.66}$	9.77${}_{-5.71}{}^{}{}_{}{}^{+19.39}$
28.792+0.141	2.08(0.52)	1.53(1.53)	16.40${}_{-7.00}{}^{}{}_{}{}^{+18.10}$	56.40${}_{-23.70}{}^{}{}_{}{}^{+61.80}$	0.74(0.74)	7.90${}_{-3.90}{}^{}{}_{}{}^{+8.94}$	27.17${}_{-13.27}{}^{}{}_{}{}^{+30.53}$	9.07${}_{-5.36}{}^{}{}_{}{}^{+18.45}$	20.06${}_{-12.07}{}^{}{}_{}{}^{+40.49}$
30.357-0.837	1.42(0.26)	2.87(1.85)	1.90${}_{-0.80}{}^{}{}_{}{}^{+2.30}$	9.00${}_{-3.60}{}^{}{}_{}{}^{+9.20}$	2.02(1.35)	1.34${}_{-0.61}{}^{}{}_{}{}^{+1.63}$	6.33${}_{-2.58}{}^{}{}_{}{}^{+6.57}$	3.00${}_{-1.70}{}^{}{}_{}{}^{+6.54}$	8.69${}_{-5.24}{}^{}{}_{}{}^{+17.33}$
31.946+0.076	3.66(0.78)	11.28(0.45)	17.30${}_{-8.00}{}^{}{}_{}{}^{+20.00}$	66.40${}_{-29.10}{}^{}{}_{}{}^{+79.10}$	3.08(0.67)	4.73${}_{-2.41}{}^{}{}_{}{}^{+5.56}$	18.15${}_{-8.84}{}^{}{}_{}{}^{+21.96}$	7.62${}_{-4.25}{}^{}{}_{}{}^{+15.89}$	18.40${}_{-10.78}{}^{}{}_{}{}^{+38.58}$
32.006-0.51	1.57(0.26)	14.78(1.49)	3.20${}_{-1.40}{}^{}{}_{}{}^{+3.40}$	10.90${}_{-4.30}{}^{}{}_{}{}^{+11.80}$	9.39(1.81)	2.03${}_{-0.95}{}^{}{}_{}{}^{+2.19}$	6.93${}_{-2.96}{}^{}{}_{}{}^{+7.58}$	6.76${}_{-3.92}{}^{}{}_{}{}^{+13.82}$	14.58${}_{-9.00}{}^{}{}_{}{}^{+29.68}$
34.131+0.075	2.73(0.52)	11.81(6.01)	1.70${}_{-0.60}{}^{}{}_{}{}^{+2.10}$	10.00${}_{-4.20}{}^{}{}_{}{}^{+11.80}$	4.33(2.35)	0.62${}_{-0.25}{}^{}{}_{}{}^{+0.78}$	3.67${}_{-1.69}{}^{}{}_{}{}^{+4.38}$	1.60${}_{-0.96}{}^{}{}_{}{}^{+3.60}$	5.71${}_{-3.35}{}^{}{}_{}{}^{+12.15}$

Table 8. Table 8: Skewness parameter δ 𝛿 \delta v evaluated in each clump according to the definition of Mardones et al . ( 1997 ) . Col. 1: clump name; Col. 2-3: skewness parameter evaluated for the HCO +  and HNC spectra respectively.

Clump	$\delta$v (HCO+)	$\delta$v (HNC)
15.631-0.837	-0.06	0.06
18.787-0.286	0.71	0.58
19.281-0.387	0.70	0.75
22.53-0.192	-0.28	-0.18
22.756-0.284	1.25	-0.50
23.271-0.263	0.28	0.31
24.013+0.488	-0.70	-0.58
25.609+0.228	-0.64	-0.64
25.982-0.056	0.48	0.36
28.178-0.091	-0.82	0.19
28.537-0.277	0.89	0.63
28.792+0.141	-0.01	0.08
30.357-0.837	-0.43	-0.32
31.946+0.076	-0.95	-0.23
32.006-0.51	1.42	1.53
34.131+0.075	1.36	-0.11

Table 9. Table 9: Parameters used to estimate infall velocities adopting the Myers et al . ( 1996 ) model. The parameters have been derived fitting two Gaussians at each HCO +  spectrum showing infall signatures with the mpfitfun routine (Markwardt, 2009 ) . The uncertainties comes from the fit for all but T dip for which we assume the same uncertainties of the corresponding T red and T blue . Col. 1: Clump name; Cols. 2 − - 3: velocities of the red- and blue-shifted peak respectively; Cols. 4 − - 5: temperatures of the red- and blue-shifted peak respectively; Col. 6: temperature of the dip between the red- and blue-shifed peaks. Col. 7: line spectrum used to fit the Gaussians.

Clump	vred	vblue	Tred	Tblue	Tdip	Line
	(km s-1)	(km s-1)	K	K	K
22.53-0.192	79.86(0.03)	75.43(0.02)	0.58(0.01)	1.54(0.03)	0.52(0.01)	HCO+
22.756-0.284	106.66(0.04)	103.89(0.02)	0.31(0.01)	0.44(0.01)	0.24(0.01)	HNC
24.013+0.488	96.50(0.01)	93.00(0.01)	0.46(0.01)	0.94(0.02)	0.31(0.01)	HCO+
25.609+0.228	114.87(0.05)	112.00(0.01)	0.46(0.01)	0.53(0.01)	0.45(0.01)	HNC
28.178-0.091	99.50(0.01)	96.12(0.02)	0.65(0.01)	0.85(0.02)	0.51(0.01)	HCO+
30.357-0.837	79.75(0.03)	78.26(0.03)	0.22(0.01)	0.30(0.01)	0.17(0.01)	HNC
31.946+0.076	98.67(0.02)	93.71(0.02)	0.85(0.02)	1.22(0.02)	0.25(0.01)	HCO+

Table 10. Table 10: Infall parameters of the seven clumps with blue-shifted spectra. Col.1: Clump name; Col. 2: Infall velocity; Col. 3: Mass accretion rate derived from the infall velocity; Col. 4: Estimated free-fall time; Col. 5: Mass accreted within 1 free-fall time; Col. 6: Presence (or absence) of a 24  μ 𝜇 \mu m counterpart.

Clump	Infall vel.	Accr. rate	tff	Mff	24 $\mu$m
	(km s-1)	(10-3 M⊙ yr-1)	(${10}^5$ yr)	(M⊙)
22.53-0.192	0.34(0.01)	2.04(0.66)	3.0(0.6)	610(237)	1
22.756-0.284	0.07(0.01)	0.25(0.08)	2.6(0.6)	65(24)	0
24.013+0.488	0.21(0.01)	1.52(0.46)	2.7(0.6)	415(152)	1
25.609+0.228	0.02(0.01)	0.17(0.06)	2.8(0.6)	50(19)	0
28.178-0.091	0.05(0.01)	0.37(0.11)	2.8(0.6)	105(39)	1
30.357-0.837	0.04(0.01)	0.07(0.02)	4.7(1.0)	31(12)	1
31.946+0.076	0.37(0.01)	1.96(0.63)	3.3(0.7)	637(247)	1

Table 11. Table 11: Average values of various parameters and the associated dispersion for the three classes of 70  μ 𝜇 \mu m quiet clumps: objects with no 24  μ 𝜇 \mu m counterparts and low skewness parameter ( N ​ 24 ​ _ ​ L 𝑁 24 _ 𝐿 N24\_L ); clumps with no 24  μ 𝜇 \mu m counterparts but significant value of the skewness parameter ( N ​ 24 ​ _ ​ S 𝑁 24 _ 𝑆 N24\_S ); clumps with 24  μ 𝜇 \mu m counterparts ( Y ​ 24 𝑌 24 Y24 ). Col. 1: clump phase; Col.2: L/M ratio; Col. 3: dust temperature; Col. 4: surface density; Col. 5: velocity dispersion; Col. 6: N 2 H +  abundance relative to the H 2 ; Col. 7: mass accretion rate.

Clump	L/M	T	$\Sigma$	$\sigma$	$X($N2H+)	$\dot{M}$
group	(L⊙/M⊙)	(K)	(g cm-2)	(km s-1)	(10-10)	(10-3 M⊙ yr-1)
$N24\_L$	0.11$\pm$0.07	10.1$\pm$0.7	0.07$\pm$0.01	0.65$\pm$0.49	1.92$\pm$1.73	$-$
$N24\_S$	0.18$\pm$0.08	11.0$\pm$1.1	0.14$\pm$0.06	0.79$\pm$0.27	3.56$\pm$2.91	0.21$\pm$0.06
$Y24$	0.20$\pm$0.07	11.6$\pm$0.9	0.15$\pm$0.07	0.91$\pm$0.29	2.66$\pm$1.04	1.19$\pm$0.92

Equations

$$S_{\nu}=\frac{\mathrm{M}\kappa_{0}}{d^{2}}\bigg{(}\frac{\nu}{\nu_{0}}\bigg{)}^{\beta}B_{\nu}(T)\Omega$$

$${\mathrm{T}}_{ex}=\frac{{\mathrm{T}}_0}{\ln ({\mathrm{A}}^{-1}+1)}\mathrm{with}$$

$$A=\frac{{\mathrm{T}}_b}{{\mathrm{T}}_0(1-e^{-\tau })}+\frac{1}{e^{{\mathrm{T}}_0/{\mathrm{T}}_{bg}}-1}$$

$$\delta \mathrm{v}=\frac{{\mathrm{v}}_{\mathrm{thick}}-{\mathrm{v}}_{\mathrm{thin}}}{\Delta {\mathrm{v}}_{\mathrm{thin}}}$$

$$\mathrm{v}_{in}=\frac{\sigma^{2}}{\mathrm{v}_{red}-\mathrm{v}_{blue}}\ \mathrm{ln}\bigg{(}\frac{1+e^{\mathrm{(T_{blue}-T_{dip})/T_{dip}}}}{1+e^{\mathrm{(T_{red}-T_{dip})/T_{dip}}}}\bigg{)}.$$

Full text

1 Introduction

Massive stars play a crucial role in the formation and gas enrichment of the hosting Galaxy, and yet the formation mechanism of these extreme objects is unclear (Beuther et al., 2007; Zinnecker & Yorke, 2007; Tan et al., 2014). The massive star formation begins in molecular clouds with sufficient density to form massive objects (Tan et al., 2014). Some of these regions have been detected in absorption against the strong 8 and 24  $\mu$m background, the so-called IRDCs (Perault et al., 1996; Carey et al., 1998). Massive stars form in the densest part of the natal molecular cloud, within condensations that are called clumps (Blitz, 1993; Zinnecker & Yorke, 2007; Tan et al., 2014), objects with size $\simeq 0.5-2$ pc (Urquhart et al., 2014; Traficante et al., 2015a). The most massive of these regions, with surface density in excess of $\Sigma=0.05$ g cm*-2* (Urquhart et al., 2014), and with no signatures of on-going star formation activity are ideal massive starless clump candidates.

These extremely young clumps may be the precursors of massive starless cores, an initial condition required in core accretion models of star formation (Tan et al., 2014). It is however not well determined if massive clumps embed massive starless cores, or if they fragment into a number of low-mass cores. In a sample of 9 high-mass infared-quiet cores in Cygnus X, Duarte-Cabral et al. (2013) found that 8 out of 9 of these cores are driving outflows, therefore must be protostellar. The remaining one has only a tentative outflow detection, and could potentially be in a prestellar phase.

Starless clump candidates are hard to find, in particular for their short lifetime, of the order t$\simeq 10^{4}$ yr (Motte et al., 2007). The identification of a statistically significant number of these candidates in the Galaxy requires unbiased surveys of the Galactic Plane at wavelengths which allows us to trace the cold dust envelopes of these star forming regions, which emit principally in the far-infrared (FIR)/sub-mm.

The ATLASGAL survey (Schuller et al., 2009) observed a wide portion of the I and IV quadrant of the Galactic plane at 870 $\mu$m and produced a survey of starless clumps in the region $10^{\circ}\leq l\leq 20^{\circ}$, $|b|\leq 1^{\circ}$ (Tackenberg et al., 2012). This survey identified 210 starless clumps, but only 14 which may form stars more massive than 20 M*⊙. The search for young massive cluster (YMC) precursors in the range $20^{\circ}\geq l\geq 280^{\circ}$, combining ATLASGAL data with methanol emission, found only 7 potential YMC candidates (Urquhart et al., 2013). In the characterization of the properties of cluster progenitors combining the MALT90 (Jackson et al., 2013) and the ATLASGAL surveys, Contreras et al. (2017) identified 24 over 1244 sources as potential starless candidates and only 1 clump with properties consistent with a YMC precursor. A recent search for young massive cluster progenitors in the Galactic center using ATLASGAL and the H2O southern Galactic plane survey (HOPS) found that all YMC candidates are already forming stars (Longmore et al., 2017). These results are in agreement with the finding of Ginsburg et al. (2012) using the Bolocam Galactic Plane survey (BGPS, Aguirre et al., 2011). These authors searched for massive clumps in the first quadrant and found that none was effectively starless. More recently Svoboda et al. (2016) identified over 2000 starless clump candidates in the entire BGPS survey and a lack of candidates with masses in excess of $10^{4}$ M⊙*.

A major contribution to the field comes from the Herschel survey of the Galactic Plane, Hi-GAL (Molinari et al., 2010). Hi-GAL observed the entire Plane in the wavelength range $70\leq\lambda\leq 500$ $\mu$m allowing a direct estimation of temperature, mass and luminosity of the clumps with known distances. In particular, since the presence of a 70 $\mu$m source is interpreted as a signpost of protostellar activity (Dunham et al., 2008), the Hi-GAL survey allows a characterization of hundreds of 70 $\mu$m quiet clumps, i.e. starless clump candidates (Veneziani et al., 2013; Elia et al., 2013; Elia & et al., 2017). A first search of starless clumps with Hi-GAL was carried out by Veneziani et al. (2013). These authors analyzed Hi-GAL science verification data taken in two $2^{\circ}\times 2^{\circ}$ wide regions centered on $l=30^{\circ}$ and $l=59^{\circ}$ and found hundreds of starless clump candidates. In a $\simeq 10^{\circ}$ wide region of the outer Galaxy Elia et al. (2013) identified 688 starless clumps, the majority of them gravitationally bound sources.

Combining the Hi-GAL data with the comprehensive catalogue of IRDCs of Peretto & Fuller (2009), Traficante et al. (2015a) have performed a survey of starless and protostellar clumps associated with IRDCs with known distances ($\simeq 3500$) in the Galactic range $15^{\circ}\leq l\leq 55^{\circ}$. These authors found 667 starless clump candidates with masses up to $10^{4}$ M*⊙, $\simeq 240$ of which with surface density $\Sigma\geq 0.05$ g cm-2*, so potentially forming massive stars.

In this paper we present a detailed study of a sample of 70 $\mu$m quiet clumps mostly extracted from the Traficante et al. (2015a) catalogue for which we made follow-up observations in the dense molecular tracers N2H*+* (1$-$0), HNC (1$-$0) and HCO*+* (1$-$0) with IRAM 30m telescope111IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).. A second paper in this series (Traficante et al., 2017, submitted; hereafter, Paper II) is dedicated to the study of the properties of the non-thermal motions of these clumps.

This paper is divided as follow: in Section 2 we present the observations of the dust continuum and the line emission; in Section 3 we describe the photometry steps that we follow to obtain the fluxes of these clumps combining the Hi-GAL, ATLASGAL and BGPS datasets. In Section 4 we analyze the spectral energy distribution (SED) of the clumps and derive the main properties of their dust emission. In the same Section we also analyze the Mid-infrared (MIR) counterparts to identify clumps with faint 24 $\mu$m emission and GLIMPSE counterparts. In Section 5 we analyze the spectra of the dense gas tracers and we derive gas column densities and abundances. In Section 6 we identify clumps with evidence of infalling motions and explore the relations between dust and dense gas tracer properties, comparing the properties of the 24 $\mu$m dark and 24 $\mu$m bright sources; In Section 7 we summarize our results.

2 Observations

A sample of 17 starless clump candidates has been selected from the Traficante et al. (2015a) as objects with $\Sigma\geq 0.05$ g cm*-2*, mass $\mathrm{M}\geq 300$ M*⊙, bolometric luminosity over envelope mass ratio $\mathrm{L/M}\leq 0.3$ and very low dust temperature (T$<$15 K), indicative of very young stage of evolution (see Section 4), no (or faint) emission at 70 $\mu$m after visual inspection of each source and no counterparts in the MSX and WISE catalogues in correspondence of the Herschel dust column density peak. In addition, we checked for different masers emission associated with these clumps, as they are an indication of on-going star formation activity. We searched in the methanol multibeam survey (MMB, Green et al., 2009) and found no Class II CH3OH masers in the sources of our sample (Breen et al., 2015); from the MMB survey we also searched for hydroxyl (OH) masers at 6035 MHz (Avison et al., 2016), a transition often associated with high-mass star forming regions, and also found no associations. We searched for CH3OH and OH masers using also the Arecibo surveys of Olmi et al. (2014), which is more sensitive than the MMB survey and it is targeted to identify weak masers associated with Hi-GAL high-mass objects. We found no CH3OH masers at distances less then 100″ from the source centroids. We found one source, 34.131+0.075, with a weak OH maser (peak emission of 20 mJy, $\simeq 3\sigma$ above the r.m.s. of the observations made with the Arecibo telescope, Olmi et al., 2014). Finally, we checked several surveys of water masers in the first quadrant (Merello & et al., 2017, and references therein) and found that only one source, 23.271-0.263, has a H2*O maser association (at $\simeq 3$″ from the source centroid), identified in the survey of Svoboda et al. (2016). Note that the source 18.787-0.286 is classified as starless in Traficante et al. (2015a) catalogue and has no maser associations, although it shows a 70 $\mu$m counterpart. The 70 $\mu$m source however is faint, with a peak emission of $\simeq 60$ mJy/pixel compared to a background of $\simeq 130$ mJy/pixel. The clump follows all the other selection criteria and has very low dust temperature (T=10.6 K, see Section 4), so we include it in the analysis. The clump embedded in the cloud SDC19.281-0.387, which follows the same criteria but it is not in the Traficante et al. (2015a) catalogue, was also included in our selection. The final sample of 18 clumps presented here contains some of the most massive 70 $\mu$m quiet clumps observed in the Galaxy.

2.1 Dust continuum datasets

The dust continuum properties of the clumps have been evaluated from the Hi-GAL fluxes at 160, 250 350 and 500 $\mu$m. We combined these data with fluxes at 870 $\mu$m from the ATLASGAL survey (Schuller et al., 2009) and at 1.1 mm from the BGPS survey (Aguirre et al., 2011).

The Hi-GAL survey (Molinari et al., 2010) observed the whole Galactic plane ($|b|\leq 1^{\circ}$, and following the Galactic warp) at wavelengths of 70, 160, 250, 350 and 500 $\mu$m using both PACS (Griffin et al., 2010) and SPIRE (Poglitsch et al., 2010) instruments in parallel mode. The nominal Hi-GAL spatial resolution is $\simeq[5,10.2,18,24,34.5]$″ at [70, 160, 250, 350, 500] $\mu$m respectively. However, due to the fast scan speed mode adopted in the parallel mode, the 70 and 160 $\mu$m beams are degraded down to $\simeq 10.2$″ and $13.5$″ respectively. The sensitivity is $\simeq[27,70,16,6,6]$ MJy/sr at [70, 160, 250, 350, 500] $\mu$m respectively (Traficante et al., 2011). The data reduction follows the standard Hi-GAL data reduction pipeline (Traficante et al., 2011), and the final maps have been corrected following the weighted-GLS procedure described in Piazzo et al. (2011). The maps have been calibrated in comparison with IRAS and Planck data as described in Bernard et al. (2010).

The ATLASGAL survey (Schuller et al., 2009) covers a total of $\simeq 420$ square degrees of the Galactic Plane in the longitude range $-80^{\circ}\leq l\leq 60^{\circ}$ and has been carried out with the LABOCA instrument installed in the APEX 12m telescope. The survey has a spatial resolution of $\simeq 19.2\arcsec$ and a sensitivity of $\simeq 70$ mJy/beam in the $|l|\leq 60^{\circ}$ longitude region (Csengeri et al., 2014).

The BGPS survey has covered $\simeq 170$ square degrees of the inner Galaxy in the range $-10.5^{\circ}\leq l\leq 90.5^{\circ}$, $|b|\leq 1^{\circ}$ and has mapped the emission at 1.1 mm with a spatial resolution of 33″ and a sensitivity of 30-100 mJy/beam (Aguirre et al., 2011; Ginsburg et al., 2013).

We made dedicated photometry measurements for each clump directly on the maps instead of using the existent catalogues in order to minimize the uncertainties arising from the combination of different surveys. The method is described in in Section 3.

2.1.1 Mid-infrared sources

We searched for sources associated with each clump in the mid-infrared (MIR) using the Spitzer surveys of the Galactic Plane at 24 $\mu$m (MIPSGAL, Carey et al., 2009), and in the range $3.6-8$ $\mu$m (GLIMPSE, Benjamin et al., 2003). The MIPSGAL sensitivity is $\simeq$2.5 mJy/beam, while the GLIMPSE sensitivity is [0.5,0.5,2.0,5.0] at [3.6,4.5,5.8,8.0] $\mu$m respectively (Carey et al., 2009, and references therein). Details of the data reduction are in Benjamin et al. (2003) and (Carey et al., 2009) for GLIMPSE and MIPSGAL respectively.

2.2 Line data

Molecular line data were acquired at the IRAM 30m telescope in June 2014 under the project 034-14. The observations have been carried out with the On the Fly observing mode to map a $2\arcmin\times 2\arcmin$ wide region which covers the entire extension of each clump. Off-positions has been chosen within 30$\arcmin$ from the source centroids and checked with single pointings to verify that they were emission-free. The EMIR receiver at 3 mm was tuned at the N2H*+* (1–0) central frequency (93.17346 GHz). This tuning includes the simultaneous observations of the HNC (1$-$0) and HCO*+* 1($-$0) emission lines. The VESPA backend was tuned at the maximum spectral resolution, 20 kHz ($\simeq 0.06$ km/s), to resolve the N2H*+* hyperfine components and covers only the N2H*+* (1$-$0) emission line. The Fast Fourier Transform Spectrometer (FTS) was tuned to cover a wider range of frequencies and was used to trace the HNC and HCO*+*emission with a spectral resolution of 50 kHz ($\simeq 0.17$ km/s). The system temperature varied in the range $92\leq\mathrm{T_{sys}}\leq 162$ K. The data have been reduced with the standard GILDAS CLASS222http://www.iram.fr/IRAMFR/GILDAS software. The average sensitivity per channel of the reduced spectra has been evaluated after smoothing the data to $\simeq 0.2$ km/s and measuring the r.m.s. in 20 emission-free channels for each source. The 1-sigma r.m.s. per $\simeq 0.2$ km/s channel varies in the range $0.13\leq\mathrm{\sigma}\leq 0.32$ K. The beam FWHM at this frequency is $\simeq 27\arcsec$.

3 Far infrared dust photometry

The far infrared (FIR) dust photometry at wavelengths $160\leq\lambda\leq 1100$ $\mu$m has been done using Hyper, an enhanced aperture photometry algorithm specifically designed for crowded regions, blended sources and multi-wavelength analysis (Traficante et al., 2015b). The photometry process is the same adopted in Traficante et al. (2015a). For each source, a 2d-Gaussian fit at 250 $\mu$m defines the clumps. The fit can vary to encompass a region of at least 1 FWHM at 250 $\mu$m ($18\arcsec$) and can be up to twice the 250 $\mu$m FWHM in each direction, to avoid the contribution from underlying filamentary structures. The FWHMs of the Gaussian fit define the aperture radius. This definition of the aperture region includes at least one 500 $\mu$m beam. The aperture region is used to estimate the flux at 160, 250, 350 and 500 $\mu$m for Hi-GAL and it is also used to estimate the flux of the ATLASGAL and BGPS counterparts directly from the maps. With this choice we estimate the flux consistently at all wavelengths.

3.1 Hi-GAL clumps

The Hi-GAL fluxes of the clumps in the Traficante et al. (2015a) catalogue have been re-evaluated with Hyper parameters tuned specifically for each clump, in order to maximize the photometry accuracy of these highly confused regions. For 6 sources, the adapted photometry coincides with the photometry of the Traficante et al. (2015a) catalogue. For 9 clumps we perform a different source deblending with respect to the source catalogue. In 5 cases we deblended more sources to account for faint sources not identified in Traficante et al. (2015a), and for the other 4 cases we did not include any companion subtraction since the regions are highly confused and the background estimation dominates the emission surrounding the clump. Due to the complexity of the local background emission associated with each source, in most cases the fit reaches the maximum allowed size (FWHM=36″) along one direction. However, the Hyper fit converges for all sources but one, 23.271-0.263. We manually forced the source aperture for 28.792+0.141 and 23.271-0.263. In 28.792+0.141 the Hyper fit converges but the region is heavily confused and we forced the aperture to be circular. In 23.271-0.263 the automatic fit did not converge and we manually force the aperture region in order to encompass at least a $\simeq 30\arcsec$ region in one direction. The fluxes we estimate differ for $\simeq 25\%$ on average with the fluxes in Traficante et al. (2015a).

The Hi-GAL fluxes have been corrected for both aperture and colour corrections as described in Traficante et al. (2015a). For the colour correction we consider a clump temperature of T=11 K, the average temperature of the clumps (see Section 4).

The coordinates and photometry for all the clumps are in Table 1.

3.2 ATLASGAL clump counterparts

We evaluated the 870 $\mu$m fluxes from the ATLASGAL calibrated map for each source. For consistency, we compare our photometry with the fluxes presented in the ATLASGAL compact sources catalogue (Csengeri et al., 2014). This catalogue contains $\simeq 10000$ centrally condensed, compact objects (Csengeri et al., 2014). Thirteen clumps of our sample have been identified in the ATLASGAL catalogue. The mean difference in the flux estimation of these sources between the Hyper integrated fluxes and the fluxes in the ATLASGAL catalogue is $\simeq 40\%$, most likely for the different approaches used to evaluate the flux. The extraction method adopted by Csengeri et al. (2014), a multi-scale wavelet filtering of the large scale structures, preserves the compact dust condensations but filters out the emission arising from scales larger than $\simeq 50$″. The aperture size chosen for our clumps is up to $72\arcsec$, so it is likely that part of the flux in the ATLASGAL catalogue has been filtered out. To test for this effect and also to check the reliability of Hyper on the ATLASGAL maps we compare the Hyper photometry with the photometry of the ATLASGAL catalogue in a random region of the Galactic Plane. We chose a 3 degree wide region, $21^{\circ}\leq l\leq 24^{\circ}$, and we measure with Hyper the integrated flux in a circular region of radius R=25″ and R=30″, similar to the median of the source size of the ATLASGAL sources (27″, Csengeri et al., 2014) and of our 18 clumps (30.8″) respectively. We identified 122 sources in common in the region $21^{\circ}\leq l\leq 24^{\circ}$. The peak fluxes, less sensitive to the chosen algorithm, are in excellent agreement between Hyper and the ATLASGAL catalogue (Figure 1, upper panel). The average flux difference is only $\simeq 2\%$. The integrated fluxes are similar using an aperture radius of 25″ and slightly different using an aperture radius of 30″ (Figure 1, central and lower panel, respectively), with an average difference of $\simeq 17\%$ and $\simeq 32\%$ respectively. These tests show that the Hyper photometry on the ATLASGAL maps is reliable and that part of the large-scale flux may be filtered out in the ATLASGAL catalogue. The source fluxes at 870 $\mu$m are in Table 1.

3.3 BGPS clump counterparts

The 1100 $\mu$m flux has been evaluated from the BGPS maps, which cover 15 out of 18 sources. Three sources are not covered by the BGPS survey (22.53-0.192, 25.609+0.228 and 25.982-0.056). We corrected the Hyper photometry for a factor 1.46, the suggested aperture correction for an aperture radius of 20 ″ (Aguirre et al., 2011), since the aperture radii vary between 20″ and 35″(see Table 1). Ten of the fifteen sources have been identified in the most recent version of the catalogue, the BGPSv2.1, which contains 8594 compact sources (Ginsburg et al., 2013). The Hyper fluxes of these 10 sources in common are in good agreement with the fluxes in the BGPS catalogue, with a mean difference of $\simeq 13\%$. The good agreement between the Hyper photometry and the BGPS catalogue has been also showed in Traficante et al. (2015b), with mean flux differences of $\simeq 9\%$. The Hyper fluxes at 1.1 mm are in Table 1.

3.4 Photometry sources of uncertainties

One of the major sources of uncertainties in the flux estimation arises from the background emission associated with each clump. The background identification and removal is particularly critical for Hi-GAL data, where the cold dust emission associated with the background structures has the peak of its emission (e.g. Peretto et al., 2010). The uncertainties in the background estimation are less important in the ATLASGAL and BGPS data, since most of the extended emission is filtered out in these ground-based experiments. Furthermore, these clumps are not isolated but found in proximity of other contaminating sources. Hyper does a 2d-Gaussian modeling of the source companions which are then subtracted prior to the flux evaluation. As indicated in Table 1, we performed the source deblending in seven clumps. Finally, another source of uncertainty arises from the comparison of surveys with different sensitivities and spatial resolutions, as showed in Table 2, which could potentially affect the background estimation and the flux estimation of the clumps. Also, in ground-based experiments part of the extended emission may be filtered-out, resulting in a potential underestimation of the source fluxes. We have assumed an error on the Hi-GAL fluxes of 20% (as in Traficante et al., 2015a) and of 40% on the ATLASGAL and BGPS fluxes to account for the uncertainties arising from the combination of different surveys, namely due to the different beam responses and filtering associated with space-based and ground-based surveys.

4 Dust properties

We adopted a single-temperature greybody model to fit the source SEDs. The model assumes that the temperature gradient across the clump is small due to the absence of significance internal sources and strong interstellar radiation fields outside the clump. However, the error associated with the mass estimation assuming a single-temperature greybody model instead of solving a complete radiative transfer model is negligible for starless clumps (Wilcock et al., 2011).

The source flux $S_{\nu}$ at frequency $\nu$ is:

[TABLE]

where M is the source mass, $d$ the heliocentric distance, $\kappa_{0}$ the dust mass absorption coefficient at reference frequency $\nu_{0}=230$ GHz, fixed to 0.5 g cm*-2* and assuming a gas/dust mass ratio of 100 (Preibisch et al., 1993). $B_{\nu}(T)$ is the blackbody value at temperature $T$ and frequency $\nu$, and $\Omega$ is the solid angle of the source. The free parameters of the fit are mass and temperature. We fixed the spectral index $\beta$ for all the sources to $\beta=2.0$, appropriate for dense, cold clumps (e.g. Sadavoy et al., 2013).

We derived the physical parameters of each clump using Hi-GAL plus ATLASGAL and (when available) BGPS fluxes. The fits to the SEDs have been performed with the mpfit routine (Markwardt, 2009). The SEDs are shown in Figure 2 and the physical parameters are summarized in Table 3.

The mean temperature of these clumps is $<\mathrm{T}>=11.2\pm 1.0$ K. This temperature is significantly lower than the average temperature of starless clump candidates ($<\mathrm{T}>\simeq 15.5$ K, Traficante et al., 2015a). Massive clumps prior to the ignition of any protostars are mostly warmed-up by the external radiation field and show a temperature gradient, from about 18-28 K at the edges towards about 8-11 K at the center, where the region is shielded by external radiation fields (Peretto et al., 2010; Wilcock et al., 2011). The adopted greybody model does not take into account this variation, however the very low temperatures we measure are compatible with a cold central core and relatively low contribution of any external radiation field.

The mass of the clumps covers the range $269\leq$M$\leq 3098$ M*⊙, with a bolometric luminosity (evaluated in the range $24\,\leq\lambda\leq 1100$ $\mu$m) $18\leq\mathrm{L}\leq 669$ L⊙. The mean mass is $\simeq 1.2\times 10^{3}$ M⊙* for a mean bolometric luminosity of $\simeq$200 L*⊙*. The average L/M ratio, a good indicator of the evolution of massive regions (Molinari et al., 2016; Cesaroni et al., 2015), is only $<\mathrm{L/M}>\simeq 0.17$, $\simeq 5$ times lower than the mean L/M of starless clump candidates (1.1, Traficante et al., 2015a). $\mathrm{L/M}<<1$ is a strong indication that the regions are still quiescent (Molinari et al., 2016). In Figure 3 we show the clump distribution in a L$-$M diagram against the sample of starless clump candidates in Traficante et al. (2015a). The green tracks are the Molinari et al. (2008) evolutionary tracks for single high-mass cores. This model describes the evolution of massive cores from the beginning of the star formation process prior to the formation of a zero-age main sequence star in the L$-$M diagram, following the McKee & Tan (2003) accretion model. The high-mass cores follow a two-phases model. During the initial accretion phase the luminosity increases sustained by the collapse and the mass slightly decreases due to outflows (the analogous of a Class 0 object in the low-mass regime). The sources follow an almost vertical path in the diagram, up to a turnover point (the analogous of a class I object) after which the sources follow an almost horizontal path. This second phase begins with the formation of an HII region, the luminosity remains roughly constant while the mass is expelled through radiation and molecular outflows. The dotted line in Figure 3 corresponds to the best fit of the analogous of Class 0 objects for massive cores (Molinari et al., 2008). This model assumes that the accretion rate onto the central star increases with time. For comparison, we also show the empirical border between Class 0 and Class I sources (discussed in, e.g. André, Ward-Thompson & Barsony, 2000; Duarte-Cabral et al., 2013), which instead considers a decreasing accretion rate (blue-dotted line). In both cases, the clumps distribution lie well below the Class 0 regime, in a region characterized by extremely young objects.

4.1 Massive stars from massive clumps

In order to investigate if our clumps are likely going to form massive stars we first explore their mass-radius relationship. An empirical high-mass star formation threshold in this diagram has been proposed by Kauffmann & Pillai (2010, KP), which identified as potentially high-mass star forming regions all the clusters with $\mathrm{M(r)}>870$ M*⊙* ($\mathrm{r/pc})^{1.33}$. Recently, Baldeschi et al. (2017) analyzed the bias in the estimation of the physical parameters of massive clumps and found a relationship $\mathrm{M(r)}>1282$ M*⊙* ($\mathrm{r/pc})^{1.42}$, which is more stringent then the KP threshold. As shown in Figure 4, following the KP criterion all but three clumps are above the threshold, 15.631-0.377, 30.357-0.837 and 32.006-0.51. Conversely, following the Baldeschi et al. (2017) criterion 7 clumps may not form high-mass stars: 15.631-0.377, 30.357-0.837 and 32.006-0.51 plus 19.281-0.387, 25.982-0.056, 28.792+0.141 and 34.131+0.075.

Another criterion to identify massive star forming regions is assuming that there is a mass surface density threshold $\Sigma_{t}$ below which clumps may not form massive stars. The value of $\Sigma_{t}$ is still debated. Tan et al. (2014) assumes $0.1\leq\Sigma_{t}\leq 1$ g cm*-2* as the range of values for the threshold, while e.g. Urquhart et al. (2014) identified $\Sigma_{t}=0.05$ g cm*-2* based on the analysis of massive clumps in ATLASGAL. In our sample there are no sources with $\Sigma\leq 0.05$ g cm*-2*, and 7 sources with $\Sigma\leq 0.1$ g cm*-2* (Table 3). In Figure 4 we also show lines of constant surface density. The clumps 15.631-0.377, 30.357-0.837 and 32.006-0.51 are the sources with the lowest mass and surface density. Four sources, 18.787-0.286, 24.013+0.488, 24.528-0.136 and 25.609+0.228 have $\Sigma\geq 0.2$ g cm*-2* and $\mathrm{M}\geq 2\times 10^{3}$ M*⊙*. These clumps are among the most massive 70 $\mu$m quiet clumps observed to date and they will potentially form stars with mass comparable with the most massive protostars observed in the Galaxy (Peretto et al., 2013; Avison et al., 2015).

Combining these two criterion the majority of these clumps will likely produce massive stars.

4.2 MIR counterparts

Although 70 $\mu$m quiet clumps are good candidates to be starless (Dunham et al., 2008; Giannini et al., 2012; Veneziani et al., 2013), it may happen that some of these clumps already embed 24 $\mu$m sources (Elia & et al., 2017), identified thanks to the better MIPSGAL sensitivity compared with Hi-GAL (Table 2). In the Gutermuth & Heyer (2015) 24 $\mu$m source catalogue, 5 clumps have indeed a 24 $\mu$m counterpart. Two of them (18.787-0.286 and 30.454-0.135) however are likely to be foreground sources, as they are the only to have a 2MASS counterpart and they are not located at the Hi-GAL column density peak.

In order to look for faint 24 $\mu$m sources not identified in the Gutermuth & Heyer (2015) catalogue, we visually inspected the MIPSGAL counterparts of each clump. This inspection reveals that in 50$\%$ of the sample at least one 24 $\mu$m counterpart is present within a 250 $\mu$m beam centered in the clump centroid. Four sources (18.787-0.286, 24.013+0.488, 30.357-0.837 and 34.131+0.075) have more than one 24 $\mu$m counterpart within the clump region. The complete list of sources identified by eye, with their positions, is in Table 4.

To investigate the properties of these counterparts we performed a dedicated photometry on the MIPSGAL maps with the Aperture Photometry Tool package333http://www.aperturephotometry.org/aptool/. We estimate the photometry in a circular region with a radius of 3.5″, which includes $\simeq 1$ MIPS 24 $\mu$m beam (5.9″) in the aperture, and applied the corresponding aperture correction of 2.78, as suggested in the MIPS instrument handbook444http://irsa.ipac.caltech.edu/data/SPITZER/docs/ mips/mipsinstrumenthandbook/50/. The background has been estimated as the median value of the emission measured in a circular annulus surrounding each source. These 24 $\mu$m counterparts are very faint, with fluxes $F_{24}$ in the range $4.45\leq F_{24}\leq 29.31$ mJy. The source photometry is in Table 4.

To further characterize and classify these sources we look for counterparts in the GLIMPSE survey (Benjamin et al., 2003). Eight sources have GLIMPSE counterparts within a radius of 3″ from the 24 $\mu$m centroid (the MIPSGAL beam), associated with 5 different clumps, showed in Table 5. One source, 19.287-0.386, has 2 GLIMPSE sources associated with the 24 $\mu$m counterpart. All but 1 source have GLIMPSE counterparts at all the four IRAC bands: 3.6, 4.5, 5.8 and 8.0 $\mu$m. 18.787-0.286_2 does not have a counterpart at 8.0 $\mu$m. We use the GLIMPSE fluxes to classify these clumps according to the prescriptions of Lada (1987) and Gutermuth et al. (2009). The Lada (1987) classification scheme is based on the slope $\alpha$ of spectral index in the IRAC bands: $0\leq\alpha\leq 3$ for Class I, $-2\leq\alpha\leq 0$ for Class II and $-3\leq\alpha\leq-2$ for Class III objects. The scheme proposed by Gutermuth et al. (2009) is a IRAC colour-colour classification (Phase 1) plus a refining using a JHKS[3.6][4.5] YSO classification (Phase 2). The two classification schemes agree well and we found 6 Class I sources and 3 Class II sources. 18.787-0.286_2 has a Class I and Class II source associated with the same 24 $\mu$m counterpart. The classification is in Table 5. At least five clumps embed Class I or Class II sources but are 70 $\mu$m quiet in the Hi-GAL maps.

The MIR fluxes can be used to estimate the properties of the central stars with the Robitaille et al. (2006) SED fitting tool. This tool computes radiative transfer models of young stellar objects in a range of masses and evolutionary stages to model a central star, an accretion disk and an envelope (Robitaille et al., 2006, 2007). The tool provides several hundreds of models, each one describing a set of physical parameters with a $\chi^{2}$ value describing the goodness of the fit. The grid of models for each source are shown in Figure 5. In order to obtain a representative value for the physical parameters, we average the fit results for all models with $\chi^{2}\leq 5$. The representative values are obtained with a weighted mean, with the weight being the inverse of the $\chi^{2}$ value, similar to the procedure adopted in Grave & Kumar (2009). If the number of models with $\chi^{2}\leq 5$ is less then 50, or the $\chi^{2}$ values were always above 5, we average the results of the best 50 models. The weighted values of the extinction, column density and mass of the central star for each source are in Table 5. All these sources are highly extincted, with AV in the range $25.7\leq\mathrm{A}_{V}\leq 93.4$ mag, corresponding to surface densities $5.8\leq\Sigma\leq 21.0$ g cm*-2* applying the conversion factor described in Bohlin, Savage & Drake (1978). The best fit models are compatible with stars of intermediate mass, in the range $2.7\leq\mathrm{M}_{*}\leq 5.5$ M*⊙*.

The absence of a 70 $\mu$m counterpart in a FIR/submm bright clump is not a good indicator that the clumps are starless. Some of these clumps embed already formed intermediate mass stars, most of which may still be in the process of accreting.

In the next Sections we analyze the data obtained from the dense gas tracers and we explore if these data may help to distinguish between clumps with and without 24 $\mu$m counterparts.

5 Line analysis

We detected N2H*+* (1$-$0), HNC (1$-$0) ad HCO*+* (1$-$0) emission in all our clumps. We excluded from the line analysis the clump 24.528-0.136 as it shows absorption features due to the contamination of the off-position and 30.454-0.135 because the spectra show at least two components along the line of sight with similar intensities but we cannot separate the two sources in the dust continuum data. The spectra of the 16 clumps are presented in Appendix B.

5.1 N2H*+*

The N2H*+* (1$-$0) spectra were fitted using the CLASS hfs task, which takes into account all the hyperfine components. Total optical depth $\tau_{tot}$, central N2H*+* velocity vLSR and velocity dispersion $\sigma$ are estimated directly from the fit, following the prescription of the CLASS manual. These are shown in Table 6. The total optical depth varies in the range $0.3\leq\tau_{tot}\leq 12.2$. In four clumps the uncertainties associated with the estimation of the optical depth, $\sigma_{\tau}$, are such that $\tau_{tot}\leq 3\sigma_{\tau}$. In two of these clumps, 28.178-0.091 and 28.792+0.141, the fit gives $\tau_{tot}=0.10$ which is the lowest allowed value permitted in the hfs routine. For these four clumps we evaluated the r.m.s. of the spectrum in CLASS before and after the subtraction of the line fit. These values differ for less than 10%, suggesting that the fits are well constrained, so we accepted the parameters estimated from the fit. The optical depth of the main component can be recovered from $\tau_{tot}$ as $\tau_{main}=r_{i}*\tau_{tot}$, with $r_{i}$ the relative intensity of the main hyperfine component (0.259). The majority of the clumps for which the fit converges are moderately optically thin with $\tau_{main}\lesssim 1$. The average value is $<\tau_{main}>=0.61$, in agreement with the average optical depth of quiescent clumps associated with IRDCs (Sanhueza et al., 2012). Two sources have $\tau_{main}\geq 1$. There is an indication that the N2H*+* $(1-0)$ optical depth decreases as the clump evolves from a quiescent to an active phase (from $\overline{\tau}_{quiescent}=0.8$ to $\overline{\tau}_{active}=0.5$, Sanhueza et al., 2012), to increase again when the protostellar core forms an active HII region identified by a bright 8 $\mu$m emission ($\overline{\tau}_{red}=1.8$, Sanhueza et al., 2012). Lower optical depths in protostellar cores with respect to prestellar cores are also observed in low-mass star forming regions with comparable column densities (Crapsi et al., 2005; Emprechtinger et al., 2009).

The excitation temperature Tex is derived assuming LTE:

[TABLE]

[TABLE]

where $\mathrm{T}_{0}=h\nu/k$, $\mathrm{T}_{bg}=2.7$ K and $\mathrm{T}_{b}=\mathrm{T}^{*}_{A}\ F_{eff}/B_{eff}$. $F_{eff}$ and $B_{eff}$ are the telescope forward and beam efficiency respectively. We assume $F_{eff}=0.98$ and $B_{eff}=0.78$ (Rygl et al., 2013).

Tex varies in the range $3.0\lesssim\mathrm{T}_{ex}\lesssim 22.7$ K, as shown in Table 6. Assuming a kinetic temperature T${}_{k}=10$ K, similar to the average dust temperature of these sources, all but three clumps (28.178-0.091, 28.792+0.141, 31.946+0.076) are not thermalised, with Tex$\lesssim$Tk. However, we assumed a filling factor of 1 for all the observations which may overestimate the region of N2H*+* emission and underestimate Tex.

The N2H*+* column density has been derived as in Caselli et al. (2002) and Pineda et al. (2013), with a dipole magnetic moment $\mu_{D}=3.4$ Debye and a rigid rotor rotational constant of $B=46.58687$ GHz (Sanhueza et al., 2012, and references therein). The N2H*+* column densities are in Table 7. The average N2H*+* column density is 9.2$\times 10^{12}$ cm*-2*, with a maximum of $\simeq 1.7\times 10^{13}$ cm*-2* in 24.013+0.488. The N2H*+* column density estimates are similar to what is observed in other star forming regions (Sanhueza et al., 2012; Rygl et al., 2013).

The N2H*+* abundance with respect to the H2, $X($N2H*+), the latter derived from the dust surface density and assuming a mean molecular weight of 2.33, is in the range $0.5\times 10^{-10}\leq X($N2H+) $\leq 9.4\times 10^{-10}$ (Table 7), with an average value of $<X$(N2H+$)>=2.9\times 10^{-10}$. These values are likely to be underestimated as we are assuming a filling factor of 1, however they are in agreement with the findings in massive dense cores/clumps (Pirogov et al., 2007; Rygl et al., 2013) and in massive clumps associated with IRDCs (Ragan et al., 2006; Sanhueza et al., 2012). The source with the highest N2H+* abundance, $X($N2H*+)$\simeq 10^{-9}$, is 32.006-0.51. This is one of the sources with the lowest mass surface density, is below the KP threshold for the formation of massive stars (see Section 4.1) and does not embed any 24 $\mu$m source. However 15.631-0.377, one of the other source below the KP threshold and with no 24 $\mu$m counterparts, has a N2H+* abundance close to the average. There is no clear indication that $X($N2H*+*) is different between sources with and without a 24 $\mu$m counterpart.

5.2 HNC and HCO*+*

The HNC and HCO*+* optical depth and gas column density cannot be directly estimated from the data as we have only observed the HNC and HCO*+* ($1-0$) line, which we expect to be optically thick within these cold, dense regions (e.g Sanhueza et al., 2012).

We have estimated the HCO*+* and HNC column densities with RADEX (van der Tak et al., 2007) in a similar fashion to Peretto et al. (2013). We have well constrained measurements of the dust column density, gas temperature (assumed equal to the dust temperature) and velocity dispersion (from N2H*+* emission) for our clumps, and the only unknown variable is the gas column density. We run RADEX iteratively assuming different values of the HCO*+* and HNC column densities until the evaluated radiation temperature matched the measured peak temperature. We consider this temperature as the temperature of the main peak in the HCO*+* and HNC spectra. RADEX allows also the estimation of the optical depth of the lines. As showed in Table 7, both the HCO*+* and HNC lines are, within the uncertainties, optically thick in all clumps. Therefore, the measured temperature is a lower limit for the temperature of the main peak in both HCO*+* and HNC spectra, which gives a lower limit to the estimated gas column density. The HCO*+* and HNC column densities with the uncertainties and the abundances relative to the dust and to the N2H*+* are in Table 7. Although with the strong caveat that this procedure gives only an estimate of the gas column density using a single optically thick line observation, the values we found are comparable to those found in massive clumps (e.g. Miettinen, 2014; Zhang et al., 2016). We do not find any significant differences in the HNC or HCO*+* abundances between clumps that host a 24 $\mu$m counterpart and clumps without any MIR counterparts. Due to the uncertainties in these measurements however, we cannot give definitive conclusions on the observed trends.

5.3 Skewness in HNC and HCO*+* spectra

Optically thick line profiles can be used to identify signatures of dynamical activity in star forming regions. The difference in velocity between the brightest peaks of an optically thick line and an optically thin line can be used to compute the the skewness parameter $\delta$v (Mardones et al., 1997), defined as

[TABLE]

where $\mathrm{v_{thick}}$ is the LSR velocity of the brightest HCO*+* or HNC peak, $\mathrm{v_{thin}}$ and $\Delta\mathrm{v_{thin}}$ are respectively the LSR velocity and FWHM of N2H*+*  assumed as an optically thin line. We estimate $\mathrm{v_{thick}}$ fitting 2 Gaussians to each HCO*+*  and HNC spectrum with the mpfitfun IDL routine (Markwardt, 2009). In clumps without a well defined double peak in the HCO*+* or HNC spectra, $\mathrm{v_{thick}}$ was estimated with a single Gaussian. The skewness parameters are in Table 8. Mardones et al. (1997) define $|\delta\mathrm{v}|>0.25$ as a significant detection of skewness.

A positive skewness parameter indicates a red asymmetry in the spectrum, which could be interpreted as signature of outflows activity (e.g. Peretto, Andre & Belloche, 2006). A negative $\delta\mathrm{v}$ conversely is indicative of blue-asymmetric spectrum, a signature of infall motions in both HCO*+* (Fuller, Williams & Sridharan, 2005) and HNC (Kirk et al., 2013) spectra. Simulations of infalling high-mass star forming regions showed however that red-asymmetric HCO*+* $(1-0)$ spectra may be observed also in absence of outflows activity (Smith et al., 2013). Outflows are more reliably traced by looking for high-velocity wing emission away from the systemic velocity of the cloud. Since this is outside the scope of this paper, for the purpose of this work, we simply interpret a significant value of $\delta\mathrm{v}$ in the HCO*+* spectra as an indication of significant dynamical activity, regardless of their origin. This is the case for all but 2 clumps. Five HNC spectra have $|\delta\mathrm{v}|\leq 0.25$, an indication that HCO*+* emission traces more dynamically active gas as further explored in Paper II. Note that the clump 25.982-0.056 has symmetric HCO*+* and HNC line profiles with peaks shifted from the N2H*+* central velocity. This clumps has a skewness parameter higher than in some asymmetric, double-peaked spectra (e.g. 28.178-0.091, see spectra in Appendix B).

In the following, we restrict the analysis to clumps with asymmetric line profiles and significant blue-shifted peaks (i.e. with $\delta\mathrm{v}\leq-0.25$), as these profiles are the more unambiguously associated with infall motions. We identify 7 clumps with significant infall signatures in either HCO*+* and HNC spectra, or both. The HCO*+* spectrum of 22.756-0.284 has the blue- and red-shifted peaks of the same intensity. The HNC spectra on the contrary has a blue-shifted peak, compatible with infall signatures. Three of these clumps (22.53-0.192, 24.013+0.488 and 30.357-0.837) embed a Class I or II source (Tables 5 and 7). It is very likely that these clumps are already collapsing and these central stars are still accreting from the surrounding clumps.

5.4 Infalling properties

The spectra of the seven clumps with blue-shifted peaks can be used to determine their infall velocities vinand mass accretion rates $\dot{\mathrm{M}}$.

We calculated the infall velocities following the “two layers” model of Myers et al. (1996). According to this model, vin is:

[TABLE]

$\mathrm{v}_{red}$ and $\mathrm{v}_{blue}$ are the velocities of the red and blue peaks respectively, and $T_{red}$ and $T_{blue}$ their main beam temperatures. $T_{dip}$ is the main beam temperature of the valley between the two peaks. The parameters for each source are in Table 9. These parameters have been obtained fitting 2 Gaussians to either the HCO*+* or HNC spectrum, as indicated in Table 9, with the mpfitfun IDL routine (Markwardt, 2009). The spectra with the Gaussian fits are in Figure 8.

The infall velocities are listed in Table 10. They are in the range $0.02\leq\mathrm{v}_{in}\leq 0.37$ km s*-1*, with an average $\mathrm{\overline{v}}_{in}$=0.16 km s*-1*. Infall velocities of massive star forming regions are in the range $0.1\leq\mathrm{v_{in,mass}}\leq 1$ km s*-1* (Fuller, Williams & Sridharan, 2005), and similar velocities have been observed in massive collapsing clouds (Kirk et al., 2013; Peretto et al., 2013).

The infall velocities allow us to evaluate the mass accretion rate $\dot{\mathrm{M}}=4\pi\mathrm{R}^{2}n_{\mathrm{H}_{2}}\mu m_{\mathrm{H}}\mathrm{v}_{in}$ (Myers et al., 1996), where $m_{\mathrm{H}}$ is the hydrogen mass and $n_{\mathrm{H}_{2}}$ the volume density obtained from the dust mass and assuming spherical clumps for simplicity. The accretion rates we obtain are $0.07\leq\dot{\mathrm{M}}\leq 2.04\times 10^{-3}$ M*⊙*/yr (see Table 10). These values are comparable with the predicted accretion rates onto massive protostellar cores (McKee & Tan, 2003) and with values observed in high-mass star forming regions (Fuller, Williams & Sridharan, 2005; Rygl et al., 2013; Peretto et al., 2013) and individual protostellar sources (Duarte-Cabral et al., 2013).

The free-fall time t${}_{ff}=(3\pi/(32\mathrm{G}n_{\mathrm{H}_{2}}))^{1/2}$, with G gravitational constant and $n_{\mathrm{H}_{2}}$ gas column density, is $2.6-4.7\times 10^{5}$ yr (see Table 10), significantly higher than the massive starless clump candidates lifetime ($\lesssim 10^{4}$ yr, Motte et al., 2007; Tackenberg et al., 2012; Svoboda et al., 2016) but consistent with the accretion timescales of Duarte-Cabral et al. (2013). Within 1 free-fall time, assuming constant accretion rate equal to the value measured today, the clumps accrete a mass $31\lesssim\mathrm{M}_{accr.}\lesssim 637$ M*⊙. A clump such as 22.53-0.192 already embeds a core with a central star of $\simeq 5$ M⊙* and has the potential to accrete mass comparable with or even higher than the most massive core in SDC335 (Peretto et al., 2013; Avison et al., 2015) within one free-fall time.

Two clumps have infall signatures but no visible 24 $\mu$m counterparts, 22.756-0.284 and 25.609+0.228. We can estimate an upper limit to the mass of an embedded object assuming a constant accretion rate over time into a single protostar equal to the actual clump accretion rate and a clump lifetime prior to the formation of a detectable protostar ($t\sim 10^{4}$ yr, the expected lifetime of infrared-quiet high-mass protostars, Motte et al., 2007). This is an upper limit as the clump accretion rate measured today may have increased since the start of the collapse, and the clump may fragment in several protostellar cores. We obtain $\textrm{M}<2.5$ M*⊙* and $\textrm{M}<1.7$ M*⊙* for 22.756-0.284 and 25.609+0.228 respectively. Also, the mass upper limits are below the range of masses estimated for stars with MIPSGAL and GLIMPSE counterparts. We cannot exclude that low-mass stars may have already formed but are deeply embedded in these clumps, and likely not yet visible at 24 $\mu$m.

6 Evolutionary indicators in 70 $\mu$m quiet clumps

In this Section we compare various evolutionary indicators (L/M ratio, dust temperature, surface density, linewidth, N2H*+* abundance and mass accretion rate) to look for differences between clumps with or without a 24 $\mu$m source. We first divide our 16 clumps with well defined gas emission spectra in two groups: clumps without 24 $\mu$m counterparts ($N24$, 8 clumps) and clumps with a 24 $\mu$m counterpart ($Y24$, 8 clumps). We further divide the first group in two sub-samples: clumps with low values of the skewness parameter (2 clumps, $N24\_L$), which may be considered as the less evolved, and clumps with a significant value of the skewness parameter ($N24\_S$). The properties for the three groups are summarized in Table 11.

• $\bullet$

L/M ratio: The first indicator, the L/M ratio, is a well-identified indicator of clumps evolution (Molinari et al., 2008, 2016; Cesaroni et al., 2015). The average values of the three groups are L/M${}_{N24\_L}=0.11\pm 0.07$, L/M${}_{N24\_S}=0.18\pm 0.08$ and L/M${}_{Y24}=0.20\pm 0.07$. The L/M ratio is very low in each group and, within the dispersion of the measurements, they all exhibit a very similar value.

• $\bullet$

Dust temperature: The dust temperature is also thought to increase as the clump evolves, and the inner cores warm-up the dust envelope. The average temperatures of the three groups are T${}_{N24\_L}=10.1\pm 0.7$, T${}_{N24\_S}=11.0\pm 1.1$ and T${}_{Y24}=11.6\pm 0.9$. As for the L/M indicator, although on average the ${N24\_L}$ clumps are slightly colder than the ${N24\_S}$ and $Y24$ clumps, within the dispersion there is no clear indication of a trend among these three groups. The observed 24 $\mu$m sources may be too young to significantly alter the properties of the surrounding dust on clump scales.

• $\bullet$

Surface density: A trend of increasing surface density from more quiescent to more evolved clumps has been observed in previous surveys of star forming clumps (e.g. Urquhart et al., 2014; Svoboda et al., 2016), although it is not well established (e.g. Rathborne et al., 2010). Combining a large sample of massive clumps in the Galaxy taken from the Hi-GAL survey, Merello & et al. (2017) showed that there is no evidence of increasing $\Sigma$ with the clumps evolution. Here, we find $\Sigma_{N24\_L}=0.07\pm 0.01$, $\Sigma_{N24\_S}=0.14\pm 0.06$ and $\Sigma_{Y24}=0.15\pm 0.07$ g cm*-2*. The ${N24\_L}$ clumps have the lowest values of surface densities on average, but there is no differences between ${N24\_L}$ and $Y24$ clumps.

• $\bullet$

Velocity dispersion: An indicator of the evolution of massive clumps and cores that can be derived from the gas properties is the expected increase of the linewidth as the region evolves (Smith et al., 2013). Average values for the three groups are $\sigma_{N24\_L}=0.65\pm 0.49$, $\sigma_{N24\_S}=0.79\pm 0.27$ and $\sigma_{Y24}=0.91\pm 0.29$. There is a slight increase going from the first to the third group, however the average values of the velocity dispersion are consistent within the dispersion.

• $\bullet$

$\Sigma$* vs. $\sigma$*: In Figure 9 we show the relation between gas velocity dispersion and mass surface density of our clumps to explore if there is an evolutionary trend. The Pearson correlation coefficient is significant, 0.71. However the clumps with a 24 $\mu$m source do not occupy a specific locus of points in this plot, supporting the hypothesis that this correlation may be due to the dynamical properties of the star forming regions, and not with the clumps evolution, as suggested in some star formation models (Ballesteros-Paredes et al., 2011).

• $\bullet$

N2H+ abundance: There is evidence that the N2H*+* abundance increases as the clump evolves (Sanhueza et al., 2012). The N2H*+* abundance for the three groups are respectively $X($N2H*+$)_{N24\_L}=1.92\pm 1.73\times 10^{-10}$, $X($N2H+$)_{N24\_S}=3.56\pm 2.91\times 10^{-10}$ and $X($N2H+*$)_{Y24}=2.66\pm 1.04\times 10^{-10}$. Again, there is a weak indication that the abundance of the ${N24\_L}$ clumps is lower than in the other two groups, but consistent within the dispersion.

• $\bullet$

Mass accretion rate: In some star formation models the accretion rate is expected to increase with time (e.g. McKee & Tan, 2003). We have this information available for only 7 clumps, 2 ${N24\_S}$ and 5 $Y24$ clumps. The mean values are $\mathrm{\dot{M}}_{N24\_S}=0.21\pm 0.06\times 10^{-3}$ and $\mathrm{\dot{M}}_{Y24}=1.19\pm 0.92\times 10^{-3}$ M*⊙* yr*-1* respectively. Despite using only 7 clumps,this is suggestive that the mass accretion rate is higher in clumps with a detectable 24 $\mu$m source than in regions with still no observable inner cores (but dynamically active at the clump scale). If the accretion rate is increasing with time, this increase must be very rapid based on these (few) points, as all the other indicators do not yet show evidences of evolution among these groups of clumps.

The best candidates to embed massive pre-stellar cores are the two ${N24\_L}$ clumps, 15.631-0-377 and 28.792-0.141. These clumps have on average slightly different values of the evolutionary indicators compared to the other two groups. However, the values are compatible among the three groups within the dispersion of the measurements.

These results suggest that 70 $\mu$m quiet clumps are all at a very similar (and very early) stage of evolution. The early rise of a visible 24 $\mu$m source does not alter the properties of star forming regions at the clump scales.

7 Summary

We investigated the gas and dust properties of a sample of 18 massive clumps selected to be in a very early stage of massive star formation and 70 $\mu$m quiet.

The dust properties have been constrained combining data from the Hi-GAL, ATLASGAL and BGPS surveys. The clumps have mass of $1.2\times 10^{3}$ M*⊙* on average with 2 clumps that exceed $2\times 10^{3}$ M*⊙, and mass surface densities $\Sigma\geq 0.05$ g cm-2*. Based on the analysis of the mass surface density and the KP criterion to identify high-mass stars precursor in IRDCs, the majority of these clumps have the potential to form high-mass stars. The dust temperatures are $\mathrm{T}<13$ K, lower than the average dust temperatures of starless clump candidates ($\mathrm{T}\simeq 15$ K Traficante et al., 2015a). The luminosity is on average $\simeq 2\times 10^{2}$ L*⊙* with $\mathrm{L/M}\simeq 0.17$, significantly lower than the L/M ratio below which clumps are thought to be quiescent (L/M=1, Molinari et al., 2016). These values in clumps selected to be 70 $\mu$m quiet suggest that these massive clumps are at earliest stages of star formation.

The inspection of the 24 $\mu$m maps shows that half of these clumps have at least one faint 24 $\mu$m counterpart. Eight 24 $\mu$m sources associated with 5 different clumps have at least one GLIMPSE counterpart. We used the SED fitter tool of Robitaille et al. (2006) to get an estimate of the properties of these MIR sources and found that they all have central stars deeply embedded in the clumps with $25\lesssim\mathrm{A}_{V}\lesssim 93$ mag. These are sources with masses $2.7\lesssim\mathrm{M}_{*}\lesssim 5.5$ M*⊙* and the equivalent of low-mass Class I and Class II sources.

The gas dynamics has been studied analyzing the emission of the dense gas tracers N2H*+* ($1-0$), HNC ($1-0$) and HCO*+* ($1-0$) in the 16 clumps for which we have well defined spectral line emission. The N2H*+* emission is moderately optically thin ($<\tau>_{main}=0.6$), in line with previous observations of regions at the early stages of star formation. Blue asymmetries in HNC and HCO*+* spectra have been used to identify infall signatures. Seven clumps have blue-shifted spectra with skewness parameter $\delta_{\mathrm{v}}\leq-0.25$. Two clumps with no visible 24 $\mu$m sources have signs of infall, suggesting that they are in a dynamical state at the clump scale prior to the formation of an intermediate/high mass core.

The infall velocities are $\simeq 0.16$ km s*-1* on average, similar to what is observed in other high-mass star forming regions with hints of protostellar activity. Similarly the mass accretion rate, $0.04\leq\dot{\mathrm{M}}\leq 2.0\times 10^{-3}$ M*⊙*/yr, is comparable with other massive star forming regions. With these accretion rates, a clump such as 22.53-0.192 has the potential to form massive stars comparable with the most massive protostellar cores observed in the Galaxy to date within one free-fall time, t${}_{ff}\simeq 2.5-4.5\times 10^{5}$ yr. Assuming a lifetime of $10^{4}$ yr, clumps with infall signatures and no 24 $\mu$m sources may embed faint, low-mass protostars not detected in the MIPSGAL survey.

Finally we combine the dust properties with the gas dynamics to discuss the evolution of these clumps and to search for differences between clumps with and without 24 $\mu$m counterparts. We divided the clumps in three groups: clumps with no 24 $\mu$m counterpart and low values of the HCO*+* skewness parameter, ($N24\_L$, 2 clumps), clumps with no 24 $\mu$m counterparts but significant value of the skewness parameter ($N24\_S$, 6 clumps), and finally objects with at least one 24 $\mu$m counterpart ($Y24$, 8 clumps). We found no significant differences, within the dispersion, between these three groups from indicators as L/M ratio, dust temperature, surface density, N2H*+* velocity dispersion and gas abundance. The only evidence is that the accretion rate increases from 24 $\mu$m dark to 24 $\mu$m bright clumps. This increase of the accretion rate may be the first sign of evolution in massive clumps, as all the other indicators do not show any significant difference between clumps with and without 24 $\mu$m counterparts.

We conclude that massive starless clumps are extremely rare. The lack of 70 $\mu$m (and possibly 24 $\mu$m) emission is a necessary, but not sufficient condition to identify massive starless clumps. Massive condensations may quickly form deeply embedded protostars, and the majority, if not all of these massive clumps may already harbor low-mass fragments. High resolution observations are needed to reveal the embedded content of these high density regions.

Appendix A 24 $\mu$m and Hi-GAL images

24 $\mu$m and Higal 70-500 $\mu$m images of the 18 clumps. The blue cross in each map is the position of the 250 $\mu$m source centre.

Appendix B N2H*+* (1$-$0), HNC (1$-$0) and HCO*+* (1$-$0) spectra