A Model for Protostellar Cluster Luminosities and the Impact on the CO–H₂ Conversion Factor

We summarize the protostellar luminosity function (PLF) formalism from Offner & McKee (2011) here for completeness and discuss our extensions to the work.

The PLF is derived by adopting an accretion model, which in turn prescribes the underlying distribution of protostellar masses assuming that the final masses of the protostars obey a specified stellar initial mass function (IMF). In this framework, the accretion rate of a particular protostar, $\dot{m}$, is solely a function of its current mass, m, and its final mass, ${m}_{f}$. The protostellar mass function (PMF) describes the distribution of current protostellar masses, i.e., the present-day PMF. McKee & Offner (2010) define the PMF as

$$\begin{eqnarray}&&{\psi }_{p}(m)={\int }_{{m}_{f,{\ell }}}^{{m}_{u}}{\psi }_{p2}(m,{m}_{f})d\mathrm{ln}{m}_{f},\end{eqnarray} \tag{ 2 }$$

where ${\psi }_{p2}(m,{m}_{f})$ is the bivariate PMF that defines the fraction of protostars in a star-forming region with current masses in the range dm and final masses in the range ${{dm}}_{f}$. The bivariate PMF is related to the bivariate number density, dN²_p, within a cluster by

$$\begin{eqnarray}&&{{dN}}_{p}^{2}={N}_{p}{\psi }_{p2}(m,{m}_{f})d\mathrm{ln}{md}\mathrm{ln}{m}_{f},\end{eqnarray} \tag{ 3 }$$

where N_p is the number of protostars in the cluster. We denote the stellar IMF as ${\rm{\Psi }}({m}_{f})$. For a steady SFR,

$$\begin{eqnarray}&&{\psi }_{p2}(m,{m}_{f})=\displaystyle \frac{m{\rm{\Psi }}({m}_{f})}{\dot{m}\langle {t}_{f}\rangle },\end{eqnarray} \tag{ 4 }$$

where ${\rm{\Psi }}({m}_{f})$ is the stellar IMF, t_f is the time it takes to form a star with mass ${m}_{f}$, and $\langle {t}_{f}\rangle$ is the average time to form a star:

$$\begin{eqnarray}&&\langle {t}_{f}\rangle ={\int }_{{m}_{l}}^{{m}_{u}}d\mathrm{ln}{m}_{f}{\rm{\Psi }}({m}_{f}){t}_{f}({m}_{f}).\end{eqnarray} \tag{ 5 }$$

Following McKee & Offner (2010), we assume that ${\rm{\Psi }}({m}_{f})$ is a Chabrier IMF (Chabrier 2005) truncated at some maximum mass, m_u.

Offner & McKee (2011) parameterize the accretion model as

$$\begin{eqnarray}&&\dot{m}={\dot{m}}_{1}{\left(\displaystyle \frac{m}{{m}_{f}}\right)}^{j}{{m}_{f}}^{{j}_{f}}{\left[1-{\delta }_{n1}{\left(\displaystyle \frac{m}{{m}_{f}}\right)}^{1-j}\right]}^{1/2},\end{eqnarray} \tag{ 6 }$$

where ${\dot{m}}_{1}$ is a constant, j and j_f are model parameters, and ${\delta }_{n1}$ is a parameter determining whether the accretion rate limits to zero at t_f ("tapered"). In this study, we consider three different accretion histories:

• 1.  
  Inside-out collapse of an isothermal sphere (IS; Shu 1977), which gives

  $$\begin{eqnarray}&&\dot{m}={\dot{m}}_{\mathrm{IS}}=1.54\times {10}^{-6}{(T/10{\rm{K}})}^{3/2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 7 }$$

  where T is the gas temperature. In this model, the accretion rate is constant for a given temperature and is independent of stellar mass.
• 2.  
  Turbulent core (TC) model (McKee & Tan 2003), in which the turbulent pressure exceeds thermal pressure. The accretion rate is

  $$\begin{eqnarray}&&{\dot{m}}_{\mathrm{TC}}=3.6\times {10}^{-5}{{\rm{\Sigma }}}_{\mathrm{cl}}^{3/4}{\left(\displaystyle \frac{m}{{m}_{f}}\right)}^{j}{m}_{f}{}^{3/4}\,{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 8 }$$

  where ${{\rm{\Sigma }}}_{\mathrm{cl}}$ is the surface mass density, given in units of ${\rm{g}}\,{\mathrm{cm}}^{-2}$, and m and m_f are defined above. ${\dot{m}}_{1}\,={\dot{m}}_{\mathrm{TC}}=3.6\times {10}^{-5}{{\rm{\Sigma }}}_{\mathrm{cl}}^{3/4}$. Following McKee & Tan (2003), we use $j=\tfrac{1}{2}$. In this model, higher-mass stars accrete at higher rates.
• 3.  
  Tapered turbulent core (TTC) model (Offner & McKee 2011),

  $$\begin{eqnarray}&&{\dot{m}}_{\mathrm{TTC}}={\dot{m}}_{\mathrm{TC}}{\left[1-{\left(\displaystyle \frac{m}{{m}_{f}}\right)}^{1-j}\right]}^{1/2}\,{M}_{\odot }\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 9 }$$

  where the parameters are taken to be the same as the turbulent core model but ${\delta }_{n1}=1$. The tapered accretion rate produces smaller luminosities in later stages of protostellar evolution.

For accretion histories formulated in this way, the formation time of an individual star is

$$\begin{eqnarray}&&{t}_{f}={t}_{f1}{m}_{f}^{1-{j}_{f}}(1+{\delta }_{n1}),\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{t}_{f1}=\displaystyle \frac{1}{(1-j){\dot{m}}_{1}}\end{eqnarray} \tag{ 11 }$$

and t_f1 is the time to form a star of $1\,{M}_{\odot }$. We discuss the impact of adopting a different tapering model in Appendix B.

The PLF, ${\psi }_{p}(L)$, is defined such that ${\psi }_{p}(L)d\mathrm{ln}L$ is the fraction of protostars within the luminosity range dL. Offner & McKee (2011) showed that the bivariate PLF is related to the bivariate PMF by

$$\begin{eqnarray}&&{\psi }_{p2}(L,{m}_{f})d\mathrm{ln}L\,d\mathrm{ln}{m}_{f}={\psi }_{p2}(m,{m}_{f})d\mathrm{ln}m\,d\mathrm{ln}{m}_{f},\end{eqnarray} \tag{ 12 }$$

such that the PLF is defined as

$$\begin{eqnarray}&&{\psi }_{p}(L)=\int d\mathrm{ln}m\cdot {\psi }_{p2}(L,m).\end{eqnarray} \tag{ 13 }$$

Offner & McKee (2011) calculate the PLF by transforming Equation (13) to

$$\begin{eqnarray}&&{\psi }_{p}(L)={\int }_{{m}_{f,l}(L)}^{{m}_{u}}d\mathrm{ln}{m}_{f}\displaystyle \frac{{\psi }_{p2}(m(L),{m}_{f})}{\left|\tfrac{\partial L}{\partial m}\right|},\end{eqnarray} \tag{ 14 }$$

where ${m}_{(f,l)}(L)$ = max(m_l, m(L)).

To calculate the luminosities, we adopt the model in Offner et al. (2009), which is based on McKee & Tan (2003). This model represents the protostellar luminosity as the sum of two parts, $L={L}_{\mathrm{acc}}+{L}_{\mathrm{int}}$, where ${L}_{\mathrm{acc}}$ is the accretion luminosity and ${L}_{\mathrm{int}}$ is the internal protostellar luminosity, including Kelvin–Helmholz contraction and nuclear burning. The total accretion luminosity is defined by

$$\begin{eqnarray}&&{L}_{\mathrm{acc}}={f}_{\mathrm{acc}}\displaystyle \frac{{Gm}\dot{m}}{r},\end{eqnarray} \tag{ 15 }$$

where ${f}_{\mathrm{acc}}$ is the efficiency at which mechanical energy is converted to radiation, G is the gravitational constant, m is the protostar mass, $\dot{m}$ is the accretion rate (given by Equation (6)), and r is the protostar radius calculated following Offner et al. (2009). Following Offner & McKee (2011), we use ${f}_{\mathrm{acc}}=0.75$. The total internal luminosity is approximated by the MS mass–luminosity relationship given in Tout et al. (1996):

$$\begin{eqnarray}&&{L}_{\mathrm{int}}=\displaystyle \frac{\alpha {M}^{5.5}+\beta {M}^{11}}{\gamma +{M}^{3}+\delta {M}^{5}+\epsilon {M}^{7}+\zeta {M}^{8}+\eta {M}^{9.5}}.\end{eqnarray} \tag{ 16 }$$

Coupling this model to our PDR calculation requires some assumption about the shape of the protostellar spectrum. We assume that each luminosity component is a blackbody as described by the Planck function, such that the luminosity in a given energy range is

$$\begin{eqnarray}&&{L}_{{\rm{\Delta }}E}={f}_{{\rm{\Delta }}E}({L}_{\mathrm{acc}})\times {L}_{\mathrm{acc}}+{f}_{{\rm{\Delta }}E}({L}_{\mathrm{int}})\times {L}_{\mathrm{int}},\end{eqnarray} \tag{ 17 }$$

where f_i(L) is the fraction of the Planck function within the given energy range of interest. The blackbody temperature is derived using the Stefan–Boltzmann law with the protostar radius and luminosity from the Offner et al. (2009) model. In the limiting case where ${\rm{\Delta }}E\to \infty$, ${L}_{{\rm{\Delta }}E}\to L$.

This model is intended to represent relatively young clusters, whose membership is dominated by protostars. Appendix C discusses the results for clusters that include a secondary population of MS stars.

The PMF describes the likelihood that a cluster contains a protostar with a specific instantaneous mass and final mass. Because lower-mass stars are much more numerous than higher-mass stars, small clusters are statistically unlikely to include any high-mass stars. Under the assumption of a perfectly sampled PMF, the number of stars in a cluster, N_p, can be related to the final mass of the highest-mass star within the cluster, m_u. McKee & Offner (2010) show that the cluster size, the highest-mass star in the cluster, and the maximum possible stellar mass, ${m}_{\max }$, are related by

$$\begin{eqnarray}&&\displaystyle \frac{1}{{N}_{p}({m}_{u})}={\int }_{{m}_{u}}^{{m}_{\max }}{\psi }_{p}(m)d\mathrm{ln}m.\end{eqnarray} \tag{ 18 }$$

From an observational standpoint, the maximum mass, ${m}_{\max }$, is highly uncertain owing to a variety of factors. Crowding in clusters and unresolved binarity make measurements of individual high-mass stars challenging (Tan et al. 2014). Furthermore, constraining ${m}_{\max }$ requires measuring the populations of very young massive clusters, which are rare and distant. This work focuses mainly on small to intermediate-sized clusters (${N}_{p}\sim 10\mbox{--}{10}^{5}$), so we adopt ${m}_{\max }=100\,{M}_{\odot }$. The total cluster mass is then ${M}_{\mathrm{cl}}={N}_{p}\times \langle m\rangle$, where $\langle m\rangle \,={\int }_{{m}_{l}}^{{m}_{\max }}d\mathrm{ln}{mm}{\rm{\Psi }}(m)$. Figure 1 shows ${N}_{p}({m}_{u})$ as a function of the highest-mass star in the cluster. We adopt a minimum mass, ${m}_{\min }=0.033\,{M}_{\odot }$. For the TTC model, the average mass $\langle m\rangle \approx 0.2\,{M}_{\odot }$.

Zoom In Zoom Out Reset image size

Equation (14) can be numerically integrated given a protostellar model for L(m) and $r(m,{m}_{f})$. This approach allows the distribution to be calculated exactly, i.e., direct integration produces perfect sampling of the underlying function. However, the stellar radius undergoes several discontinuous jumps owing to changes in the nuclear state (e.g., Figure 5 in Offner & Arce 2014) and is consequently difficult to invert. Moreover, the mass functions of small clusters are subject to Poisson statistics and, thus, not perfectly sampled. We therefore adopt a statistical approach to compute the PLF and cluster properties.

We calculate the PLF and PMF of a cluster using the conditional probability method. The first step of the method is to marginalize the bivariate PMF (4) over the protostar final mass, m_f. This one-dimensional distribution function is then sampled for a protostar mass, m, using the inversion method numerically. We then calculate the conditional probability distribution for the final mass given the current mass, $\psi ({m}_{f}| m)=({\psi }_{{\rm{p}}2}(m=m,{m}_{f}))/({\psi }_{p}(m=m))$.The conditional probability is then sampled using the inversion method again to obtain the final mass, m_f. This procedure is done for as many protostars as in each cluster. The protostellar masses drawn this way converge to the analytic PMF with a sample of 10⁵ protostars. Figure 2 shows the convergence of the PMF distribution to the analytic result for the isothermal sphere accretion model as a function of the number of stars included in the distribution. We find that the distribution converges well to the analytic distribution by ${N}_{* }\approx {10}^{5}$.

Zoom In Zoom Out Reset image size

To calculate cluster statistics, we draw N_* protostars for a number of mock clusters, ${N}_{\mathrm{cl}}$, using the procedure described above. For each mock cluster, we calculate the bolometric, FUV, and ionizing luminosities for each protostar using Equation (17). The total luminosities and masses are calculated for the mock cluster. After drawing ${N}_{\mathrm{cl}}$ mock clusters, we calculate the average and spread of the different total luminosities and the mass. When we compare to observations in Section 3, we use ${N}_{\mathrm{cl}}=100$ to achieve statistical robustness for the mean and the spread. For the chemistry, we use the average of the total cluster luminosities and masses. As such, we optimize the procedure by calculating the running mean of the total bolometric luminosity and drawing clusters until the running mean converges to 0.1% relative error. We find that the running mean converges in ${N}_{{cl}}\approx 15\mbox{--}20$ across 4 dex of N_*.

We use the PDR code 3d-pdr³ (Bisbas et al. 2012) to model the chemistry of the molecular gas in our models. 3d-pdr obtains the gas temperature and abundance distributions for a given input density distribution by balancing the heating and cooling. Cooling mainly occurs owing to [C i], [O i], and [C ii] forbidden line emission. 3d-pdr includes four heating mechanisms: (i) photoelectric heating of dust grains due to FUV radiation, (ii) de-excitation of vibrationally excited H₂, (iii) cosmic-ray heating of the gas, and (iv) heating due to turbulent dissipation. 3d-pdr also requires the strength of the incident radiation field, information about any embedded sources, the cosmic ionization rate, and the gas velocity dispersion. See Bisbas et al. (2012) for further technical details. We adopt the umist12 chemical reaction network (McElroy et al. 2013), which uses 215 species and follows approximately 3000 reactions. We use the initial atomic abundances in Table 1 from Sembach et al. (2000). By construction, the gas is initially entirely atomic and neutral.

Table 1.  Atomic Abundances

Species	Abundance Relative to H
H	1.0
He	0.1
C	$1.41\times {10}^{-4}$
N	$7.59\times {10}^{-5}$
O	$3.16\times {10}^{-4}$
S	$1.17\times {10}^{-5}$
Si	$1.51\times {10}^{-5}$
Mg	$1.45\times {10}^{-5}$
Fe	$1.62\times {10}^{-5}$

Note. Atomic abundances adopted from Sembach et al. (2000).

Download table as:  ASCIITypeset image

Each molecular cloud is represented by a one-dimensional slab of constant density. The depth of the cloud is determined by the total molecular gas mass,

$$\begin{eqnarray}&&{R}_{c}={\left(\displaystyle \frac{3{M}_{\mathrm{gas}}}{4\pi n\mu {{\rm{m}}}_{p}}\right)}^{1/3},\end{eqnarray} \tag{ 19 }$$

where m_p is the mass of a proton, n is the gas number density, and μ is the mean molecular weight, taken to be μ = 1.4 since the cloud is assumed to be initially atomic and neutral (see below). The total gas mass is set according to two different gas models as described below.

In these models, there are two FUV components: an external field, ${F}_{\mathrm{ext}}=1$ Draine (Draine 1978), and an internal field, ${F}_{\mathrm{src}}$, from embedded sources as given by the average cluster FUV luminosity from the mock clusters. We scale the latter to the Draine field by renormalizing the units by ${\chi }_{0}=1.7{G}_{0}$, where ${G}_{0}=1.6\times {10}^{-3}$ erg s⁻¹ cm⁻² is the Habing field (Habing 1968). We adopt the fiducial cosmic-ray ionization rate from Bell et al. (2006) of ${\xi }_{0}=1.3\times {10}^{-17}$ s⁻¹ per H₂ molecule. Previous studies of ${{\rm{H}}}_{3}^{+}$ chemistry (Indriolo et al. 2007) and H_nO⁺ chemistry (Indriolo et al. 2015) toward diffuse clouds find larger cosmic-ray ionization rates on the order of 10⁻¹⁶. However, there is a large spread in observed values, and the cosmic-ray ionization rate appears to decrease toward clouds with higher column density (Padovani et al. 2009). We study the implications of cosmic-ray ionization rates higher than the fiducial value in Section 3.

Figure 3 displays a schematic of our cloud model. The field from embedded sources is indicated by blue arrows, and the external field is represented by green arrows. We define A_V, the dust extinction through the cloud, such that the surface has A_V = 0 and the stars are located at high A_V in the cloud center. We place the cluster within an evacuated bubble to approximate the effects of feedback mechanisms. The bubble has a size ${R}_{\mathrm{bubble}}$ given by

$$\begin{eqnarray}&&{R}_{\mathrm{bubble}}=\max (1\,\mathrm{pc},{R}_{s}),\end{eqnarray} \tag{ 20 }$$

where R_s is the Strömgren sphere radius for the given density and ionizing luminosity from the cluster model

$$\begin{eqnarray}&&{R}_{s}={\left(\displaystyle \frac{3{Q}_{0}}{4\pi {\alpha }_{{\rm{B}}}{n}_{e}^{2}}\right)}^{\tfrac{1}{3}},\end{eqnarray} \tag{ 21 }$$

where we approximate ${Q}_{0}=\tfrac{{L}_{\mathrm{Ionizing}}}{18\,\mathrm{eV}}$ following Draine (2011) for first-order computation, ${\alpha }_{{\rm{B}}}$ is the recombination case B coefficient, and we assume ${n}_{e}\approx {n}_{H}$.

Zoom In Zoom Out Reset image size

In addition to the abundances, 3d-pdr computes the line emissivities for C, C⁺, and CO assuming non-local thermodynamic equilibrium (NLTE) and accounting for optical depth. We use these line emissivities to calculate the CO (1–0) integrated line intensity following Röllig et al. (2007):

$$\begin{eqnarray}&&I=\displaystyle \frac{1}{2\pi }{\int }_{0}^{R}j\,{dz}\,\,\,(\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}),\end{eqnarray} \tag{ 22 }$$

with

$$\begin{eqnarray}&&W=\displaystyle \frac{{c}^{3}}{2{k}_{b}{\nu }^{3}}\,I\,\,\,({\rm{K}}\,\mathrm{km}\,{{\rm{s}}}^{-1}),\end{eqnarray} \tag{ 23 }$$

where c is the speed of light, k_b is the Boltzmann constant, and $\nu =115.3\,\mathrm{GHz}$ is the frequency of the CO (1–0) line. We calculate the H₂ column density directly from the 3d-pdr abundances

$$\begin{eqnarray}&&N({{\rm{H}}}_{2})={\int }_{0}^{R}n({{\rm{H}}}_{2})\,{dz}\,\,\,({\mathrm{cm}}^{-2}),\end{eqnarray} \tag{ 24 }$$

where $n({{\rm{H}}}_{2})={n}_{\mathrm{gas}}X({{\rm{H}}}_{2})$ and $X({{\rm{H}}}_{2})$ is the H₂ abundance.

We use two different models for the total gas mass and cloud velocity dispersion. The first, denoted by CM, is a constant-mass model where the total gas mass is ${M}_{\mathrm{gas}}={10}^{4}\,{M}_{\odot }$. This model also assumes a constant velocity dispersion of 1 km s⁻¹, making it slightly subvirial. The second model, denoted by CE, is a constant-efficiency model where the total gas mass depends on the stellar mass: ${M}_{\mathrm{gas}}=\tfrac{{M}_{* }}{{\varepsilon }_{g}}$, where ${\varepsilon }_{g}$ is related to the star formation efficiency,

$$\begin{eqnarray}&&{\varepsilon }_{\mathrm{tot}}=\displaystyle \frac{{M}_{* }}{{M}_{\mathrm{gas}}+{M}_{* }}=\displaystyle \frac{{\varepsilon }_{g}}{{\varepsilon }_{g}+1}.\end{eqnarray} \tag{ 25 }$$

We vary ${\varepsilon }_{g}$ between 0.01 and 0.2, or ${\varepsilon }_{\mathrm{tot}}$ between 0.01 and 0.166. This produces total gas masses from 10³ to 10⁸ M_⊙. We calculate the velocity dispersion for the constant-efficiency models assuming that the clouds are in virial equilibrium, such that

$$\begin{eqnarray}&&{\sigma }_{v}={\left(\displaystyle \frac{4\pi G}{15}\right)}^{1/2}R{\rho }^{1/2},\end{eqnarray} \tag{ 26 }$$

where G is the gravitational constant.

Table 2 summarizes the six models we consider. The fiducial CM model is denoted CM_1000D_1kms_1ξ, and the fiducial CE model is denoted CE_1000D_V_1ξ. We include a model with a lower number density of 500 cm⁻³ (500D), models that vary and fix the velocity dispersion (V and 1 km s⁻¹, respectively), one model with enhanced cosmic-ray ionization rates, and a model without internal sources. This last model allows us to compare the influence of stellar sources relative to the external field.

Table 2.  Chemical Models

Model Name	Constant Mass	Constant Efficiency	Velocity Dispersion	Density (cm−3)	ξ	Internal Sources
CM_1000D_1kms_1ξ	✓		1 km s−1	1000	${\xi }_{0}$	✓
CE_1000D_V_1ξ		✓	Virial	1000	${\xi }_{0}$	✓
CE_500D_V_1ξ		✓	Virial	500	${\xi }_{0}$	✓
CE_1000D_1kms_1ξ		✓	1 km s−1	1000	${\xi }_{0}$	✓
CE_1000D_V_100ξ		✓	Virial	1000	$100{\xi }_{0}$	✓
CE_1000D_V_1ξ_NS		✓	Virial	1000	${\xi }_{0}$

Note. Names and parameters for the different chemical models used. Virial denotes that the velocity is calculated using Equation (26).

Download table as:  ASCIITypeset image

We compare our PLF cluster model to data from three recent surveys of molecular clouds. Dunham et al. (2013) and Kryukova et al. (2012) each survey a number of well-studied clouds located in the Gould Belt. Kryukova et al. (2014) present a survey of the Cygnus X region, which is 1.4 kpc away and is one of the most massive star-forming complexes within 2 kpc of the Sun. Cygnus X contains multiple evolved OB associations with dozens of O stars and hundreds of B stars. It is also the largest cluster in our comparison, with nearly 2000 identified protostars.

The surveys adopt slightly different conventions for identifying protostars. Dunham et al. (2013) define protostars as point sources with at least one detection at λ ≥ 350 μm. They argue that this constraint removes older, non-protostellar sources, while including only sources that are still deeply embedded in dusty envelopes. Kryukova et al. (2012) use color–magnitude diagnostics to identify protostellar sources and do not require a submillimeter detection. Both surveys thus have their own biases: Dunham et al. (2013) likely underestimate the number of dim sources, since protostars embedded in very low-mass cores, which fall below the submillimeter detection limit, are excluded. Kryukova et al. (2012) possibly overestimates the number of low-luminosity sources, by including older, less embedded sources that would have been filtered out by requiring a submillimeter detection. In Chameleon II, however, Kryukova et al. (2012) exclude some of the objects found in Dunham et al. (2013). The net effect is that clusters reported in Kryukova et al. (2012) tend to have larger populations of low-luminosity sources (Dunham et al. 2013). Additional disagreement occurs because the two surveys assume different distances for a few of the shared clouds (i.e., for Perseus the former uses a distance of 230 pc and the latter uses 250 pc). An order-of-magnitude luminosity discrepancy is evident between the two surveys for Chameleon II because the selection criterion in Kryukova et al. (2012) only has one of the three Chameleon II objects in Dunham et al. (2013). Both of the surveys likely suffer from incompleteness at low luminosities to some degree owing to missing sources that are either very low mass (${m}_{* }\lesssim 0.2\,{M}_{\odot }$; Offner & McKee 2011) or very young and embedded ($L\leqslant 0.1\,{L}_{\odot }$; e.g., Maureira et al. 2017)

Given that a significant number of dim protostars could lie below the survey detection limits, we assume that the reported number of sources in both cases is a lower limit that underestimates the true number by up to a factor of 2. This conservative completeness assumption encompasses sources that are either very low mass, e.g., $\lesssim 0.1\,{M}_{\odot }$, or undergoing a period of low accretion. Enoch et al. (2008) cite a 50% completeness limit for the Bolocam 1.1 mm survey, and we use that as an upper limit in the error of the observed protostar number counts to account for incompleteness. Furthermore, while our derived PLF luminosity value is exact, the measured bolometric luminosities have some intrinsic uncertainties that are not reported.

Figure 4 shows the model predictions for the total cluster luminosity across three orders of magnitude in cluster size. We include predictions for the three different accretion models described above. The figure shows the mean total luminosity of the statistically sampled clusters, where dotted lines indicate the 1σ and 2σ deviations from the mean. These boundaries are slightly irregular since they are influenced somewhat by the statistical sampling of the mock clusters. For smaller clusters, a broader PMF creates a correspondingly large spread in the cluster luminosity. The spread decreases for large clusters as the PMF becomes well sampled. For all three models, the total luminosity scales superlinearly with cluster size until ${N}_{* }\approx {10}^{3}$, when it approaches a linear scaling. For the TTC model, the bolometric luminosity is fit by

$$\begin{eqnarray}\mathrm{log}{L}_{\mathrm{Bol}}=\left\{\begin{array}{ll}1.96\cdot \mathrm{log}{N}_{* }+0.18 & \mathrm{log}{N}_{* }\lt 2.78\\ \mathrm{log}{N}_{* }+5.63 & \mathrm{log}{N}_{* }\geqslant 2.78\end{array}\right..\end{eqnarray} \tag{ 27 }$$

Zoom In Zoom Out Reset image size

Inspection of Figure 4 shows that the IS model agrees poorly with the data. It fails to match both the mean and the spread of observed luminosities of the low-mass clusters. However, both the tapered and nontapered TC models are able to reproduce the observed spread quite well. This suggests that poor statistical sampling, together with a significant range of underlying accretion rates, is needed to explain the observational data. All the models appear to significantly overpredict the total luminosity of Cygnus X. However, the brightest sources are saturated in the MIPS 24 μm band, and their luminosities are underestimated in the catalog (R. Gutermuth 2018, private communication).

The TC model does a good job of representing the spread as a function of cluster size, but it overpredicts the luminosities of clusters with sizes ${N}_{* }=10\mbox{--}100$, where observed data points fall outside the 2σ statistical sampling error. The TTC model does exceptionally well in encapsulating the data from all the surveys. The majority of the observed cluster bolometric luminosities are included within the 2σ spread of the model predictions. All models overpredict the luminosities at the smallest cluster sizes. The discrepancy may be caused by several factors, such as completeness limits and differences in the physical parameters we assume, which we discuss in more detail in Appendix A. As a result of this comparison, we adopt TTC as the fiducial accretion model for the analysis in the following sections.

While the total bolometric luminosity is an observable quantity and, thus, useful for evaluating the accuracy of PLF predictions, our PDR calculations require the strength of the FUV radiation field as an input. Since protostellar radiation is heavily reprocessed by the surrounding dusty envelope, it is not possible to directly measure the FUV component of the spectrum. Instead, our PLF models provide an approach to calculate the fraction of short-wavelength radiation. We use the approximation in Equation (17) to calculate the FUV and ionizing luminosity for each protostar in a given cluster and then compute the total by summing over all protostars, i.e., ${L}_{{\rm{\Delta }}E}={\sum }_{i}{L}_{{\rm{\Delta }}E}^{i}$.

Figure 5 shows the PLF model predictions for the total FUV luminosity as a function of cluster size. The TC and TTC models exhibit significant spread in the predicted FUV for modest cluster sizes owing to stochastic sampling of intermediate- and high-mass stars, which contribute most of the FUV radiation. The IS PMF is narrower, which produces slightly better sampling. The spread is magnified in the TC and TTC models, because they assume a broad range of accretion rates. At large cluster masses the luminosity spread diminishes for all three models. All accretion histories show a superlinear trend for small clusters, with the TTC model exhibiting the steepest dependence on cluster size:

$$\begin{eqnarray}\mathrm{log}{L}_{\mathrm{FUV}}=\left\{\begin{array}{ll}3.13\cdot \mathrm{log}{N}_{* }-2.73 & \mathrm{log}{N}_{* }\lt 2.42\\ \mathrm{log}{N}_{* }+4.84 & \mathrm{log}{N}_{* }\geqslant 2.42\end{array}\right..\end{eqnarray} \tag{ 28 }$$

Zoom In Zoom Out Reset image size

Figure 6 shows the total ionizing luminosity, which exhibits a similar trend to the FUV component. Stochastic sampling of the highest-mass protostars (future O and B stars), which are the source of all ionizing radiation, creates larger scatter in the models. Because O stars dominate the budget of ionizing radiation, clusters with ${N}_{* }\lt {10}^{4}$, which do not perfectly sample the high-mass end of the PMF, continue to exhibit a large amount of statistical variation. The steeper slope is due to the strong dependence of accretion rate on stellar mass and the higher peak accretion rates. The TTC model again exhibits the steepest dependence on cluster size:

$$\begin{eqnarray}\mathrm{log}{L}_{\mathrm{ION}}=\left\{\begin{array}{ll}5.4\cdot \mathrm{log}{N}_{* }-8.29 & \mathrm{log}{N}_{* }\lt 2.42\\ \mathrm{log}{N}_{* }+4.78 & \mathrm{log}{N}_{* }\geqslant 2.42\end{array}\right..\end{eqnarray} \tag{ 29 }$$

Zoom In Zoom Out Reset image size

Overall, the models predict that once star formation commences a substantial amount of FUV radiation permeates the natal cloud. For lower-mass clusters the predicted amount of ionization is very small, while a substantial amount of ionizing luminosity is expected in the highest-mass clusters, such as the ONC complex or Cygnus X. In all cases, statistical sampling introduces significant variation, which could drive environmental differences in clouds forming clusters with similar sizes.

In this section, we use model CM_1000D_1kms_1ξ to study the effects of internal embedded sources on the chemical distribution within a cloud. Figure 7 shows the abundance of H₂ and CO as a function of cloud depth and extinction (A_V), where $x/R=0$ is the surface. At low A_V the models for all cluster sizes are similar since the chemistry is dominated by the external radiation field. The H₂ abundances converge to ∼0.5, which indicates that nearly all the H is in H₂.

Zoom In Zoom Out Reset image size

The embedded FUV sources ($x/R=1$) create a shell of H₂, which becomes progressively thinner with increasing cluster size. For ${N}_{* }={10}^{6}$, the H₂ shell is only ∼60% of the total cloud radius. In addition, the amount of CO is reduced even in the region that remains molecular. This is because the column density of material that provides self-shielding is much lower. Consequently, the embedded sources significantly alter the CO abundance profile compared to the typical 1D PDR model. Without embedded sources, the CO abundance asymptotically approaches a value around 10⁻⁴ at high A_V. However, the model predicts that CO is effectively dissociated by ${A}_{V}\geqslant 7$ for all clusters. Increasing the cluster size from N_* ~ 100 to 10⁶ causes 2 orders of magnitude difference in the CO abundance at A_V = 4.

Figure 8 shows the temperature structure of the cloud. 3d-pdr determines the temperature by balancing the heating and cooling as described above. Without embedded sources, the cloud cools to a temperature of 10 K when ${A}_{V}\geqslant 1$. The model results show that the embedded sources have a strong impact, heating the high-A_V gas to hundreds of kelvin. Comparing this temperature structure to the abundance profiles in Figure 7 indicates that there is a large amount of warm CO. These temperatures lead to higher excitation, so more emission preferentially comes from higher rotational levels.

Zoom In Zoom Out Reset image size

The right panel of Figure 8 displays the CO abundance as a function of gas temperature. This phase diagram shows a tight correlation between gas hotter than approximately 50 K and decreasing CO abundance. At cold temperatures, the phase diagram is more complicated owing to the formation and destruction of CO at various points in the cloud.

In this section we investigate how changes in chemistry due to the presence of embedded sources impact the observed CO emission. We control for other factors, including the cluster size, star formation efficiency, turbulent line width, and gas density.

3.3.1. Cluster Size with Fixed Cloud Mass

We first consider the simplest cloud model, model CM_1000D_1kms_1ξ, which holds the cloud mass fixed for all cluster sizes. Figure 9 shows the model predictions for X_CO as a function of the number of stars in the cluster. X_CO approaches 10²⁰ cm⁻² (K km s⁻¹)⁻¹ for small cluster sizes but increases steeply for large clusters. The abundance and temperature profiles in Figures 7 and 8 show the cause of the increase. For a large number of stars, the amount of FUV radiation increases superlinearly, reducing the shell of molecular gas and the column density of H₂. However, due to the high optical depth of the CO (1–0) line, the intensity is dominated by emission near the surface of the cloud. While more CO is dissociated owing to the embedded FUV radiation, the gas also exhibits higher temperatures. These competing factors cancel, producing only a factor of 2 change in X_CO over 4 dex of N_*. This insensitivity to cluster size is encouraging, since it seems to suggest that X_CO is largely invariant. However, our model assumptions break down for large clusters when the stellar mass becomes much greater than the gas mass.

Zoom In Zoom Out Reset image size

3.3.2. Cluster Size with Varying Cloud Mass

Figure 10 shows X_CO as a function of cluster mass, M_*, and the star formation efficiency, ${\varepsilon }_{g}$. Here, M_* is the total protostellar mass (${{\rm{\Sigma }}}_{i}^{{N}_{* }}{m}_{i}$). This figure shows the opposite trend to the constant-mass model shown in Figure 9. For fixed values of the efficiency, X_CO drops by a factor of a few as the cluster mass increases by 4 dex. This mainly occurs as a result of assuming that the molecular cloud is virialized. The corresponding larger line widths increase the integrated CO intensity, causing X_CO to decline. This model is more physically motivated than the simpler constant-mass model; however, it shows that the gas-to-star conversion has a significant impact on the relationship between the column density and CO emission.

Zoom In Zoom Out Reset image size

In Figures 10–14, the solid white contour indicates the average X_CO measured in the MW, X_CO = 2 × 10²⁰ cm⁻² (K km s⁻¹)⁻¹, and the dotted white contours indicate the ±30% error (Bolatto et al. 2013). Our predicted X_CO are consistent with the measured MW values for a large fraction of the parameter space. If we further constrain to look at the region of parameter space encompassing measured star formation efficiencies (see below), the model is consistent for clusters between N_* ≈ 20 and 10⁴. The protostar surveys, mentioned above, span that range of cluster sizes for local star-forming regions where X_CO measurements are best measured.

3.3.3. Star Formation Efficiency

3.3.4. Mean Gas Density

Molecular clouds have a range of mean densities. In this section, we explore the impact of the mean gas density on X_CO. Model CE_500D_V_1ξ is the same model as the fiducial, except with n_H = 500 cm⁻³. Figure 11 shows the same parameter space as the fiducial model shown in Figure 10. The lower density causes X_CO to increase over much of the parameter space. Lowering the density also reduces the column density (and thus the dust extinction), making photochemistry more important. However, changes in the amount of molecular hydrogen and the CO (1–0) emission compete and partially cancel. If there is a reduction in both, X_CO may increase, but not by a large factor. In Figure 11 the overall trend remains the same as the fiducial model but is amplified for moderate and smaller clusters. There is no change for the largest clusters since they have sufficient mass such that changes in the interior abundances occur after the line has become opaque.

Zoom In Zoom Out Reset image size

X_CO for small clusters is greatly amplified owing to the lower dust extinction within the cloud. These clouds have significantly less CO emission compared to their H₂ column density, i.e., they have a larger fraction of "CO-dark" gas (i.e., Wolfire et al. 2010). Typically, CO-dark clouds are assumed to have faint CO emission owing to their low densities, but the gas here is CO deficient owing to dissociation caused by the embedded sources. Consequently, the MW average values occupy only a narrow band across the parameter space. Within the range of typical star formation efficiencies, X_CO is only consistent for clusters with masses between 10³ and 10⁴ M_⊙.

Prior work has also found that X_CO is sensitive to the gas density (Bell et al. 2006; Shetty et al. 2011). Densities higher than 10³ cm⁻³ make dust extinction more efficient, while lower densities enhance the effect of embedded clusters since more of the cloud is influenced by photochemistry. Our models predict that X_CO is most sensitive to density for clouds forming small clusters with high star formation efficiencies. However, this trend may not be evident in observations since diffuse clouds are less likely to form stars with high efficiencies.

3.3.5. Turbulent Velocity Dispersion

The turbulent velocity dispersion is an important factor in the calculated CO emission owing to its influence on the line width and, hence, the optical depth. The line optical depth for a given line of sight is inversely proportional to the velocity dispersion. There is ample evidence that higher-mass clouds have greater velocity dispersions and that many clouds are close to virial equilibrium (Heyer & Dame 2015). Although a constant line width model is unphysical, it is useful to examine the importance of velocity information. In this section, we study the effects of the turbulent velocity dispersion by comparing model CE_1000D_V_1ξ with model CE_1000D_1kms_1ξ.

Figure 12 shows X_CO across the parameter space assuming a constant turbulent line width of 1 km s⁻¹. X_CO exhibits a similar trend to that of the constant-mass model shown in Figure 9. A smaller turbulent velocity dispersion increases the line optical depth, which decreases the overall integrated line flux. The decline in flux, for the same H₂ distribution, increases X_CO. This completely reverses the trend illustrated in Figure 10. Thus, increasing velocity dispersion accounts for much of the decline in X_CO with increasing cloud mass, and the local velocity dispersion is essential to understanding trends in X_CO.

Zoom In Zoom Out Reset image size

3.3.6. Cosmic-ray Ionization Rate

We investigate the effect of the cosmic ionization rate, ξ, which prior work indicates strongly influences X_CO (e.g., Bisbas et al. 2015). X_CO increases with ξ owing to the increased destruction of CO and overall decline in the emission. Higher cosmic-ray fluxes also lead to higher gas temperatures, which in principle could cause X_CO to decline. However, a value of ξ = 100 is not high enough to make cosmic-ray heating the dominant heating mechanism throughout the whole cloud (Bell et al. 2006). Model CE_1000D_V_100ξ adopts a cosmic ionization rate enhanced by a factor of 100 compared to the other models. An increase in the cosmic-ray ionization rate is observed in environments with more star formation, such as those in ULIRGs (Papadopoulos 2010) and toward the central molecular zone of the MW (Le Petit et al. 2016).

Figure 13 shows X_CO for the enhanced cosmic-ray ionization rate. The higher rate increases X_CO by a nearly constant value for all stellar masses. However, the overall trend remains, and the total spread is similar. Since X_CO increases, the fraction of the parameter space consistent with the measured MW values declines.

Zoom In Zoom Out Reset image size

3.3.7. Impact of Internal Sources

To constrain the impact of embedded sources, specifically, on X_CO, model CE_1000D_V_1ξ_NS excludes the star cluster FUV. Figure 14 shows model CE_1000D_V_1ξ_NS with an external field only. Clusters with a mass greater than a few thousand solar masses show almost no difference in X_CO compared to the fiducial model with the inclusion on internal fields. Toward smaller clusters, the model without internal radiation shows an opposite trend. Without the internal FUV radiation, X_CO decreases toward small efficient clusters. Furthermore, X_CO decreases enough that the average MW value is no longer represented in the parameter space. The model values of X_CO within the local star formation efficiency band are only consistent with the lowest measured values.

Zoom In Zoom Out Reset image size

The inclusion of internal FUV radiation increases X_CO for clusters within the sizes indicated in Figure 4 toward MW average values. Large clusters are relatively unaffected because the turbulent line width dominates over chemical effects. Embedded photochemistry only affects the smaller clusters since CO emission is dominated by flux emitted closer to the surface.

Figure 15 shows the linear ratio of X_CO with embedded sources and without them. For most of the parameter space, the embedded sources increase X_CO by 30%–50%. For small efficient clusters, the change is up to a factor of 8, increasing rapidly toward clouds with smaller gas mass. For these clouds, X_CO would likely be time dependent, evolving with the protostellar population.

Zoom In Zoom Out Reset image size

Measuring molecular gas mass in extragalactic sources relies on two dominant methods: dust observations in the infrared and submillimeter and CO emission. In the latter case, the common procedure is to use some approximate conversion factor to calculate the total molecular gas mass within a galaxy. Molecular gas measurements for local galaxies have resolutions of tens to hundreds of parsecs (e.g., Corbelli et al. 2017; Kamenetzky et al. 2017). Furthermore, many of the galaxies targeted are actively star-forming. Since the measured CO integrated flux is an average over the spatially larger star-forming regions, our results suggest that embedded star formation must be taken into account.

Our model results show that for the largest clusters there is little impact from the embedded FUV radiation, and the conversion factor is instead dominated by the turbulent line width. However, for smaller clusters in the range of hundreds to thousands of stars, the embedded radiation has a clear effect. These smaller clusters have X_CO values factors of 3–10 larger than otherwise assumed without the embedded clusters. Furthermore, excluding embedded FUV sources will bias chemical models toward lower densities, higher external radiation, or higher cosmic-ray fluxes. Our X_CO factors presented here are lower limits since our models are one-dimensional constant-density slabs. Real clouds have significant structure and porosity, and the embedded protostars are not tightly grouped into a central cluster, but distributed throughout the cloud. Both of these effects would serve to increase the embedded FUV throughout the cloud, amplifying these trends.

Carbon monoxide has been measured in galaxies out to high redshifts using large single-dish integrated line measurements (e.g., Yun et al. 2015). At these redshifts, the SFR densities are typically much greater than present-day values (Madau & Dickinson 2014). These measurements are dominated by the brightest CO regions, which we predict inherently correspond to lower values of X_CO.

Many molecular gas surveys use dense gas tracers to more directly measure the molecular gas undergoing star formation. Tracers such as HCN and HCO⁺ are the most common alternatives to CO owing to their relatively high abundances (e.g., Gao et al. 2007; Chen et al. 2017). Ammonia (NH₃) is also readily observed in local galaxies (e.g., Lebrón et al. 2011), even though it is associated with dense gas and has a low abundance, because it has a low critical density. Optically thin isotopologues of CO such as ¹³CO and C¹⁸O are often used for line ratio diagnostics (Bolatto et al. 2013). Because they are optically thin, emission from these tracers is sensitive to the conditions of the high-A_V gas, especially molecules such as NH₃ and HCN. Strong FUV radiation from embedded forming star clusters not only dissociates the molecules but heats the gas in the vicinity to hundreds of degrees. Therefore, our work underscores the importance of considering the embedded FUV when modeling optically thin emission from regions expected to have accreting protostars.

In some ways, this is not a novel conclusion. A variety of prior observational work has studied the evolution of gas chemistry near protostars, and observations of high-mass protostars in particular show significant chemical time variation with protostellar evolution (i.e., Hofner et al. 1996; Krumholz et al. 2014). Our results build on the previously acknowledged importance of protostellar feedback to provide a framework for quantifying the impact of feedback on chemistry at cloud scales in addition to the well-studied smaller scales of individual protostars.

A large number of surveys have studied X_CO variation within the MW. For example, Sodroski et al. (1995) and Strong et al. (2004) investigate the radial dependences of X_CO. Both results, although using different methods, conclude that X_CO increases with radius. There are various possible explanations for this trend. Of particular import is the role of the turbulent line width and the mass surface density of clouds. In the center of the galaxy molecular clouds have not only larger masses on average but also larger column densities compared to clouds in the outer regions of the galaxy (Roman-Duval et al. 2010). Furthermore, the overall galactic mass surface density decreases with radius in the MW except for a slight increase around 4 kpc. The SFR surface density also generally decreases except for the same bump at 4 kpc (Kennicutt & Evans 2012). Sodroski et al. (1995) measure X_CO in the center of the MW to be $\mathrm{log}\tfrac{{X}_{\mathrm{CO}}}{{10}^{20}}\approx -0.5$. Our work produces this value for large, inefficient star-forming molecular clouds, especially those subject to a strong cosmic-ray flux. The outer galactic values in Sodroski et al. (1995) are in the range of $\tfrac{{X}_{\mathrm{CO}}}{{10}^{20}}\approx 0.6\mbox{--}1.0$, which is represented in our parameter space by small to intermediate molecular clouds for a mean density of 10³ cm⁻³ or smaller star-forming clouds with a mean density of 500 cm⁻³.

Similar trends are observed in other nearby galaxies. Sandstrom et al. (2013) measured X_CO using high-resolution Herschel maps of 26 nearby disk galaxies. The survey showed that nearly all galaxies exhibit a decrease in X_CO in their inner regions. For example, NGC 6946 is a nearby disk galaxy around 7 Mpc away and one of the galaxies included in Sandstrom et al. (2013). NGC 6946 has also been shown to have a radially decreasing molecular gas and SFR surface density (Kennicutt & Evans 2012). The model parameter space used in this work has a mass surface density that increases from the upper left corner to the lower right corner. Furthermore, Sandstrom et al. (2013) find the average central $\mathrm{log}\tfrac{{X}_{\mathrm{CO}}}{{10}^{20}}\,\approx -0.5-0.2$ consistent with large star-forming molecular clouds in this work. We find that when embedded sources are included we replicate the trends between molecular gas surface density and X_CO in NGC 6947. We stress, however, that this trend does not appear without the FUV radiation from embedded star formation.

Bell et al. (2006) used the ucl-pdr code⁴ (Bell et al. 2005) to perform a parameter study of X_CO as a function of A_V. This cannot be directly related to a cloud-integrated X_CO, but rather to the cloud average A_V, but the trends are similar. Bell et al. (2006) show that in high-density environments X_CO is only weakly affected by the impinging UV field, with X_CO decreasing slightly over 5 orders of magnitude of increasing field. The trend reverses at low A_V, where slight increases in the FUV field significantly increase X_CO. In our study, the external field is fixed while the internal field is increased, and we find that low-A_V gas is relatively unaffected by the embedded sources.

Narayanan et al. (2012) studied the effect of galaxy mergers on X_CO. They compared simulations of quiescent star-forming disks with merging starburst systems and found that the local variation within the galaxy is a smooth function of metallicity. However, starburst systems, which have much higher SFRs, exhibit a lower X_CO factor. They attributed the lower value to an increase in temperature caused by heating from young high-mass stars and larger gas velocity dispersions. This is in good agreement with our model, where the inclusion of FUV radiation from embedded star formation systematically increases the temperature locally, while more massive, turbulent clouds have lower X_CO.

Recently, Clark & Glover (2015) performed simulations to study how X_CO varied with SFR. They fixed bulk properties such as the total mass and initial turbulent field and varied environmental factors that are thought to correlate with SFR. They linearly increased the external FUV field and cosmic-ray ionization rate with the assumed star formation (Papadopoulos 2010). They found that X_CO increases with SFR, contrary to other studies. In fact, their models are similar to our constant-mass model, CM_1000D_1kms_1ξ, (Figure 9), and constant-velocity model, CE_1000D_1ks_1ξ (Figure 12), in which X_CO also increases with cluster mass. Here, this is due to the rapid photodissociation of CO, such that the clouds become CO deficient or rather "CO faint." The constant-mass model is also represented in our constant-efficiency models by using a fixed cluster mass and increasing the efficiency. X_CO increases as a function of star formation efficiency in all of our CE models. Our results show the same qualitative trends as Clark & Glover (2015) when considering models that keep bulk hydrodynamic properties fixed. Keeping the velocity dispersion constant as the star formation activity increases leads to an increasing X_CO as a function of star formation activity.

Previous theoretical studies probed the SFR by changing the external environment. Higher SFR clouds are bathed in stronger FUV fields and in some cases (Lagos et al. 2012; Clark & Glover 2015) experience higher cosmic-ray ionization. The cosmic-ray ionization rate also correlates with the supernova rate and thus the SFR. The result of the scaling between the supernova rate and the SFR creates the linear scalings, $\chi \sim {\chi }_{0}\times \mathrm{SFR}$ and $\xi \sim {\xi }_{0}\times \mathrm{SFR}$, where ${\chi }_{0}$ and ${\xi }_{0}$ are the MW ISRF and cosmic-ray ionization rate, absolutely. The scalings of the impinging FUV radiation and cosmic-ray flux with the SFR apply for galactic-wide studies where a galaxy-averaged SFR is used. On smaller scales these correlations do not hold owing to the star formation activity becoming more stochastic.