New constraints on the initial parameters of low-mass star formation from chemical modeling

TL;DR

This study uses chemical modeling and observational data to identify chemical tracers that constrain the initial physical parameters of low-mass star formation, providing new insights into core properties and collapse conditions.

Contribution

It introduces a method to link chemical tracers with initial star-forming core parameters using a combination of chemical modeling and observational comparison.

Findings

Gas phase CH3CN, NS, and OCS trace initial temperature.

H2CS traces initial density and free-fall time.

Low-mass protostars form in large cores with specific size, mass, temperature, and density ranges.

Abstract

The complexity of physico-chemical models of star formation is increasing, with models that take into account new processes and more realistic setups. These models allow astrochemists to compute the evolution of chemical species throughout star formation. Hence, comparing the outputs of such models to observations allows to bring new constraints on star formation. The work presented in this paper is based on the recent public release of a database of radiation hydrodynamical low-mass star formation models. We used this database as physical parameters to compute the time dependent chemical composition of collapsing cores with a 3-phase gas-grain model. The results are analyzed to find chemical tracers of the initial physical parameters of collapse such as the mass, radius, temperature, density, and free-fall time. They are also compared to observed molecular abundances of Class 0 protostars. We find numerous tracers of the initial parameters of collapse, except for the initial mass. More particularly, we find that gas phase CH3CN, NS and OCS trace the initial temperature, while H2CS trace the initial density and free-fall time of the parent cloud. The comparison of our results with a sample of 12 Class 0 low mass protostars allows us to constrain the initial parameters of collapse of low-mass prestellar cores. We find that low-mass protostars are preferentially formed within large cores with radii greater than 20000 au, masses between 2 and 4 Msol, temperatures lower or equal to 15 K, and densities between 6e4 and 2.5e5 part.cm-3, corresponding to free-fall times between 100 and 200 kyrs.

Tables

Table 1. Table 1: Initial physical parameters ranges for the original dataset and the reference dataset for chemical modeling

IPPC	Original dataset	Reference dataset
$M_0$ (M⊙)	[0.2 ; 8]	[0.5 ; 8]
$T_0$ (K)	[5 ; 30]	[10 ; 30]
$R_0$ (au)	[3000 ; 30000]	[3000 ; 30000]
${\rho }_0$ (part.cm-3)	[2.98$\times {10}^4$ ; 1.20$\times {10}^8$]	[5.98$\times {10}^4$ ; 1.20$\times {10}^8$]
${\rho }_c$ (part.cm-3)	[4.20$\times {10}^5$ ; 1.69$\times {10}^9$]	[8.43$\times {10}^5$ ; 1.69$\times {10}^9$]
$t_{\text{ff}}$ (kyrs)	[4.31 ; 272.83]	[4.31 ; 192.92]

Table 2. Table 2: Initial abundances and cosmic-ray ionization rate ζ 𝜁 \zeta (in s -1 ). * a ​ ( b ) 𝑎 𝑏 a(b) stands for a × 10 b 𝑎 superscript 10 𝑏 a\times 10^{b} .

Element	$n_i/n_H$*	References
H2	0.5
He	0.09	1
N	6.2(-5)	2
O	2.4(-4)	3
C+	1.7(-4)	2
S+	1.5(-5)	2
Si+	8.0(-9)	4
Fe+	3.0(-9)	4
Na+	2.0(-9)	4
Mg+	7.0(-9)	4
P+	2.0(-10)	4
Cl+	1.0(-9)	4
F	6.7(-9)	5
$\zeta$	$1.3\times {10}^{-17}$

Table 3. Table 3: Summary of the calculated Spearman’s coefficients between the abundances at final time of the 482 gaseous species of the chemical network and the IPPC, using the 138 cells of the Hot Corino Region.

Initial physical parameters	${\left|{\rho }_S\right|}_{max}$	Number of species with $\left|{\rho }_S\right|>0.6$
$M_0$	0.389	0
$T_0$	0.978	137
$R_0$	0.909	88
${\rho }_0$	0.983	117
$t_{\text{ff}}$	0.983	117

Table 4. Table 5: Summary of the application our method to Schöier et al. ( 2002 ) ’s observations of IRAS 16293-2422. IC stands for inconclusive and ND means that no data were available. Upper and lower limits are only defined within their corresponding IPPC range of values (see table 1 ). In red are contradictory results, and in blue the final constraints on the IPPC of the source. e 𝑒 e is the uncertainty factor.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	15	IC	$\ge 8.96E+04$	$\le 157.52$	0.2	3	3
HCO+	$\ge 2$	$\le 20$	$\ge 20000$	$\le$2.02E+05	$\ge 105.01$	10.9	3	10
CN	IC	IC	IC	IC	IC	1.2	4	3
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge 4.98$	2.3	3	3
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge 4.98$	0.2	3	3
HC3N	IC	IC	IC	$\le$ 2.99E+07	$\ge 8.63$	0.1	3	3
CH3CN	$\ge 2$	IC	IC	$\le$ 8.96E+07	$\ge 4.98$	0.9	7	3
CCH	IC	IC	IC	IC	IC	0.3	4	3
c-C3H2	IC	IC	IC	IC	IC	5.7	6	3
CH3CCH	$\ge 1$	IC	IC	[1.20E+05 ; 8.96E+07]	[4.98 ; 136.41]	1.7	6	3
H2CO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	1.9	7	3
CH3OH	IC	IC	IC	IC	IC	1.2	23	3
CS	IC	IC	$\ge 5000$	$\le$ 3.82E+06	$\ge 24.11$	0.5	3	3
SO	$\ge 1$	$\le 25$	IC	$\le$ 8.96E+07	$\ge 4.98$	1.8	9	3
SO2	$\ge 1$	IC	IC	[1.20E+05 ; 8.96E+07]	[4.98 ; 136.41]	0.5	10	3
OCS	IC	$\le 25$	IC	IC	IC	1.5	7	3
HCS+	IC	IC	$\ge 5000$	$\le$ 3.23E+06	$\ge 26.25$	7.0	2	3
H2CS	IC	IC	IC	$\le$ 5.97E+07	$\ge 6.10$	1.4	6,4	3
SiO	$\ge 1$	IC	$\le 20000$	[3.02E+05 ; 8.96E+07]	[4.98 ; 85.74]	0.8	8	3
N2H+	$\ge 1$	IC	IC	IC	IC	ND	1	10
Constraints	[2 ; 8]	15	[5000 ; 20000]	[3.02E+05 ; 3.23E+06]	[26.25 ; 85.74]

Table 5. Table 7: Summary of the application of our method to Jørgensen et al. ( 2004 ) ’s observations of LDN1448-IRAS2.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 20	IC	IC	IC	ND	4	3
CS	IC	IC	IC	IC	IC	ND	3	10
SO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HCO+	IC	IC	IC	$\le$2.58E+07	$\ge$ 9.28	ND	1	10
N2H+	$\ge$1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HNC	IC	IC	IC	IC	IC	ND	2	3
CN	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	0.13	4	3
HC3N	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	ND	1	10
Constraints	[1 ; 8]	[10 ; 20]	[3000 ; 30000]	[8.96E+04 ; 2.58E+07]	[9.28 ; 157.52]

Table 6. Table 8: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of LDN1448-C.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	15	IC	$\ge$ 8.96E+04	$\le$ 157.52	2.1	4	3
CS	IC	IC	$\ge$ 5000	$\le$ 5.74E+06	$\ge$ 19.69	2.6	2	3
SO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HCO+	IC	IC	IC	IC	IC	0.21	4	3
N2H+	IC	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	$\ge$ 1	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	2.5	3	3
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
CN	IC	IC	IC	$\le$ 1.94E+07	$\ge$ 10.72	2.4	4	3
HC3N	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	ND	1	10
H2CO	IC	IC	$\ge$ 5000	$\le$ 7.65E+06	$\ge$ 17.05	1.2	3	3
CH3OH	IC	IC	IC	IC	IC	ND	1	10
Constraints	[1 ; 8]	15	[5000 ; 30000]	[8.96E+04 ; 5.74E+06]	[19.69 ; 157.52]

Table 7. Table 9: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of NGC1333-IRAS2A.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 15	IC	IC	IC	3.9	4	3
CS	IC	IC	$\ge$ 5000	$\le$ 3.82E+06	$\ge$ 24.11	0.91	2	3
SO	$\ge$ 1	$\le$ 25	IC	$\le$ 8.96E+07	$\ge$ 4.98	1.3	2	3
HCO+	4	10	30000	5.98E+04	192.92	0.59	4	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	7.1	3	3
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
CN	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	2.3	4	3
HC3N	IC	IC	IC	$\le$ 2.98E+07	$\ge$ 8.63	ND	1	10
H2CO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	0.63	3	3
CH3OH	IC	IC	IC	IC	IC	20.4	8	10
CH3CN	$\ge$ 2	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	8.5	4	3
Constraints	4	10	30000	5.98E+04	192.92

Table 8. Table 10: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of NGC1333-IRAS4A.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 20	IC	IC	IC	0.7	4	3
CS	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	0.019	2	3
SO	$\ge$ 1	$\le$ 25	IC	$\le$ 8.96E+07	$\ge$ 4.98	0.66	2	3
HCO+	$\ge$ 2	$\le$ 20	$\ge$ 20000	$\le$ 2.01E+05	$\ge$ 105.01	1.2	3	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	18	3	10
HNC	IC	IC	IC	IC	IC	ND	2	3
CN	IC	IC	IC	IC	IC	1.4	4	3
HC3N	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
H2CO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	0.59	3	3
CH3OH	$\ge$ 2	10	$\ge$ 5000	$\le$ 2.58E+07	$\ge$ 9.28	4.5	7	3
Constraints	[2 ; 8]	10	[20000 ; 30000]	[5.98E+04 ; 2.01E+05]	[105.01 ; 192.92]

Table 9. Table 11: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of NGC1333-IRAS4B.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 20	IC	IC	IC	0.3	4	3
CS	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	1.1	2	3
SO	$\ge$ 1	$\le$ 25	IC	$\le$ 8.96E+07	$\ge$ 4.98	0.82	2	3
HCO+	$\ge$ 2	$\le$ 20	$\ge$ 20000	$\le$ 1.20E+05	$\ge$ 136.41	2.7	2	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	1.9	3	3
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
CN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	1.5	4	3
HC3N	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
H2CO	$\ge$ 1	IC	$\ge$ 7500	$\le$ 9.56E+05	$\ge$ 48.23	6.2	3	3
CH3OH	2	10	20000	1.01E+05	148.51	11.4	7	10
Constraints	2	10	20000	1.01E+05	148.51

Table 10. Table 12: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of LDN1527.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 15	IC	$\ge$ 8.96E+04	$\le$ 157.52	1.0	4	3
CS	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	2.9	4	3
SO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	1.4	2	3
HCO+	$\ge$ 2	$\le$ 20	$\ge$ 20000	$\le$ 1.20E+05	$\ge$ 136.41	0.41	5	3
N2H+	$\ge$ 1	IC	IC	IC	IC	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
CN	IC	IC	IC	$\le$ 1.94E+07	$\ge$ 10.72	0.93	4	3
HC3N	$\ge$ 1	IC	IC	$\le$ 2.99E+07	$\ge$ 8.62	ND	1	10
H2CO	IC	IC	$\ge$ 5000	$\le$ 7.65E+06	$\ge$ 17.05	7.7	3	3
Constraints	[2 ; 8]	[10 ; 15]	[20000 ; 30000]	[8.96E+04 ; 1.20E+05]	[136.41;157.52]

Table 11. Table 13: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of VLA1623.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	15	IC	$\ge$ 8.96E+04	$\le$ 157.52	0.01	4	3
CS	IC	IC	$\ge$ 5000	$\le$ 3.82E+06	$\ge$ 24.11	2.6	2	3
SO	$\ge$ 1	$\le$ 25	$\le$ 20000	[3.02E+05 ; 8.96E+07]	[4.98 ; 85.74]	1.6	3	3
HCO+	IC	IC	IC	IC	IC	2.0	3	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
CN	IC	IC	IC	$\le$ 2.98E+07	$\ge$ 8.62	2.0	4	3
HC3N	IC	IC	IC	$\le$ 2.98E+07	$\ge$ 8.62	ND	1	10
H2CO	IC	IC	$\ge$ 5000	$\le$ 7.65E+06	$\ge$ 17.05	0.061	3	3
CH3OH	IC	IC	IC	IC	IC	ND	1	10
Constraints	[1 ; 8]	15	[5000 ; 20000]	[3.02E+05 ; 7.65E+06]	[17.05 ; 85.74]

Table 12. Table 14: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of LDN483.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 20	IC	IC	IC	1.1	4	3
CS	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	1.7	2	3
SO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HCO+	[2 ; 4]	10	$\ge$ 20000	$\le$ 1.01E+05	$\ge$ 148.51	0.63	4	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	1.0	3	3
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
CN	IC	IC	IC	$\le$ 2.99E+07	$\ge$ 8.63	2.3	4	3
HC3N	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	ND	1	10
H2CO	IC	IC	$\ge$ 5000	$\le$ 1.61E+06	$\ge$ 37.13	ND	1	10
Constraints	[2 ; 4]	10	[20000 ; 30000]	[5.98E+04 ; 1.01E+05]	[148.51 ; 192.92]

Table 13. Table 15: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of LDN723.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 20	IC	IC	IC	2.5	4	3
CS	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	13	2	10
SO	$\ge$ 1	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HCO+	4	10	30000	5.98E+04	192.92	0.61	3	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
CN	IC	IC	IC	$\le$ 2.99E+07	$\ge$ 8.63	1.8	4	3
HC3N	$\ge$ 1	IC	IC	$\le$ 2.99E+07	$\ge$ 8.63	ND	1	10
H2CO	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
Constraints	4	10	30000	5.98E+04	192.92

Table 14. Table 16: Summary of the application of our method to Jørgensen et al. ( 2004 , 2005 ) ’s observations of LDN1157.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 20	IC	IC	IC	3.3	4	3
CS	IC	IC	IC	IC	IC	ND	1	10
SO	$\ge$ 1	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	0.25	3	3
HCO+	$\ge$ 2	$\le$ 20	$\ge$ 20000	$\le$ 1.20E+05	$\ge$ 136.41	0.29	2	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HCN	IC	IC	IC	IC	IC	ND	1	10
HNC	IC	IC	IC	IC	IC	ND	2	3
CN	IC	IC	IC	IC	IC	1.1	4	3
HC3N	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	1	10
H2CO	IC	IC	IC	IC	IC	3.9	3	3
Constraints	[2 ; 8]	[10 ; 20]	[20000 ; 30000]	[8.96E+04 ; 1.20E+05]	[136.41 ; 157.52]

Table 15. Table 17: Summary of the application of our method to Jørgensen et al. ( 2004 ) ’s observations of LDN1551-IRAS5.

Species	$M_0$ ($M_{\odot }$)	$T_0$ (K)	$R_0$ (au)	${\rho }_0$ (part.cm-3)	$t_{\text{ff}}$ (kyrs)	${\chi }^2$	$N_L$	$e$
CO	IC	$\le$ 15	IC	IC	IC	1.3	4	3
CS	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	0.64	4	3
SO	$\ge$ 1	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	0.95	2	3
HCO+	[2 ; 4]	10	$\ge$ 20000	$\le$ 1.01E+05	$\ge$ 148.51	1.3	2	3
N2H+	$\ge$ 1	$\ge$ 20	IC	$\ge$ 8.96E+04	$\le$ 157.52	ND	1	10
HNC	IC	IC	IC	$\le$ 8.96E+07	$\ge$ 4.98	ND	2	3
HC3N	IC	IC	IC	$\le$ 5.97E+07	$\ge$ 6.10	ND	1	10
Constraints	[2 ; 4]	10	[20000 ; 30000]	[8.96E+04 ; 1.01E+05]	[148.51 ; 157.52]

Equations

$$\Delta R_i=(1+\alpha )\Delta R_{i+1}$$

$${\rho }_c$$

$$M_{\text{BE}}$$

$$\displaystyle\rho_{c}=5.78M_{0}\Big{(}\frac{4\pi}{3}R_{0}^{3}\Big{)}^{-1}$$

$$t_{\text{ff}}=\sqrt{\frac{3\pi }{32G{\rho }_c}}$$

$$\text{CN}+\text{CH3}\to \text{H}+\text{H2CCN}$$

$$\text{s-CN}+\text{s-CH3}\to \text{s-CH3CN}$$

$$\text{CH3CN}+\text{H3O+}\to \text{H2O}+\text{CH3CNH+}$$

$$\text{s-HCS}+\text{s-H}\to \text{s-H2CS}$$

$$\displaystyle\begin{split}\text{s-H${}_{2}$CS}+\text{s-H}&\to&\text{s-CH${}_{2}$SH}\\ &\to&\text{s-CH${}_{3}$S}\end{split}$$

$$\text{s-N}+\text{s-HS}\to \text{s-NS}+\text{s-H}$$

$$\text{b-O}+\text{b-NS}\to \text{b-NO}+\text{b-S}$$

$$\text{CH2}+\text{NS}\to \text{H}+\text{HCNS}$$

$$\text{CS}+\text{OH}\to \text{H}+\text{OCS}$$

$$\text{s-S}+\text{s-CO}\to \text{s-OCS}$$

$$\displaystyle\begin{split}\text{s-OH}+\text{s-CO}&\to&\text{s-CO${}_{2}$}+\text{s-H}\\ &\to&\text{s-HOCO}\end{split}$$

$$P(y_i)=\frac{1}{N_s}\sum _{s=1}^{N_s}{\delta }_{is}\frac{1}{N_y(s)}$$

Full text

The work presented in this paper is based on the recent public release of a database of radiation hydrodynamical low-mass star formation models. We used this database as physical parameters to compute the time dependent chemical composition of collapsing cores with a 3-phase gas-grain model. The results are analyzed to find chemical tracers of the initial physical parameters of collapse such as the mass, radius, temperature, density, and free-fall time. They are also compared to observed molecular abundances of Class 0 protostars.

We find numerous tracers of the initial parameters of collapse, except for the initial mass. More particularly, we find that gas phase CH3CN, NS and OCS trace the initial temperature, while H2CS trace the initial density and free-fall time of the parent cloud. The comparison of our results with a sample of 12 Class 0 low mass protostars allows us to constrain the initial parameters of collapse of low-mass prestellar cores. We find that low-mass protostars are preferentially formed within large cores with radii greater than 20000 au, masses between 2 and 4 M*⊙, temperatures lower or equal to 15 K, and densities between 6$\times 10^{4}$ and 2.5$\times 10^{5}$ part.cm-3*, corresponding to free-fall times between 100 and 200 kyrs.

††pubyear: 2017††pagerange: New constraints on the initial parameters of low-mass star formation from chemical modeling–A

1 Introduction

It is now commonly accepted that low-mass stars form from gravitational collapses within dense molecular clouds (see McKee & Ostriker, 2007, for a review). Indeed, the physical processes that govern the dynamic of these dense clouds, such as self-gravity, radiation field, and magnetic field, cause density fluctuations that allow the formation of denser regions within the clouds, called prestellar cores. When the thermal pressure within a prestellar core can no longer compensate its internal gravitational force, it starts collapsing onto itself. The prestellar core is initially optically thin and collapses isothermally, with compressional heating evacuated through radiation. The efficiency of this radiative cooling drops as the density of the core increases at its center. When it can no longer equilibrate with the compressive heating, the evolution of the matter at the center of the core becomes adiabatic, and a hydrostatic core, known as the first hydrostatic core (FHSC), or first Larson’s core (Larson, 1969), is formed at the center. The protostar is born. During a short period of time (< 1000 years, see for example Vaytet & Haugbolle, 2017), the FHSC will continue to accrete surrounding material from its envelope, consequently increasing its mass, density and temperature. When its temperature exceeds $\approx$2000 K, the dissociation of H2 molecules inside the core is triggered. Because of the endothermic nature of this process, a second phase of collapse begins. When most of the H2 molecules have been split, the second collapse ceases and a much more dense and compact core forms, called the second Larson’s core (Larson, 1969; Masunaga & Inutsuka, 2000; Vaytet et al., 2013). The core continues to slowly accrete material for a few hundred thousand years, until its temperature is high enough to trigger nuclear fusion reactions and the protostar becomes a young star.

In order to put constraints on the chemical evolution of the interstellar matter, from the dense molecular clouds to the multiple components of stellar systems, chemical and physical odels of prestellar collapse need to be associated. Such studies have already been conducted, from simple one-dimensional models (Ceccarelli et al., 1996; Doty et al., 2004; Aikawa et al., 2008), to multi dimensional ones (Visser et al., 2011; Furuya et al., 2012; Hincelin et al., 2016). Moreover, the complexity of these models keep on increasing. Indeed, physical models of collapse can now include several physical processes in increasingly realistic setups (see for instance Commerçon et al., 2010; Li et al., 2014; Tomida et al., 2015; Matsumoto et al., 2015). Chemical models are also constantly improved, notably via the implementation of multi-phases modeling (Taquet et al., 2012; Ruaud et al., 2016) (introducing a difference between the most external layers of the ices on the grains and the rest of the bulk ice) and new chemical processes (Garrod et al., 2007; Ruaud et al., 2015), but also thanks to public comprehensive chemical database such as UMIST (McElroy et al., 2013) and KIDA (Wakelam et al., 2015a). Such refinement of physico-chemical model allows to converge towards a more realistic representation of star formation as well as a better interpretation of the observations of star forming regions.

In this context, the work presented in this paper was initially motivated by the recent public release of a database of radiation hydrodynamical (RHD) low-mass star formation models by Vaytet & Haugbolle (2017) (http://starformation.hpc.ku.dk/grid-of-protostars). Each model simulates the 1D collapse of a spherical dark cloud described by a given set of initial physical parameters, and stops shortly after the formation of the second Larson’s core, when the collapse is believed to cease. The idea is to combine the physical models of this database with the Nautilus chemical model in order to derive useful information on the physical parameters of low-mass star formation. More particularly, it aims to challenge the possibility of putting constraints on the initial parameters of collapse of dark clouds using chemical modeling. In the first section, we detail the dataset of Vaytet & Haugbolle (2017) and the chemical model. In the second section, we detail how both were combined to obtain a set of physico-chemical models. From the results of these models, we highlight species that could be possible tracers of initial physical parameters of collapse (hereafter IPPC). Additionally, we developed a simple method to find constraints on the IPPC of observed low-mass Class 0 protostars. These results are presented respectively in the third and fourth sections. We will then discuss these results to finally conclude.

2 The physical and chemical models

In order to compute the time dependent chemical composition of collapsing cores, we have used the physical structures computed by a 1D RHD model by Vaytet & Haugbolle (2017). This model computes the collapse of a spherical cloud in a Lagrangian approach. The radius of the sphere is split in independent cells of materials that will move towards the center during the collapse, while their respective temperatures and densities will increase. Hence, for each cell and at each time step of the model, its radius, temperature, and density are known. These parameters are then used as inputs for the chemical model. In the following are described the RHD models database as well as the chemical model we used in this paper.

2.1 The RHD models database

The dataset of 1D low-mass star formation models from Vaytet & Haugbolle (2017) is composed of 143 models that follow the formation of a protostar, from the initial isothermal collapse to the FHSC, and then to the second collapse and the formation of the second Larson’s core. To solve the equation of RHD for each model, Vaytet & Haugbolle (2017) used a 1D fully implicit Godunov Lagrangian code with a grid comprising 4068 cells. Figure 1 displays for a given core a representation of the mesh used in the model and its evolution throughout the collapse. In the initial state ($t=0$, top panel), the mesh is superimposed on the radius of the core. Hence the outer radius of the outermost cell is equal to the radius R0 of the initial BE sphere, and the inner radius of the innermost cell is $R_{c}=0$. Moreover, in order to obtain a higher resolution towards the center of collapse C (black cercle), the cell sizes are decreased progressively with decreasing radius such as:

[TABLE]

where $\Delta R_{i}$ is the width of cell $i$ and $\alpha=8\times 10^{-4}$. During the collapse (t, middle panel), the mesh moves toward the center of the core, each of the cells enclosing the mass of the parcel of material it was initially matched with, notably by decreasing in size. Indeed, since the mass of each cell stays constant during the collapse, their respective densities can increase only if their sizes decrease. In the following, I will therefore use the term "cell" to refer to its corresponding parcel of material.

Each of the models is defined by a set of initial parameters and considers the initial parent cloud as a Bonnor-Ebert (BE) sphere (Bonnor, 1956; Ebert, 1955). For all models the value of the density contrast of the initial BE density profile is $\gamma=\rho_{c}/\rho_{0}=14.1$, which is the maximal value to obtain a marginally stable sphere, and which corresponds to a dimensionless radius $\xi=6.45$. In order to fully define the gas density of each initial BE sphere, fixed values for its temperature $T_{0}$ and radius $R_{0}$ were chosen. For these parameters, the central density and mass of the corresponding critical BE sphere, i.e. which defines a stable cloud, can be found via:

[TABLE]

where $k_{B}$ is the Boltzmann’s constant, $\mu$ the mean atomic weight (=2.31) and $m_{H}$ the hydrogen atom mass. The initial cloud mass $M_{0}$ is then chosen for the ratio $\epsilon=\frac{M_{\text{BE}}}{M_{0}}$ to have values below 1, corresponding to a cloud that would undergo gravitational collapse. The central density of the cloud can now be expressed as:

[TABLE]

And its corresponding free-fall time as:

[TABLE]

The initial parameter space (cf Fig. 1 of Vaytet & Haugbolle (2017)) was chosen for $\epsilon$ to cover a wide range of values between 0 and 1 while using a set of different initial cloud masses, temperatures and radii, respectively ranging from 0.2 to 8 M*⊙*, 5 to 30 K and 3000 to 30000 au. Such a wide parameter space allows to work with models covering a large diversity of molecular clouds in which low-mass protostars can form. The second column of table 1 (labelled "Original dataset") summarizes the range of values for each IPPC of interest in the present study. For the sake of clarity, it should be noted that we include the free-fall time in the term "IPPC", even if it is not a physical parameter.

2.2 The chemical model

The chemical model we use is the Nautilus 1.1 chemical model, which computes the evolution of chemical abundances for a given set of physical and chemical parameters. It can simulate a three-phase chemistry including gas-phase, grain-surface and grain bulk chemistry, along with the possible exchanges between the different phases (Ruaud et al., 2016). These exchanges are: the adsorption of gas-phase species onto grain surfaces, the thermal and non-thermal desorption of species from the grain surface into the gas-phase, and the surface-bulk and bulk-surface exchange of species. The chemical desorption process used in the model is the one depicted in Garrod et al. (2007) where we consider a 1% efficiency evaporation for all species. Moreover, the grain chemistry takes into account the standard direct photodissociation by photons along with the photodissociation induced by secondary UV photons (Prasad & Tarafdar, 1983). These processes are effective on the surface as well as in the bulk of the grains. The model also takes into account the newly implemented competition between reaction, diffusion, and evaporation as suggested by Chang et al. (2007) and Garrod & Pauly (2011). The diffusion energies of each species are computed as a fraction of their binding energies. We take for the surface a value of this ratio of 0.4 as suggested by experiments and theoretical work made on H (see Ruaud et al., 2016, and reference therein), CO and CO2 (see Karssemeijer & Cuppen, 2014). This value is then extrapolated to every species on the surface. For the bulk, we take a value of 0.8 (see also Ruaud et al., 2016). Given the high temperature regimes encountered in this study, we use the ad hoc formation mechanism for H2 described in Harada et al. (2010), as depicted in Vidal & Wakelam (2018).

The reference chemical network is kida.uva.2014 (see Wakelam et al., 2015a) for the gas-phase and the one described in Ruaud et al. (2016) for the grains. To this was added the sulphur network detailed in Vidal et al. (2017) (including the reactions given in Druard & Wakelam, 2012), as well as the chemical schemes for carbon chains proposed in Wakelam et al. (2015b), Loison et al. (2016), Hickson et al. (2016) and Loison et al. (2017). Note that all abundances in this paper are expressed with respect to the total H density.

3 Pre-treatment of the data

3.1 Selection of the reference dataset of physical models

For one set of initial parameters, i.e. one model, the RHD code gives the physical conditions of the 4096 cells of the Lagrangian grid at every time step. Detailed examples of the cell evolution in the density-temperature plane are presented in Figure 3 of Vaytet & Haugbolle (2017). We first select the models and their respective cells in order to satisfy the temperature validity domain of [10 K - 800 K] of the chemical network. With this first selection, 110 models out of 143 are useable for the chemical model application. Plus, for the computing time of Nautilus to be acceptable, we selected a number of 16 cells per RHD model between the outermost cell of the grid and the last inner cell which temperatures stay below 800 K, selecting one cell every $4096/16=256$ cells, including the outermost one. As such the reference dataset for chemical modeling includes 110 RHD models defined by a set of initial parameters, each of them comprising 16 cells.

In the third column of table 1 (labeled "Reference dataset") are displayed the resulting ranges of the selected IPPC. By comparing them with those of the original dataset, we can see that the selection process has only a small impact on the ranges of the IPPC. Indeed, except for $T_{0}$, the differences of range for $\rho_{0}$ and $t_{\text{ff}}$ are due to only two models that are ruled out because of their initial temperatures. Hence, the reference dataset still accounts for a significant range of possible IPPC.

3.2 Initial parameters of the chemical modeling

After selecting the reference dataset of physical models, the idea is now to run for each cell of the selected clouds a 0D dynamic chemical model while separately saving its radius at each time step.

As dark clouds are believed to evolve during a period of time which is not negligible for the chemistry before they can become gravitationally unstable, we first had to consider the chemical evolution of each cloud before its collapse to obtain chemical input data for the corresponding chemical model. To that effect, a standard 0D model of a dark cloud was run for each of the 110 selected clouds, using as temperature their respective initial temperatures $T_{0}$ and as proton total density the initial external density of the collapsing cloud ($n_{H}=\rho_{0}$). The choice to use the initial external density and not the central density can be justified by the fact that, given the selection of cells explained in the previous section, most of the selected cells have a initial density closer to the external density of their respective clouds. Other parameters of interest for the chemical modeling such as visual extinction and cosmic ionization rate were set to their commonly used values for dark cloud models, respectively 15 mag and $1.3\times 10^{-17}$ s*-1*. The set of initial abundances used for all pre-collapse dark cloud runs is summarized in table 2. The resulting models were run for 106 years, and their computed chemical compositions used as inputs for the chemical modeling of the corresponding cells of collapsing clouds.

In the following, we treat all the selected cells independently from each other. Hence, the outputs of the chemical modeling of the collapsing clouds describe for each of the $110\times 16=1760$ cells their respective chemical evolution during the collapse. Because most of the youngest low-mass stellar objects studied in the literature are Class 0 protostars, we focus the study of these outputs on the final time step of the chemical models, when both Larson’s cores are formed, and the object is therefore a Class 0 protostar.

3.3 Identification of two chemical regimes

In order to find possible tracers of the IPPC, we first studied the correlation plots between the abundance distribution of each species at the final time (when the protostar is formed) and the IPPC ($M_{0}$, $T_{0}$, $R_{0}$, $\rho_{0}$ and $t_{\text{ff}}$). We did not find any obvious correlations with this method. However, two chemical regimes seemed to appear in the abundance distributions of most of the species: a low temperature regime (T $\leq$ 100 K) where most species are depleted onto grains, and a high temperature regime (T $\geq$ 145 K), dominated by the evaporation of species from the grain mantles. In order to visualize these temperature intervals, we plot in figure 2 the abundance of water for all cells as a function of their final temperatures and radii. Cells within the first temperature interval are located in the outer part of the clouds (further out than a few tens of au), i.e. their Envelope Region, while cells within the two remaining intervals are located in the inner part of the clouds, i.e. their Hot Corino Region. we can thus define regions corresponding to each temperature range:

An Envelopes Region (ER), which corresponds to the temperature interval of 10 to 100 K and contains 1590 cells, 2. 2.

A Transitional Region (TR), which corresponds to the temperature interval of 100 to 145 K and contains 32 cells, 3. 3.

A Hot Corinos Region (HCR), which corresponds to the temperature interval above 145 K and contains 138 cells.

It should be noted that the separation of the cells regarding final temperatures studied here for water is common to most of the species of the chemical network. However, for the remaining species, the defined temperature intervals will not necessarily fit their abundance distributions. Indeed, TPD experiments on more volatile species ($E_{D}<E_{D}^{H_{2}O}$) show that they begin to desorb in the gas phase before water, but have their main desorption peak at the water desorption temperature as they co-desorb with water which is the main ices constituant (see Fayolle et al., 2011; Ruaud et al., 2016). Plus, more refractory species ($E_{D}>E_{D}^{H_{2}O}$) partly co-desorb with water, but are expected to desorb mainly at temperature higher than water (see for example Chaabouni et al., 2018). However, as the defined intervals would be consistent in first approximation with more volatile species, and because more refractory species only account for less than 15% of the total number of species in the chemical network, we choose to apply the region differentiation as defined above to all species regardless of their desorption energies. For the following correlation study, we decide to ignore the 33 points from the TR because they represent less than 2% of the total number of cells. Hence, the resulting loss of information is negligible.

4 Search for tracers of the IPPC

We study the correlations between species abundances and IPPC using the Spearman’s correlation coefficient ($\rho_{S}$). It assesses how well the relationship between two variables can be described using a monotonic function. Hence, it is more suitable to this study than the commonly used Pearson’s coefficient, which evaluates the linear relationship between two variables. Indeed, we first used the latter and found that linear relationships between the abundance of a given species and the IPPC were scarce. Finally we restrain our study to observed or observable species.

Note that in the ER, the results of correlations are inconclusive since the majority of the resulting correlation coefficients are too low (< 0.2 in absolute value) or concern unobservable species. This result is not conflicting with the use of the ER dataset to constrain the IPPC as detailed in section 5. Therefore, we only present here the correlation results for the HCR.

Table 3 displays maximum values of the Spearman’s coefficient obtained for each IPPC, $|\rho_{S}|_{max}$, as well as the number of species linked with them, i.e. with $|\rho_{S}|>0.6$. In the HCR, it appears that numerous correlations can be found between species abundances and IPPC, except for M0 for which we find only low Spearman coefficients ( < 0.4 in absolute value).

Among the species linked with one of the IPPC, 11 are detected toward the well studied Class 0 protostar IRAS 16293-2422. Here we present the results of the correlation for 4 of these species, for which the abundances of the HCR cells span over more than 3 orders of magnitude: CH3CN, H2CS, NS, and OCS. For the other species for which the abundances span over a smaller range of values, the results would be difficult to compare with observations. In that case, the variation of abundances would indeed not be significant enough for the corresponding species to be clear tracers, especially given the common uncertainties on observed abundances.

Figures 3, 4, 5, and 6 display the respective abundances of the 4 selected species as a function of their respective best correlated IPPC. H2CS is best correlated with the free-fall time and the density of the collapsing clouds, when CH3CN, NS and OCS are best correlated with the initial temperature. Here we try to understand how these correlations can be explained with chemistry.

4.1 CH3CN

Figure 3 shows that the abundance of acetonitrile (CH3CN) in the HCR is positively correlated with the initial temperature of the collapsing clouds with a Spearman’s coefficient of 0.713. This means that the more abundant CH3CN is in a hot corino, the higher its parent cloud’s initial temperature. More precisely, an observer that would detect CH3CN in a hot corino with an observed abundance $[CH_{3}CN]_{\text{obs}}$ could interpret this correlation such as:

• •

If $[CH_{3}CN]_{\text{obs}}\leq 10^{-8}$, $T_{0}\leq 20$ K,

• •

If $10^{-8}\leq[CH_{3}CN]_{\text{obs}}<10^{-6}$, $T_{0}\approx 25$ K,

• •

If $[CH_{3}CN]_{\text{obs}}\geq 10^{-6}$, $T_{0}\geq 30$ K.

Indeed, in the chemical model, CH3CN is known to efficiently form on the grain surfaces during the dark cloud lifetime (see Andron et al., 2018) through the hydrogenation of adsorbed H2CCN, the latter being preferentially formed in the gas phase following:

[TABLE]

At the end of the parent dark cloud runs, CH3CN abundance reaches higher values when the temperature of the cloud is higher. This increase can go up to 2 orders of magnitude in the interval of initial temperatures considered, which would explain the shape of the abundance distribution of CH3CN seen on figure 3, and the resulting correlation. This increase in CH3CN abundance with the temperature of the dark clouds is due to the fact that higher temperatures favor the diffusion of grains surface species, which allows reactions inefficient at lower temperatures to become efficient. In the case of CH3CN, it is the reaction:

[TABLE]

that becomes efficient because of the increased diffusivity of solid CN and CH3, resulting in the increase in CH3CN abundance.

Afterwards during the collapse, CH3CN thermally desorbs in the gas phase where its abundance remains stable because its main destruction reaction is:

[TABLE]

which is not very efficient because of the low abundance of H3O+. Moreover, its product CH3CNH+ forms CH3CN again via electronic recombination.

4.2 H2CS

Figure 4 shows that the abundance of thioformaldehyde (H2CS) in the HCR is positively correlated with the free-fall time with a Spearman’s coefficient of 0.633 (left panel), as well as anti-correlated with the initial density with a Spearman’s coefficient of -0.630 (right panel). Indeed, because $t_{\text{ff}}\propto\frac{1}{\sqrt{\rho_{c}}}$ (see equation 5), parent clouds with a low initial density have long free fall time and vice-versa. Given an observed abundance $[H_{2}CS]_{\text{obs}}$ in a hot corino, an observer could interpret this correlation such as:

• •

If $[H_{2}CS]_{\text{obs}}$ is close to 10*-7*, then the free-fall time of the parent cloud of the hot corino would be rather long, around 102 kyrs. Respectively, the initial density would be around 105 part.cm*-3*,

• •

If $[H_{2}CS]_{\text{obs}}<10^{-9}$, then the free-fall time would be rather short, below 50 kyrs. Respectively, the initial density would be above 106 part.cm*-3*

We find that in parent clouds with high density, H2CS forms rapidly on the grains surface through hydrogenation of HCS:

[TABLE]

but is ultimately hydrogenated in either CH2SH or CH3S via:

[TABLE]

In low density clouds however, H2CS is formed much slower on grains surface. Hence at the end of the initial cloud runs, the differences in H2CS abundances between a slow and a fast collapsing parent cloud can reach up to 2 orders of magnitude. These differences can afterwards vary during the collapse, since fast collapsing cells will reach high temperatures faster than slow collapsing cells. Indeed, once thermally desorbed in the gas phase at high temperatures, H2CS abundance can endure significant variations (see Vidal & Wakelam, 2018), which could explain the range of abundances that corresponds to fast collapsing cells ($t_{\text{ff}}<50$ kyrs) in figure 4. For slow collapsing cells, H2CS can continue to form on grains surface at low density before its thermal desorption. At that time however, the reactants of its main formation and destruction reactions at high temperature regime, such as atomic sulphur and CH3, are scarce in the gas phase, which could explain the smaller range of abundances for slow collapsing cells ($t_{\text{ff}}>50$ kyrs).

4.3 NS

Figure 5 displays the abundance distribution in the HCR of mononitrogen monosulfide (NS) and its good anti-correlation with the initial temperature of the collapsing clouds with a Spearman’s coefficient of -0.835. The corresponding Pearson’s correlation coefficient is also low (-0.843), which implies a nearly linear relationship between NS abundance and the initial temperature. As such, this relationship defines NS as a precise tracer of the initial temperature of parent clouds. Hence, for an observed value of NS in a hot corino $[NS]_{\text{obs}}$, one could deduce the following:

• •

If $[NS]_{\text{obs}}$ is close to $10^{-6}$, $T_{0}\leq 15$ K,

• •

If $10^{-8}<[NS]_{\text{obs}}\leq 10^{-7}$, $T_{0}\approx 20$ K,

• •

If $10^{-9}<[NS]_{\text{obs}}\leq 10^{-8}$, $T_{0}\approx 25$ K,

• •

If $[NS]_{\text{obs}}<10^{-10}$, $T_{0}\geq 30$ K.

In the chemical model, NS is at first efficiently formed in the initial cloud on the grains surface from the main reservoir of sulphur HS (see Vidal et al., 2017) via:

[TABLE]

This reaction is efficient in all the initial runs, regardless of the initial temperature, and the NS formed on the surfaces tends to sink in the grains bulk. However, it appears that the higher the temperature of the cloud, the more atomic oxygen can diffuse in the grains bulk. This has for main effect to destroy NS in the bulk via:

[TABLE]

This destruction of NS is therefore more efficient at high temperatures, which results in a difference between its abundance in initial clouds at 10 K and 30 K that can go up to four orders of magnitude. Moreover, once NS thermally desorbs in the hot gas, it is consumed by the following reaction:

[TABLE]

Where the CH2 comes from gas phase photodissociation by secondary UV photons of CH4, which itself forms from the CH3 produced by the destruction of thermally desorbed methanol (CH3OH) and methyl formate (HCOOCH3). This reaction is believed to be more efficient for cells that spend more time in a high temperature regime, which explains why those coming from parent clouds with $T_{0}=30$ K have such a low NS abundances compared to those coming from parents clouds with $T_{0}=10$ K.

However, this result should be taken with caution since work is still ongoing to complete the NS chemical network, as discussed in Vidal et al. (2017).

4.4 OCS

Figure 6 displays the abundances distribution in the HCR of carbonyl sulfide (OCS), which is similar to the NS one. Indeed, as for NS, a nearly-linear relationship exists between OCS abundances in the HCR and the initial temperatures of the collapsing clouds, with a Spearman’s coefficient of -0.869. Hence, OCS appears to be a very good tracer of the initial temperature. For an observed value of OCS in a hot corino $[OCS]_{\text{obs}}$, one could deduce the following:

• •

If $[OCS]_{\text{obs}}$ is close to $10^{-6}$, $T_{0}\leq 20$ K,

• •

If $[OCS]_{\text{obs}}$ is close to $10^{-7}$, $T_{0}\approx 25$ K,

• •

If $[OCS]_{\text{obs}}<10^{-9}$, $T_{0}\geq 30$ K.

It should be noted that this characterization is less precise than the one deduced from NS correlation, however we are much more confident in the completeness of the OCS chemical network than in the NS one, and therefore more confident about OCS correlation results.

From a chemical point of view, as for the previous species, the difference in OCS abundance as a function of the initial temperature can be explained with the chemistry of the initial cloud runs. Indeed, in all of them OCS is at first mainly produced in the gas phase via:

[TABLE]

It is also produced on the grains surface via:

[TABLE]

However, the higher the temperature of the initial cloud is, the less OCS is formed on the grains surface. Indeed, at low temperature, reaction (15) stays efficient until the end of the parent cloud model because CO remains abundant on the grains surface. When the temperature rises however, the surface OH originated from the photodissociation by CR induced UV photons of water and methanol, which are both abundant on the grains surface, seems to efficiently consume the remaining CO through:

[TABLE]

This reaction is believed to become efficient at high temperatures because of the increase diffusivity of OH on grains surface. Indeed, this mechanism causes differences between OCS abundances in initial clouds at 10 K and 30 K that can go up to five orders of magnitude, which would explain the corresponding correlation plot in figure 6. After thermal desorption, gas phase abundance of OCS should not vary more than one order of magnitude until the end of the collapse, as studied in Vidal & Wakelam (2018).

5 Comparison with observations

We now present the part of the study which focuses on the ER dataset, i.e. the cells that do not exceed a temperature of 100 K at the final time of their respective models (see section 3.3). Indeed, because the correlations between species abundances and IPPC in this region are not conclusive, we decided to exploit this part of the dataset in another way. Hence, we developed a simple method to efficiently constrain the IPPC from comparison with observations of Class 0 protostar envelopes. In the following, we first describe the method and the observational dataset, then continue on the presentation of a comprehensive exemple of the application of the method to the envelope of IRAS 16293-2422. Finally, we present the summary of the method applied to a sample of 12 sources.

5.1 Presentation of the method

The set of observational data used for the comparison with observations includes 12 Class 0 protostar envelopes toward which 7 to 20 species included in the chemical network had been observed. This set of observational data comprises the single -dish observations of Schöier et al. (2002), Jørgensen et al. (2004), and Jørgensen et al. (2005) and is summarized in table 4. Note that in the case of the observations of Schöier et al. (2002), we use the abundances derived for the outer part of the envelope thanks to a jump model (third column of table 7 in Schöier et al. (2002)). In each of these papers, the ability of the model to reproduce the observational constraints is quantified using the $\chi^{2}$ statistics.

Given the aforementioned dataset of observed species and their respective abundances in a Class 0 protostellar envelope, the principle of the method is to compare for each species its observed abundance to its corresponding abundance distributions in the ER for each of the IPPC considered in this study: $M_{0}$, $T_{0}$, $R_{0}$, $\rho_{0}$ and $t_{\text{ff}}$. As showed in figure 7, the results of such comparison can then either be one of the following:

• •

The observed abundance is higher than all the abundances in the ER distribution, and no constraints can be found. We note this case TH (Too High),

• •

The observed abundance intersects its respective ER distributions such as only a part of the studied IPPC range is compatible with the observations. Hence, the resulting constraint on the IPPC values will either be a single value, an interval, an upper limit or a lower limit. It should be noted that the latter are in fact also intervals since they are only defined within the studied IPPC range of values (see Reference Dataset column in table 1).

• •

The observed abundance intersects its respective ER distributions such as all the studied IPPC values are compatible, which means that no constraints on the IPPC can be found within its considered range of values. we note this case IC (InConclusive) in the following.

We do not characterize the case when the observed abundance is lower than all the abundances in the ER distribution. Indeed, for a species to be observable in the gas phase of the envelope, its abundance should be at least higher than 10*-12*, and most of species ER distribution covers a range of abundances that nearly always goes well below this limit. Hence, for each species not in the TH case, the method returns an interval of possible values for each IPPC. The final possible intervals for the source are then obtained for each IPPC as the intersection of all the corresponding intervals obtained with each species.

Finally, because $\chi^{2}$ data on calculated abundances are not available for all species in the observational results used, we define a systematic procedure in order to ensure consistency while taking into account uncertainties on observed abundances. For each species we proceed as follows:

• •

If the corresponding $\chi^{2}$ value is available, we use an error margin of a factor of 3 on the observed abundance if $\chi^{2}<10$, and of a factor of 10 if $\chi^{2}\geq 10$,

• •

If the corresponding $\chi^{2}$ value is unavailable, we use the number of lines $N_{L}$ as a criteria for the error margin definition, and use an error margin of a factor of 3 if $N_{L}>1$, and of a factor of 10 otherwise.

In figure 8, we display a comprehensive example of the application of the method to one species. we consider here SO2 in IRAS 16293-2422, which has an observed abundance of $2.25\times 10^{-10}$, and search for constraints on the initial density of the parent cloud. An error margin of factor 3 is used because the $\chi^{2}$ value is 0.5 (calculated with $N_{L}=10$). The figure shows that the observed abundance can only be compatible with an initial density of the parent cloud between $1.20\times 10^{5}$ and $8.96\times 10^{7}$ part.cm*-3*.

5.2 Results on the envelope of IRAS 16293-2422

Table 5 shows the results of the application of the method to Schöier et al. (2002)’s observations of IRAS 16293-2422. This source is an ideal example because the high number of species considered allows to display all the possible types of constraint defined in section 5.1 are found. For example, the constraint on $T_{0}$ given by CO is a single value (15 K) compatible with the upper limits given by HCO+, SO and OCS, when the $\rho_{0}$ and $t_{\text{ff}}$ constraints given by CH3CCH, SO2, and SiO are intervals. However, as for the other sources, it appears that the main type of constraints is lower or upper limits.

The result, which is by far the most interesting, is that all species give a coherent set of constraints but one. In the present case it is HCO+ that shows contradictory constraints on $R_{0}$, $\rho_{0}$ and $t_{\text{ff}}$. The $\chi^{2}$ value derived from the observations for this molecule is however the largest one in our sample, indicated an uncertain observed abundance. Additional studies of this molecule by Quénard et al. (2018b), using a 3D modeling of HCO+ and its isotopologues emission in IRAS 16293-2422, showed that the contribution of the envelope of the HCO+ emission is negligible compared to the contribution of the outflows. Such detailed analysis was however not done for other molecules for which this could be the case. The disagreement of our model with HCO+ abundance may well reflect a higher cosmic-ray ionization rate as compared to the small values that we have used. In fact, Quénard et al. (2018b)’s best model was for a ten times higher cosmic-ray ionization rate.

It appears that the method allows to successfully give constraints on the IPPC of the molecular cloud in which IRAS 16293-2422 has formed (in blue in table 5). Indeed, the fact that 95% (19/20) of the species-specific constraints agree with each other, which give confidence in the obtained results. Moreover, the constrained intervals for the envelope mass and radius are compatible with their respective values derived in Schöier et al. (2002) of 5.4 M*⊙* and 6685 au. From an astrochemical point of view, these results, and in particular those regarding initial density and temperature, force to reconsider the value commonly used in chemical models to simulate the initial cloud in which IRAS 16293-2422 had formed. Indeed, to this time, the common set of values used were usually $\leq$ 10 K for the initial temperature and a few 104 part.cm*-3* for the initial density (see for example Majumdar et al., 2017; Quénard et al., 2018a; Vidal & Wakelam, 2018). The present results suggest however that the values astrochemist should use are slightly higher, namely an initial density within the interval [3.02$\times 10^{5}$ ; 3.23$\times 10^{6}$] part.cm*-3* and an initial temperature of 15 K. This result is in agreement with those of Jaber Al-Edhari et al. (2017), who inferred from cyanopolyynes observations and chemical modeling that IRAS 16293-2422 underwent a collapse with a $t_{\text{ff}}\leq 100$ kyrs, and those of Taquet et al. (2018) who derived upper limits on the density and temperature of the parent dark cloud of IRAS 16293-2422 from deep ALMA search for O2.

5.3 Summary of the results on the source sample

5.3.1 Summary per source

Table 6 presents a summary of the IPPC constraints we find for each source listed in table 4 using the method described above (for a more comprehensive description of the results, see appendix A). As the table shows, the main result of this study is that the method allows to derive constraints on all the IPPC for all sources, except for $R_{0}$ in the case of LDN1448-IRAS2. The latter is due to the fact that for 7 of the 9 species considered for this source we applied an uncertainty factor of 10 because of the low quality of the corresponding observations (see table 7 in appendix A), making it the least constrained source of the study.

Another interesting result is that at most one species per source gives contradictory constraints except for LDN1527 for which the percentage of agreement, i.e. the percentage of species giving results that agree with the final constraint, is 100%. Hence, the percentage of agreement for the sources considered is always higher than 85%. Furthermore, except in the case of IRAS 16293-2422 (see the previous section), the species that gives contradictory constraints for the concerned sources is always N2H+. Indeed, as can be seen in figure 9, N2H+ abundances distribution in the ER as a function of $T_{0}$ is such that either its observed abundance is too high (above a few 10*-10*) and corresponds to $T_{0}\geq 20$ K (e.g. for NGC1333-IRAS2A as displayed by the red horizontal line in figure 9), either it is too low, which is the case for LDN1527 (see the horizontal blue line in figure 9), and the result is inconclusive. In the first case, the resulting constraint on $T_{0}$ is for all concerned sources contradictory with at least one other species constraint. This issue appears in 10 of the 11 sources from Jørgensen et al. (2004). Lique et al. (2015) found that the simulated line intensities of N2H+ in dark clouds increase significantly when using the newly determined collisional rate coefficients of this species with H2 compared to the one calculated based on collision with He, the latter being used in Jørgensen et al. (2004, 2005). This study indicates that the derived N2H+ abundances might be overestimated. Another possibility, as for the underestimation of HCO+ in IRAS16293-2422, would be that our adopted cosmic-ray ionization rate is too small. Using a larger $\zeta$ in our chemical model, would produce larger N2H+ abundances at low temperature. Therefore, we decided to ignore the constraint of N2H+ on $T_{0}$ in the study.

Table 6 also displays that for NGC1333-IRAS2A, NGC1333-IRAS4B and LDN723, all IPPC are constrained to single values. These single values constraints should be taken with caution, since this kind of results is model-dependent, i.e. the direct consequences of the discretization of the parameters space in the collapse models database.

Another interesting feature of table 6 is that $M_{0}$ is the least constrained IPPC. Coupled with the results found on correlation in the HCR in section 4 (see table 3), this result hints the complexity of obtaining constraints on the initial mass of the parent cloud of forming star using chemistry. On the contrary, $T_{0}$ and $\rho_{0}$ (and consequently $t_{\text{ff}}$) appear to be easily constrainable with chemistry, which make sense since temperature and density variation throughout star formation are two of the main parameters affecting the chemistry. Finally, regarding $R_{0}$, it appears that we are able to derive relatively precise constraints on its values which could seem puzzling at first since the link between the chemistry of a source and the initial radius of its parent cloud is not easy to apprehend. However, equation (2), which links the $R_{0}$ to the central density $\rho_{c}$ explains the efficiency of the method on deriving constraints on $R_{0}$.

5.3.2 Summary per IPPC

In order to study the repartition of the source constraints for each IPPC, we plot for each of them their respective probability mass functions on the sample of sources. Hence the probability of finding a low-mass protostar with a particular value of the IPPC is calculated via:

[TABLE]

where $Y=\{y_{i}\}_{i\in[1,N_{y}^{\text{tot}}]}$ is the distribution of possible values of the IPPC considered, $N_{s}=12$ is the total number of sources considered, $\delta_{is}$ is the Dirac function assessing the match between the source $s$ and the IPPC value $y_{i}$, and $N_{y}(s)$ is the number of IPPC values constrained for the source $s$. The corresponding plots display in color the contribution of each of the sources to the mass function. Consequently, in the following, we express the results in terms of probabilities for the sake of clarity. However, it should be noted that these probabilities are defined for a relatively small sample of sources (12). Therefore, the probabilities presented hereafter are likely to change for a larger sample.

Figure 10 displays the results for $M_{0}$ which is, as stated previously, the least constrained of the IPPC considered. Nevertheless, it appears that tendencies can still be inferred from these constraints. First, the probability of finding sources with $M_{0}<1$ M*⊙* is null, and of only of 5% for a parent cloud of mass $M_{0}=1$ M*⊙. Second, parent clouds of mass $M_{0}=2$ M⊙* and $M_{0}=4$ M*⊙* corresponds to a total probability of 68%. Finally, higher mass parent clouds also agree with a significant fraction of the sources, but with a much lower probability than for the previous values, namely 26%. This repartition of the source constraints hints that low-mass stars are born mostly from clouds of masses higher than 1 M*⊙. Furthermore, it highlights that a large fraction of low-mass stars should form from clouds of masses within the interval [2 ; 4] M⊙*.

Figure 11 shows the results obtained for $R_{0}$. As discussed previously, $R_{0}$ appears relatively well constrained by the method, since only two of its possible values agree with more than 4 of the sources with probabilities higher than 36%. These values correspond to the largest parent clouds considered ($R_{0}\geq 20000$ au), which suggests that low-mass stars have a total 77% probability to be born from such large clouds. The probability for low-mass stars to form from clouds with radii within the interval [5000 ; 20000] au is rather low ($0.07\times 3=21\%$) compared to the previous values. However, IRAS 16293-2422 and VLA 1623 are constrained within this interval, although with a large uncertainty on their initial radii, which hints that these sources represent particular cases of low-mass star formation. It should finally be noted that the probability to find sources corresponding to the smallest parent cloud (R0 = 3000 au) is close to 0.

The results regarding $\rho_{0}$ and $t_{\text{ff}}$ are displayed in figure 12. As both these parameters are perfectly anti-correlated, the top figure appears as the mirror of the bottom one. One of the most interesting results regarding these IPPC is that clouds of density $>10^{7}$ part.cm*-3* (respectively of $t_{\text{ff}}<17$ kyrs) appear to be ruled out of the possible parent clouds of low-mass stars. Indeed, the probability of finding such sources is only of 0.8%. Moreover, as the initial density increases (respectively $t_{\text{ff}}$ decreases), the probabilities decrease. Another interesting result is that 8 of the 12 sources are constrained within the interval [5.97$\times 10^{4}$ ; 2.02$\times 10^{5}$] part.cm*-3* (respectively $t_{\text{ff}}\in$ [105.01 ; 192.92] kyrs) indicating that approximately 70% of low-mass stars are born in clouds with relatively low densities and high free-fall times. However, as found for R0, VLA 1623 and IRAS 16293-2422 appear to represent a particular type of low-mass protostars since they are incompatible with these small values and would preferentially have formed in cloud with higher density between $3.02\times 10^{5}$ and $3.23\times 10^{6}$ part.cm*-3* (respectively $t_{\text{ff}}\in[26.25;85.74]$ kyrs), interval which corresponds to a rather low probability compared to the previous one, namely 21%. This result hints that low-mass stars can form from clouds with densities (respectively free-fall times) within a rather large range of values.

Finally, figure 13 displays the results for $T_{0}$. This IPPC is the best constrained one, which results in a clearer repartition of the mass function. First, temperatures above 20 K appear to be ruled out of the possible values since none of the sources agreed with them. Second, all sources match with an initial temperature of either 10 or 15 K, corresponding to a total probability of 94%. This shows a common trend among low-mass stars to form from a low temperature parent cloud, and particularly clouds with T0 = 10 K. Finally, the two sources corresponding to $T_{0}=20$ K (LDN1448-IRAS2 and LDN1157) are both the least constrained sources regarding initial temperature. Since the resulting probability is only of 6%, it appears that low-mass stars have a really small probability to form in clouds with temperatures as high as 20 K.

6 Discussions

In the following, we discuss the dependence of the results on the model used. Another aspect that we discuss below is the applicability of the method used in the ER to the HCR dataset.

6.1 About the modeling bias

Since we used the combination of two (physical and chemical) models, all the aforementioned results are model dependent. On the one hand, regarding the physical model, the completeness of the chemical abundances distribution over the range of the IPPC is limited by the discretization of the parameter space, especially for $M_{0}$ and $R_{0}$. Additionally, it appears for a majority of sources that the constrained lowest values of the initial densities and temperatures, as well as the highest values of the initial radii could be set by the boundaries of the parameter space. A more complete representation of these parameters could affect the results by modifying the abundance distribution of the species considered. Future similar study would therefore benefit of such a completion of the parameter space.

On the other hand, the temperature validity domain of the chemical network, especially its lower limit of 10 K forces limitations on the study. Indeed, as colder clouds can also theoretically form low-mass protostars, the chemical evolution of these clouds (available in the original dataset of Vaytet & Haugbolle (2017)) could also fit the sources considered in the ER study and therefore enlarge the range of possible $T_{0}$ we found. To solve this issue, complementary studies on the modeling of the chemistry in such low-temperature regimes are still needed. Furthermore, as work remains to be conducted on the chemical network, especially for NS, the results on the NS correlation with $T_{0}$ are susceptible to change in the future. However the chemical network for CH3CN, H2CS and OCS have been completed recently (Vidal et al., 2017, Andron et al., submitted to MNRAS), which gives confidence in the correlation results presented for these species.

Finally, the most important modeling bias is due to the fact that we used 1D models. Indeed, compared to more comprehensive 2D or 3D collapse models, 1D models do not allow to access as much detailed information, especially regarding physical structure that form towards the center during collapse such as outflows or shocks at the centrifugal barrier (Sakai et al., 2014). However 2D and 3D models require long computation times and imply large amount of output data, which would not have been compatible with the number of models and sources used in this study.

6.2 About the applicability of the method on the HCR dataset

For the sake of completeness, we also tried to apply the method developed for the ER data on the joined TR and HCR data. Both sets were taken into account in order to ensure an optimal fit with the observed data. Indeed, the cells of the TR are located at radii smaller or of the same order of magnitude (tens of au, see figure 2) as the spatial resolution of the observations of hot corinos we used. These observations are the data from Schöier et al. (2002) on the hot corinos of the binary IRAS 16293-2422 obtained with their jump model, PILS data on the hot corino within IRAS 16293-2422B (Jørgensen et al., 2016; Coutens et al., 2016; Lykke et al., 2017; Drozdovskaya et al., 2018) and Taquet et al. (2015)’s COMs observations of the hot corinos NGC1333-IRAS2A and I4B. The number of species considered for each source was respectively 11, 9, 7 and 7.

Unfortunately for all sources, at least two species gave contradictory constraints on at least two IPPC, and we were unable to derive satisfactory constraints. This failure of the method on the TR+HCR dataset could be mostly explained by the following factors:

• •

An incompleteness of the gas phase high temperature network for some of the species considered, especially COMs. Moreover, a lack of understanding of high temperature chemical mechanism such as the formation of H2 (see Wakelam et al., 2017, for a review on the subject) could also put uncertainties on the obtained chemical abundances in this region.

• •

Uncertainties on the observed abundances, notably regarding beam dilution, since the sizes of the lobe obtained with the observations are bigger or of the same order of magnitude than the scale considered in the HCR. Moreover, at the considered scale, the envelopes and outflows could also alter the quality of the observations. Finally, it is also difficult to estimate the H2 column density as the dust becomes optically thick.

• •

The use of 1D models, since they do not allow to represent the complex physical structures that appear during the collapse at its center, such as outflows and the centrifugal barrier. Hence a lack of precision is expected regarding the modeling of the inner part of the collapsing cloud, and consequently the HCR.

Nevertheless, the fact that the application of the method to the TR+HCR dataset was inconclusive does not contradict the results obtained for the ER one. In fact, we are confident in the completeness of the Nautilus chemical network for the species considered in the temperatures and densities regimes of the ER. Finally, the spatial scales considered in this region limit the uncertainties due to beam dilution.

7 Conclusions

In this article we have presented a large scale chemical study of Class 0 protostars formation, focusing on the constraint of five initial physical parameters of collapse (IPPC) of interest: $M_{0}$, $T_{0}$, $R_{0}$, $\rho_{0}$ and $t_{\text{ff}}$. The study was based on the association of a dataset of 110 1D physical models of collapsing clouds with the Nautilus chemical model. As the initial idea of the study was to use such a large dataset to infer statistical results, we first tried to find possible tracers of the IPPC by studying their respective correlations with each of the abundance distribution of the species present in the chemical network. In order to do so, we defined two physical regions in the protostellar envelope from the water abundance distribution: the Envelope Region (ER) corresponding to the the outer part of the protostellar envelope (T < 100 K) and the Hot Corino Region (HCR) corresponding to the hot corino (T>145 K). As the correlation obtained for the ER were mostly inconclusive, we focused the study on the correlation within the HCR. Despite the fact that we found no satisfactory correlations with $M_{0}$, the study shows for the four remaining IPPC $T_{0}$, $R_{0}$, $\rho_{0}$ and $t_{\text{ff}}$ that numerous possible tracers were good candidates, such as CH3CN, H2CS, NS and OCS. H2CS was found to be a possible tracer for the free-fall time and the initial density, while the three remaining species are found to be possible tracers of the initial temperature of the parent cloud.

In order to work on the ER dataset, we also developed a simple method of comparison of the abundance distributions with the observations. From this method we were able to derive the IPPC for 12 Class 0 protostars from observations of 7 to 20 species towards their envelopes. The specific result obtained for each source could be of great interest for the future modeling of their chemistry. Looking at the results for each IPPC, we were also able to derive the following probabilities for the studied low-mass protostars:

• •

The probability for them to form from parent clouds with masses $\leq 1$ M*⊙* appears to be of only 5%, while the favored interval [2 ; 4] M*⊙* reaches a probability of 68%. The outermost interval [6 ; 8] M*⊙* corresponds to a total probability of 26% and could therefore correspond to a non-negligible number of low-mass protostars,

• •

The initial temperatures, which are the most constrained IPPC, are efficiently constrained between 10 and 15 K, with all sources agreeing with either or both values. Hence, the total probability for low-mass sources to form in clouds with such temperatures reaches 94%. This result hints that, even if it is theoretically possible, low-mass stars would tend not to form from hotter parent clouds ($T_{0}\geq 20$ K),

• •

The initial radii are also well constrained with the method because of their correlation with initial densities. The result shows that low-mass stars form preferentially in vast parent clouds of $R_{0}\geq 20000$ au, with a probability of 77%. This limit is however very model dependent since the chosen value for $R_{0}$ in the original dataset are highly spaced when higher than 7500 au,

• •

Finally, the most interesting result regarding the initial densities and the free-fall times is that low-mass stars appear to have a very low probability (0.8%) to form in parent clouds of density higher than $10^{7}$ part.cm*-3* ($t_{\text{ff}}<17$ kyrs). Moreover, 70% of low-mass stars appear to be born in clouds with densities within the interval [5.97$\times 10^{4}$ ; 2.02$\times 10^{5}$] part.cm*-3* (respectively $t_{\text{ff}}\in$ [105.01 ; 192.92] kyrs). Finally the fact that the possible values obtained for all sources span between a few 104 and 106 part.cm*-3* (20 to 192 kyrs for the free-fall time) hints that low-mass stars can form in parent clouds with various possible densities (or free-fall times).

These probabilities are based on a relatively small sample of 12 Class 0 low-mass protostars which could not be representative of all low-mass protostars. Therefore, the results of our study may be subjected to change by widening the sample of sources. Change in the result may also be expected from a wider parameter space as discussed in section 6.1. Despite those two limitations, we are confident our results highlight robust tendencies regarding the initial physical parameters of formation of low mass protostars.

Another interesting results of this study is that two sources (IRAS 16293-2422 and VLA1623) have constraints on R0, $\rho_{0}$ and $t_{\text{ff}}$ within intervals which are not the favored ones, hinting that they could be representatives of a specific type of protostars which forms from medium-sized clouds with higher densities (respectively smaller free-fall time) than the majority of low-mass protostars. In order to confirm this dual trend, a similar study with a higher number of sources should be conducted.

Acknowledgements

This work has been founded by the European Research Council (Starting Grant 3DICE, grant agreement 336474). This work was supported by the Programme National "Physique et Chimie du Milieu Interstellaire" (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES.

Appendix A Constraints on the IPPC for each sources

In this appendix, we display the detailed results of the application of the comparison to observations described in section 5.1 to all the sources considered in table 3 with the exception of IRAS 16293-2422, which results are displayed in table 5. Hence, the following tables are the same as table 5 but for respectively LDN1448-IRAS2, LDN1448-C, NGC1333-IRAS2A, NGC1333-IRAS4A, NGC1333-IRAS4B, LDN1527, VLA1623, LDN483, LDN723, LDN1157, and LDN1551-IRAS5. In these table, IC stands for InConclusive and ND means that no data were available. Upper and lower limits are only defined within their corresponding IPPC range of values (see table 1). In red are contradictory results, and in blue the final constraints on the IPPC of the source. $e$ is the uncertainty factor.