The Duration of Star Formation in Galactic Giant Molecular Clouds. I. The Great Nebula in Carina

This pilot study uses previously published CCCP data sets, specifically a large subset of the 4664 Chandra X-ray point sources that were positionally matched to MIR sources from the Spitzer Space Telescope Vela-Carina Survey Point Source Archive (Broos et al. 2011b; Povich et al. 2011a).⁶ The Archive provides broadband photometry data at 3.6, 4.5, 5.8, and 8.0 μm, plus JHK_S near-IR (NIR) photometry from the Two Micron All-Sky Survey (2MASS) Point Source Catalog (Skrutskie et al. 2006), which is well matched to both the 2'' resolution of the Spitzer/IRAC detector (Fazio et al. 2004) and the sensitivity limits of the shallow, Spitzer/GLIMPSE-style Legacy surveys ([3.6] ≲ 15.5 mag; Churchwell et al. 2009).

In the process of constructing the Pan-Carina YSO Catalog of Spitzer/IRAC MIR-excess sources, Povich et al. (2011a) identified 3444 IR counterparts to CCCP X-ray sources that were consistent with normally reddened stellar photospheres, plus 213 sources with "marginal" excess emission in the IRAC [5.8] or [8.0] bands. These 3657 sources form the basis of our source sample, as they were possible X-ray-detected members of the CNC but did not exhibit significant 4.5 μm excess emission above a stellar photosphere.

Broos et al. (2011a) classified 75% of X-ray point sources in CCCP as probable members of the Carina Nebula young stellar population based on proximity to observed spatial overdensities, X-ray brightness and median energy, and visual/infrared magnitudes. Kuhn et al. (2019) analyzed proper motion and parallax information from Gaia DR2 (Gaia Collaboration 2018) for 28 young Galactic clusters and associations with available lists of X-ray-selected members, including CCCP, and found that contamination from unrelated field stars was generally <15%. We perform our own cleaning of the CCCP point-source catalog to remove remaining contaminants using an analysis of the Gaia DR2 parallax distribution, similar to that of Kuhn et al. (2019). We do not consider proper motion information because we do not wish to impose a kinematic membership criterion that might exclude high-velocity stars, for example those ejected from one of the massive clusters but still found within the wide CCCP field of view.

We first identified Gaia DR2 matches to CCCP X-ray sources with 2MASS counterparts (Broos et al. 2011b). The nearest Gaia DR2 source falling within 12 of the 2MASS source was declared matched. We then analyzed the parallaxes of 4053 Gaia DR2 matches with G > 8 mag and astrometric_excess_noise < 1.0 mas. The uncertainty on the parallax for each source was adjusted from the DR2 catalog value using the "tentative external calibration" of Lindegren et al. (2018). We then computed the median of the parallax distribution, omitting 417 sources that were >1σ outliers. Our median parallax ϖ₀ = 0.40 ± 0.04 mas gives a distance of ${2.50}_{-0.23}^{+0.28}$ kpc, in agreement with Kuhn et al. (2019) within the uncertainties, which are dominated by our shared assumption of a 0.04 mas systematic uncertainty in the parallax zero point. This moves the Carina nebula to a larger distance than the 2.3 kpc assumed in the original CCCP studies, but this distance value is still consistent with the lower end of the Gaia DR2 distance estimates. Finally, we flagged 294 CCCP sources with individual parallax measurements >3σ above the median (based on their individual adjusted parallax uncertainties) as foreground stars. Of these, 119 were previously classified as probable members of the Carina complex by Broos et al. (2011a). None of the CCCP sources could be analogously flagged as a background star.

Our parallax-based cleaning of probable foreground stars thus removed 7% of all IR-bright CCCP sources. Field stars that were undetected by Gaia or detected but with sufficiently high DR2 parallax uncertainties that they could not be confidently removed may persist as residual contaminating sources.

For this and future papers in this series, we employ a set of 10⁵ "naked" pre-MS spectral energy distribution (SED) models that was first described by Povich et al. (2016, hereafter P16). These models are publicly available.⁷ The naked pre-MS model set is similar to the s-s-i model set from Robitaille (2017, hereafter R17). Both of these model sets consist of spherical stellar photospheres only, with no circumstellar dust in a disk, envelope, or ambient medium. The naked pre-MS models use Castelli & Kurucz (2004) plane-parallel LTE stellar photospheres, the same as the R17 models for T_eff ≥ 4000 K. All R17 model sets describe the central sources solely in terms of stellar radius R_eff and temperature T_eff, with no reference to a particular set of evolutionary tracks and hence no assignment of stellar mass or age. By contrast, for the naked pre-MS models, we sampled the stellar mass M_⋆, and age t_⋆ values were sampled uniformly in logarithmic space (Robitaille et al. 2006, P16) and then converted to (T_eff, R_eff) space using pre-MS evolutionary tracks (Bernasconi & Maeder 1996; Siess et al. 2000). This adds a third derived parameter variable, surface gravity, to the naked pre-MS models that is absent in the s-s-i models, but in practice broadband photometry only constrains the independent variables R_eff and T_eff.

The regions of the traditional Hertzsprung–Russell diagram (HRD; L_bol–T_eff space) and of M_⋆–t_⋆ space sampled by the naked pre-MS models are illustrated in Figure 1. For our goal of constraining age distributions of large stellar populations, both this parameter sampling and the larger size of the naked pre-MS models offer key advantages over the R17 s-s-i models. The uniform sampling in log t_⋆ maps to the familiar, highly nonuniform distribution on the HRD characterized by the thin, densely populated diagonal line of the ZAMS and the nearly vertical pre-MS overdensity. The density of models is greatly reduced between the ZAMS and the pre-MS, a region that has been called the R-C (radiative–convective) gap (Mayne et al. 2007) because it falls between intermediate-mass main-sequence stars with radiative envelopes and their pre-MS analogs that are still fully or partially convective. The region of the HRD to the lower left of the ZAMS is unphysical and hence unpopulated (Figure 1(a)). These models do not include post-main-sequence evolution, creating an unpopulated region in the upper right of the mass–age parameter space (Figure 1(b)). While much effort continues to be devoted to the further development and refinement of pre-MS evolutionary models (e.g., Somers & Pinsonneault 2015; Choi et al. 2016; Dotter 2016; Haemmerlé et al. 2019), all stellar populations synthesized from pre-MS evolutionary tracks share these same general qualitative features. The naked pre-MS SED models therefore possess a built-in, physical "prior" probability distribution of where observed young stars are most likely to be placed on the HRD. This comes at the cost of assuming a particular set of pre-MS evolutionary tracks that may not be correct. The evolutionary models used in our SED models are two decades old, but we adopt them because (1) they have been widely used in the literature, particularly in studies of isochronal ages in the CNC and other Galactic massive star-forming regions to which we will directly compare the results from our new method; and (2) most recent innovations in pre-MS evolutionary models have focused on low-mass stars. We discuss how the choice of different evolutionary models could impact our results in Section 5.2.

Zoom In Zoom Out Reset image size

Figure 2 illustrates the time evolution of the model SEDs for a 1, 2, and 3 M_⊙ star over the first 10 Myr. The 1 M_⊙ model star shows no measurable change in SED shape over this timescale, typical of low-mass pre-MS stars that slowly descend the Hayashi tracks at near-constant T_eff. The intermediate-mass model SEDs, by contrast, exhibit a rapid transition from cool, pre-MS stars to hot, A- or B-type ZAMS stars (between 3 and 5 Myr for the 2 M_⊙ evolutionary track and between 1 and 3 Myr for 3 M_⊙). This change in SED shape is measurable using broadband IR photometry. Specifically, JHK_S photometry can distinguish between the cool pre-MS and hot ZAMS SEDs, constraining T_eff, while Spitzer/IRAC photometry at ≥3.6 μm places strong constraints on the Rayleigh–Jeans tail of the SED and hence R_eff (or L_bol) for a specified T_eff.

Zoom In Zoom Out Reset image size

We fit the SED models to available broadband photometry data using the publicly available Python port (R17)⁸ of the Robitaille et al. (2007) SED fitting tool. All of these CCCP-matched IR sources had been previously fit with different SED models and were found to be consistent with stellar photospheres by Povich et al. (2011a), which imposed the requirement that all sources were detected in at least four of the seven combined 2MASS and IRAC photometry bands, including both the IRAC [3.6] and [4.5] bands. The majority of sources were detected in all five photometry bands from 1 to 4.5 μm. Because the SED declines toward longer IR wavelengths for these diskless sources and the Carina Nebula itself produces a very bright nebular IR background, detection in the IRAC [5.8] and [8.0] bands was less frequent. For the ∼10% of sources with "marginal" excess emission detected at [5.8] or [8.0], these bands were not used for SED fitting.

The extinction curve of Indebetouw et al. (2005) was applied to the SED models as part of the fitting process, which introduced a third independent parameter variable, A_V, to our fitting results. Using the standard J − H versus H − K_S color–color diagram, we determined the maximum extinction to our IR-bright CCCP sources to be A_V = 15 mag and imposed this as a hard upper limit for reddened SED models. We also constrained the scale parameter R_eff/d of the model SEDs by restricting the distance to 2.3 kpc ≤ d ≤ 2.8 kpc.

Following our well-established practice, we use the (unreduced) χ² goodness-of-fit parameter normalized to the number of photometry data points used in the fit (N_data) to identify sources that can be successfully modeled with our chosen model set (e.g., Povich et al. 2011a, 2013; Robitaille 2017). Because the naked SED models are new, we experimented with different values and found well-fit sources to be those for which the single best-fit model satisfied ${\chi }_{0}^{2}$/N_data ≤ 1. We explain our rationale for this choice of cutoff value in Section 3.1.

We define the set i of well-fit SED models to each source using

$$\begin{eqnarray}&&\displaystyle \frac{{\chi }_{i}^{2}-{\chi }_{0}^{2}}{{N}_{\mathrm{data}}}\leqslant 1.\end{eqnarray} \tag{ 1 }$$

This is a more strict cutoff than was used by P16, but still typically generates sets of hundreds or even thousands of well-fit models for a given source.

Interstellar dust makes the photometric colors of reddened, hotter stars nearly indistinguishable from those of unreddened, cooler stars, an unfortunate property of nature that greatly complicates observational stellar astronomy. This manifests itself in our SED modeling as a strong degeneracy between T_eff and A_V. To ameliorate the impact of this degeneracy on our results, we leverage external information about each source in our sample to weight the relative likelihood of each model in the well-fit set.⁹ As described in Appendix B of P16, the relative probability of each individual model i in the set of well-fit models is given by

$$\begin{eqnarray*}&&{P}_{i}({\tau }_{d},{\tau }_{X})={P}_{n}{W}_{i}({\chi }^{2}){W}_{i}({\tau }_{d}){W}_{i}({\tau }_{X})\end{eqnarray*}$$

(their Equation (5)), in which P_n is a normalization constant ensuring that $\sum {P}_{i}=1$. We compute the first two of the following weighting terms using the functions defined by P16. The likelihood of each model based on goodness of fit is W_i(χ²) (Equation (4) of P16). The likelihood that a star of a particular mass and age is diskless is W_i(τ_d), computed for each (M_⋆, t_⋆) parameter combination using Equations (6) and (7) of P16. This weighting function disfavors very young, naked SED models for which we would expect the MIR SED to include emission from a circumstellar dust disk or protostellar envelope. The W_i(τ_d) weighting term is only appropriate for sources for which we can be confident that the stellar photosphere dominates the detected MIR emission at 4.5 μm.

The final weighting term, W_i(τ_X), is the likelihood that the star is X-ray detected based on the model (M_⋆, t_⋆) parameter combination. This term attempts to quantify our knowledge that bright, coronal X-ray emission is a strong indicator of youth. P16 cautioned that their weighting function was strictly appropriate only for low-mass T Tauri stars for which coronal X-ray emission has been well characterized. However, Gregory et al. (2016) reported coronal X-ray emission for pre-MS stars of all masses in a sample that included stars up to 3 M_⊙, the progenitors of main-sequence A-type stars. These authors also demonstrated that the X-ray luminosity decreases more rapidly with age in more massive pre-MS stars, following the time elapsed since the development of a radiative core in partially convective stars. Using the L_X(t) functional relationships for various mass ranges provided by Gregory et al. (2016), we have revised the P16 weighting function based on X-ray detection as follows:

$$\begin{eqnarray}{W}_{i}({\tau }_{X})=\left\{\begin{array}{cc}{\left(\tfrac{{t}_{\star ,i}}{{\tau }_{X}}\right)}^{-{\beta }_{X}}, & {t}_{\star ,i}\gt {\tau }_{X}\\ 1, & {t}_{\star ,i}\leqslant {\tau }_{X}.\end{array}\right.\end{eqnarray} \tag{ 2 }$$

We have introduced two new mass-dependent parameters, τ_X and β_X. The critical age at which a radiative core first develops is

$$\begin{eqnarray}\displaystyle \frac{{\tau }_{X}({M}_{\star })}{\mathrm{yr}}=\left\{\begin{array}{cc}{10}^{6}{\left(\tfrac{1.494}{{M}_{\star }}\right)}^{2.364}, & 1\lt {M}_{\star }/{M}_{\odot }\leqslant 15\\ 2.583\times {10}^{6}, & {M}_{\star }/{M}_{\odot }\leqslant 1.\end{array}\right.\end{eqnarray} \tag{ 3 }$$

Once a pre-MS star evolves beyond the critical age, its X-ray emission decays as a power law governed by

$$\begin{eqnarray}{\beta }_{X}=\left\{\begin{array}{cc}1.19, & 2\lt {M}_{\star }/{M}_{\odot }\leqslant 15\\ 0.86, & 1.5\lt {M}_{\star }/{M}_{\odot }\leqslant 2\\ 0.75, & {M}_{\star }/{M}_{\odot }\leqslant 1.5.\end{array}\right.\end{eqnarray} \tag{ 4 }$$

For low-mass T Tauri stars, β_X is identical to the value reported by Preibisch & Feigelson (2005) and adopted by P16. For the IMPS, however, the much more rapid decay in X-ray emission provides significantly stronger W_i(τ_X) constraints on the SED model fitting results.

After selecting only probable members from Broos et al. (2011a), cleaning out remaining foreground contaminants using Gaia DR2 parallaxes, and imposing our goodness-of-fit criterion, our final sample contains 2269 IR-bright, diskless stellar counterparts to CCCP X-ray point sources. The spatial distribution of these probable intermediate-mass members of the CNC young stellar population is shown in Figure 3. The spatial distribution is complex and hierarchical, exhibiting subclustering on multiple size scales (Kuhn et al. 2014; Buckner et al. 2019). Feigelson et al. (2011) divided the clustered CCCP source population into three principal spatial overdensities (regions A, B, and C, from northwest to southeast) and assigned the remaining sources to a distributed population (region D). To study large-scale variations in the star formation history across the complex, we assign each source in our sample to one of these regions, shown using different-colored symbols enclosed by the lowest density contours in Figure 3. Region A contains the massive Tr 15 and 14 clusters, which are known to have very different ages (e.g., Feinstein et al. 1980; Ascenso et al. 2007), so we bisected this into two subregions, A1 and A2 (separated by the dashed line in Figure 3). This yielded ≳300 sources per subregion (1050 sources in the distributed population), sufficiently large samples for robust statistical analysis of the SED fitting results in each.

Zoom In Zoom Out Reset image size

We utilized this unprecedented, large sample of available spectroscopic classifications for intermediate-mass members of Tr 14 and 16 to benchmark the performance of our SED fitting methodology. For these tests, we included all sources fit with χ₀²/N_data ≤ 4 (expanding our comparison sample to 154 late-type and 19 early-type stars) and performed a second fitting run in which we incorporated V-band photometry from Hur et al. (2012). 

We plotted the T_eff model parameter distributions for each source in the Gaia-ESO comparison sample and compared them to the spectroscopically measured T_eff,S values reported by D17. The underlying probability distributions are typically non-Gaussian. Bimodal parameter distributions are not uncommon, reflecting the degeneracy between a hot photosphere with higher reddening or a cooler star with less reddening. We defined a set of accuracy grades from A (excellent agreement between SED fitting and spectroscopy) to F (complete failure to agree). Definitions for the accuracy grades (ABCDF) are given in Table 1. We consider grades ABC to mean the SED modeling successfully reproduced the true T_eff,S of the star,  with an accuracy of ≤10% for A grades, ≤30% for BC grades, and decreasing precision from A (≤5% relative uncertainty) to C grades (>10% relative uncertainty). DF grades, by contrast, reveal strong or irreconcilable disagreement between the SED models and spectroscopy. D17 could not report spectroscopic T_eff,S for early-type stars, but in all cases our SED models also indicated B-type stars, so we consider them the equivalent of AB grades.

Table 1.  Summary of T_eff Accuracy Grades among the Gaia-ESO Comparison Sample

Grade	NIR(%)	NV+IR(%)	Definition
A	23 (17%)	22 (18%)	Tightly constrained Teff parameter (${\sigma }_{T}/\overline{{T}_{\mathrm{eff}}}\leqslant 0.05$); $\overline{{T}_{\mathrm{eff}}}$ within 10% of Teff,S.
B	28 (21%)	40 (32%)	Well-constrained Teff parameter (${\sigma }_{T}/\overline{{T}_{\mathrm{eff}}}\leqslant 0.1$); $\overline{{T}_{\mathrm{eff}}}$ within 30% of Teff,S.
C	32 (24%)	19 (15%)	Poorly constrained Teff parameter (${\sigma }_{T}/\overline{{T}_{\mathrm{eff}}}\gt 0.1$); $\overline{{T}_{\mathrm{eff}}}$ within 30% of Teff,S.
D	32 (24%)	25 (20%)	Poorly constrained Teff parameter (${\sigma }_{T}/\overline{{T}_{\mathrm{eff}}}\gt 0.1$); $\overline{{T}_{\mathrm{eff}}}$ deviates >30% from Teff,S.
F	10 (7%)	8 (6%)	Well-constrained Teff parameter (${\sigma }_{T}/\overline{{T}_{\mathrm{eff}}}\leqslant 0.1$); $\overline{{T}_{\mathrm{eff}}}$ deviates >30% from Teff,S.
[Early]	14 (10%)	11 (9%)	Spectra and SED modeling agree on B-type star (AB grade equivalent).

Note. $\overline{{T}_{\mathrm{eff}}}$ and σ_T are the P_i(τ_d, τ_X)-weighted mean and standard deviation, respectively, of the T_eff parameter distribution for each source.

Download table as:  ASCIITypeset image

These T_eff accuracy grades guided us to refine our goodness-of-fit criteria (Section 2.4). Among the 154 late-type stars in our expanded spectroscopic comparison sample,  those with best-fit 1 < ${\chi }_{0}^{2}$/N_data ≤ 4 represented only 20% of the sample, indicating a steep decline in fit quality. These relatively poorly fit sources were dominated by DF grades. We therefore adopted ${\chi }_{0}^{2}$ ≤ N_data as the robust cutoff for reliable SED fits. This cutoff produced a final Gaia-ESO comparison sample of 139 stars. The secondary SED fitting run, which included V-band photometry, reduced the number of well-fit sources by 10%.

The distribution of accuracy grades for the IR-only and V+IR SED fitting runs are tabulated in the N_IR and N_V+IR columns, respectively, of Table 1. In the IR-only fitting of the full Gaia-ESO comparison sample, 97 sources (70%) received ABC accuracy grades (this tally includes the 14 early-type stars). Among the well-fit sources in the V+IR fitting run, 92 (74%) received ABC grades. This insignificant gain in the accuracy of the SED fitting in reproducing T_eff was produced by the omission of a small number of DF sources that were not well fit when V-band photometry was included. The inclusion of visual photometry does tend to more tightly constrain the fitting results, as evidenced by the dozen sources that were promoted from C to B grades.

Including V-band photometry makes the SED fitting far more sensitive to the shape of the reddening law than is the case with IR photometry alone. The reddening law is known to vary with location in the CNC, especially in the more obscured regions (e.g., Povich et al. 2011b) that occupy a large area on the sky outside of the lightly reddened window studied by Hur et al. (2012) and D17 (see Figure 3). Because the reddening law itself is not a free parameter in our models, any decrease in SED fit quality at visual wavelengths at higher extinction values is unlikely to be outweighed by gains in precision over fitting the IR SED alone. We therefore do not attempt to incorporate visual photometry (for example, from Gaia DR2) into our subsequent analysis of the full CCCP IR-bright sample.

We chose six example sources to represent each of the five SED accuracy grades plus one early-type star. We plot their SED fits in Figure 5 and present their T_eff distributions in Figure 6. The D17 T_eff,S values are overplotted as dashed green lines for comparison to the P_i(τ_d, τ_X)-weighted model parameters, illustrating the basis for our accuracy grade assignments. The SED fitting results weighted only by goodness of fit (W_i(χ²); red-diamond curves) fail decisively to predict T_eff,S in all of these examples except the early-type star. This is due to the high density of models at cooler temperatures that are inconsistent with the absence of circumstellar disks and envelopes, demonstrating the critical importance of including the W_i(τ_d) parameter from P16.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In nearly all cases where we assigned DF accuracy grades (including the two examples shown in Figure 6), the SED fitting results preferred a too-hot T_eff, equivalent to AB-type stars on the main sequence, while the spectroscopic T_eff,S indicated FGK stars. The simplest explanation for these discrepancies is that we have assumed an incorrect distance in the SED fitting, while the spectroscopic T_eff measurements are distance-independent. Other potential causes for DF grades include erroneous cross matches between our sample and D17 (assumed spectroscopic T_eff is incorrect), unresolved binaries (or confusion between visual and IR sources in crowded areas), or time variability across the several years separating the epochs in which the various visual and IR data sets were obtained.

In Figure 7 we plot Li EW (top panel) and Gaia parallax (bottom panel) versus T_eff,S for the late-type stars in the Gaia-ESO comparison sample, color coded by SED model accuracy grade. We find that 13 of 15 stars with Li EW < 150 mÅ have DF grades (the other two have C grades). Stars with such low Li EW would have not been selected for membership by D17 if they lacked CCCP X-ray counterparts, increasing the chances that they could be foreground contaminants. However, the parallaxes of the large majority of stars with DF grades, including those with Li EW < 150 mÅ, exhibit generally smaller error bars and less scatter about the CCCP median value compared to the distribution of sources with ABC grades, strongly suggesting that they are CNC members, not foreground contaminants.

Zoom In Zoom Out Reset image size

Stars in the cool comparison sample with 5600 K < T_eff,S <8000 K are dominated by DF grades (Figure 7), but their lower Li EW is consistent with their earlier spectral types (A6 through G5; Pecaut & Mamajek 2013), and it is not necessarily indicative of the older ages expected of field contaminants. These stars are also visually brighter than typical of the comparison sample, which may explain their relatively low Gaia parallax uncertainties. They have log L_bol/L_⊙ ≳ 1, placing them in a region of the HR diagram where the density of the naked stellar models is minimal (Figure 1, top panel). The SEDs of stars in this region of parameter space are also well fit by ZAMS models at just slightly higher A_V. Because the model density near the ZAMS is very high, the distribution of well-fit models becomes incorrectly skewed toward higher T_eff values. We suspect this effect is responsible for the majority of the DF grades in the Gaia-ESO comparison sample, and we discuss its impact on our results in the next section.

The aggregate T_eff parameter distribution for the late-type stars in the cleaned spectroscopic comparison sample is provided in Figure 8. This plot shows that our weighted SED fitting parameters (black histogram) successfully reproduce the peak in the spectroscopic temperature distribution near log T_eff/[K] = 3.7. The T_eff distribution weighted only by χ² (red histogram) is unsuccessful, producing a large peak at cooler temperatures that would indicate disk-bearing stars, and a smaller peak at the correct temperature. Irrespective of the weighting function used, SED fitting produces a spurious hump in the distribution with log T_eff/[K] > 4, corresponding to the 30% of sources with DF accuracy grades. 

Zoom In Zoom Out Reset image size

We also construct joint probability distributions of T_eff and L_bol and plot them directly on the traditional HR diagram as two-dimensional histograms. Summing the probability distributions for all of the individual stars in the spectroscopic comparison sample, we obtain the probabilistic HR diagram (pHRD; Figure 9). pHRDs give information similar to a traditional color–magnitude diagram (CMD), and indeed they appear qualitatively similar to the distribution of points on the V₀ versus (V − I)₀ CMD presented by D17 for the entire Gaia-ESO sample (their Figure 13), of which our comparison sample is a subset. The pHRD represents a philosophical shift away from using a single color and magnitude for each individual star to infer the ensemble properties of a larger population. Instead, we synthesize a model for the ensemble population, using a multidimensional set of colors and magnitudes for each constituent star.

Zoom In Zoom Out Reset image size

In Figure 9 we compare pHRDs produced using two different weighting methods for the SED fitting results: (top panel) our standard, age-based likelihood weighting (_Pi(τ_d, τ_X); Section 2.5) and (bottom panel) a Gaussian weighting factor W_i(T_eff,S) using T_eff,S and its uncertainty reported by D17. The weighting function in the latter case is

$$\begin{eqnarray}&&{P}_{i}({T}_{\mathrm{eff},{\rm{S}}})={P}_{n}{W}_{i}({\chi }^{2}){W}_{i}({T}_{\mathrm{eff},{\rm{S}}}),\end{eqnarray} \tag{ 5 }$$

which typically (and unsurprisingly) yields much tighter constraints on all model parameters for individual sources. While both pHRDs place their pre-MS loci between the 1 and 3 Myr isochrones, the shapes of these loci are noticeably different between the two plots. The peak probability in the P_i(τ_d, τ_X)-weighted case forms a narrow ridge with 3.65 <log T_eff/[K] < 3.7, which extends in the log L_bol direction upward from the 2 M_⊙ evolutionary track to nearly the 3 M_⊙ track (red part of the heat map in the left panel, with >8 stars per bin). The pre-MS locus in the _Pi(T_eff,S)-constrained case, by contrast, extends parallel to the isochrones, and its peak is much sharper, located on the 2 M_⊙ track just above the 3 Myr isochrone (right panel, ∼14 stars per bin). Two other features of the pHRD produced using our likelihood weighting (Figure 9, top panel) are missing from the T_eff,S-constrained pHRD (bottom panel): (1) a residual locus of cool (log(T_eff/[K]) < 3.6), low-mass (M_⋆ < 1 M_⊙) pre-MS stars and (2) the intermediate-mass (3–7 M_⊙) ZAMS. Our standard likelihood procedure tends to be relocated with DF grades relocated from the hot end of the pre-MS locus (log T_eff,S/[K] ≈ 3.75) to the ZAMS (see also Figure 7). Because this displacement generally parallels the isochrones, the age determination of the population is unaffected, but the masses of stars with DF grades are overestimated.

Because pre-MS tracks and isochrones are built into the SED models, we can trivially transform the pHRDs (Figure 9) into probability distributions of t_⋆ and M_⋆ (Figures 10 and 11). Compared to the pHRD (or the traditional HRD), these transformations offer more straightforward visualizations of the fundamental stellar parameter distributions. To construct these distributions, we bin the well-fit models uniformly in logarithmic intervals, sampling across the full range of each parameter dimension (Figure 1) and applying the P_i(τ_d, τ_X) weighting functions described in Section 2.5 to each source. We experimented with different numbers of bins for each parameter and chose an optimal binning for the composite histograms that provided the smallest bin size for which fluctuations between adjacent bins did not appear dominated by noise. The final binnings chosen for the 2D pHRDs and associated 1D distributions of M_⋆ and t_⋆ were based on the total number of sources included in each distribution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We create and analyze age and mass distributions (except for the P_i(T_eff,S)-constrained models) using a two-step iterative process. In the first iteration, we identify the modal bin where M_⋆ = M_C1 in the mass distribution produced with our standard P_i(τ_d, τ_X) likelihood weighting. We include only those models with M_⋆ ≥ 1.8 M_⊙ in the t_⋆ distributions (Figure 10). We adopted this cutoff mass to avoid biasing the age distributions toward the youngest pre-MS stars. It also removes the spurious locus of too-cool models with log T_eff/[K] < 3.6 (Figures 7 and 9). The relatively shallow 2MASS and Vela-Carina surveys are insensitive to ZAMS stars fainter than early F dwarfs at the distance of the CNC. Younger, low-mass pre-MS stars are larger and cooler, so they may be detected and included in our sample, while older stars of the same mass fall below our IR photometric sensitivity limits. The existence of fainter, low-mass stars can be inferred from the age distributions of intermediate-mass stars.

In the case of a star-forming event that commenced at time τ_SF (Myr) in the past and proceeded with approximately constant SFR to the present day, logarithmically binned age distributions should exhibit a linearly increasing number of stars per bin as t_⋆ increases from 0 to τ_SF, with a sharp break in the distribution for t_⋆ > τ_SF. To locate this break point in our age distributions, we identify the modal bin in the t_⋆ histogram (Figure 10). The duration of star formation τ_SF is the geometric mean of the central values for this modal bin and the adjacent bin in the direction of increasing t_⋆. The uncertainty on τ_SF is the geometric mean of the widths of these two bins.

Both t_⋆ distributions for the comparison sample give the same value for t_SF = 2.7 ± 0.8 Myr (Figure 10). While spectroscopy sharpens the age distribution, our likelihood-weighted SED fitting to X-ray-selected, diskless sources accurately measures the duration of star formation from broadband IR photometry alone.

Three model mass distributions for the spectroscopic comparison sample are plotted in Figure 11.  The distribution created with our standard age-based weighting function (top panel) follows a Salpeter power law for M_⋆ ≥ 2.7 M_⊙ = M_C, turning over because of photometric incompleteness at lower masses. However, the P_i(T_eff,S)-constrained mass distribution, which is both more precise and more accurate, appears strikingly different. The peak of the distribution narrows, shifting to a marginally lower M_C = 2.3 M_⊙ and producing a significant deficit of stars with M_⋆ > 3 M_⊙ when compared to the Salpeter slope (the area between the solid line and histogram in the bottom panel). Sample bias introduced by our X-ray selection criteria explains this observed deficit. Intermediate-mass (2–8 M_⊙) stars of types A and late B on the MS have no known mechanism for producing detectable X-ray emission and hence occupy an "X-ray desert" between cooler stars with coronal emission (Preibisch et al. 2005) and hotter OB stars producing X-rays via shocks driven by strong winds (Gagné et al. 2011). But those with X-ray-bright T Tauri companions (e.g., Evans et al. 2011) can be selected for inclusion in our IR-bright, diskless sample.

To better constrain the model mass distribution in the general case where no T_eff,S is available, we iterated the analysis of the SED fitting results, introducing a new Gaussian weighting parameter W_i(τ_SF). The complete weighting function for this second iteration is hence

$$\begin{eqnarray}&&{P}_{i}({\tau }_{d},{\tau }_{X},{\tau }_{\mathrm{SF}})={P}_{n}{W}_{i}({\chi }^{2}){W}_{i}({\tau }_{d}){W}_{i}({\tau }_{X}){W}_{i}({\tau }_{\mathrm{SF}}).\end{eqnarray} \tag{ 6 }$$

To determine the present-day mass function of a stellar population from broadband photometry, it is a common practice to construct a dereddened luminosity function and subsequently employ a mass–luminosity relation appropriate to the age of the population (Offner et al. 2014). Here we have adopted τ_SF from the first iteration as the best isochrone for converting from L_bol to M_⋆ for our ensemble population, producing the mass distribution plotted in the center panel of Figure 11.

The two most important metrics for the modeled M_⋆ distributions are the location (M_C) and normalization of the scaled Salpeter-type initial mass function (IMF). Scaling a standard IMF model (Kroupa & Weidner 2003) to the isochronally weighted distribution predicts the same number of intermediate-mass stars as scaling to the spectroscopically constrained distribution (Figure 11, center versus bottom panels) to within 5% (N_IM ≈ 150 for M_⋆ >1.8 M_⊙; Table 2). The same IMF scaled to the mass distribution produced using P_i(τ_d, τ_X) weights alone predicts ∼30% fewer stars. Poorer constraints on the mass of individual stars broadens the aggregate mass distribution, mimicking a shallower IMF slope and reducing the peak height near the cutoff mass M_C (Figure 11, top). We therefore choose to employ the second iteration to leverage isochronal weighting for M_⋆ distributions for the full CCCP IR-bright, diskless source sample.

We caution that our omission of disk-bearing YSOs and X-ray-quiet stars makes this sample formally incomplete at all masses. These mass distributions should not be used to make measurements that require completeness over some mass interval, for example circumstellar disk fractions or the total size of the CNC stellar population. They are mainly useful for drawing comparisons among the various CCCP subpopulations.

Having validated our methodology on the spectroscopic comparison sample, we applied the same analysis to each of our five CCCP spatial subregions (Figure 3). In Figures 12–16, we present the pHRDs, age, and mass distributions for the sources in each subregion. Each pHRD is annotated with the total number N of stars in the sample, and Siess et al. (2000) isochrones and evolutionary tracks are overlaid (same as in Figure 9). Vertical lines overplotted on the age and mass distributions give the duration of star formation τ_SF and break mass M_C, respectively (as in Figures 10 and 11). These and other important quantities describing each subregion are listed in Table 2.

Previous studies have found, based on multiple lines of evidence, that Tr 15 is significantly older than Tr 14 and 16, the other two massive clusters in the CNC (Tapia et al. 2003; Wang et al. 2011). We confirm a significantly earlier onset of star formation in Tr 15 compared to most of the rest of the CNC, with τ_SF = 6.5 ± 1.3 Myr, in good agreement with the measurements reported by Feinstein et al. (1980) from UBVRI photometry and spectroscopy of its brightest members. Tr 15 does not contain a large population of IMPS, as most stars with M_⋆ > 2 M_⊙ have already reached the ZAMS (Figure 12). Disk-bearing, bright YSOs are absent from Tr 15 (Povich et al. 2011a).

Zoom In Zoom Out Reset image size

We find τ_SF = 2.6 ± 0.5 Myr for subregions A2 and B, containing Tr 14 and Tr16, respectively, each presenting a very luminous pre-MS populated by 2–3 M_⊙ IMPS (Figures 13 and 14). Our results agree broadly with previous studies (e.g., Tapia et al. 2003; Ascenso et al. 2007; Hur et al. 2012). Several studies have reported even younger ages <2 Myr for the dense core of Tr 14 (Ascenso et al. 2007; Getman et al. 2014; D17), and we do not dispute these results. A very young age is consistent with the large concentration of luminous YSOs (Povich et al. 2011a) and X-ray-bright IMPS (E. H. Nuñez et al. 2019, in preparation) in Tr 14. Because we selected diskless stars, our sample was biased (by design) toward more-evolved objects. The core of Tr 14 is confused in the 2MASS and Spitzer/IRAC images, so our present analysis can only probe the outskirts of the cluster. Our M_⋆ distributions for both subregions A2 and B break from the Salpeter slope below a relatively high M_C = 3.0 M_⊙, indicating severe photometric incompleteness caused by the combination of crowding (Tr 14), proximity to the very bright IR source η Car (Tr 16), and bright, diffuse background IR nebulosity (both).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Most of the ongoing star formation traced by Spitzer YSOs in the CNC occurs in the South Pillars region, where feedback from massive stars continues to sculpt the remaining GMCs (Smith et al. 2010; Povich et al. 2011a; Roccatagliata et al. 2013). The X-ray-bright, diskless population, however, reveals that the evident embedded population intermingles with a significantly more evolved population (τ_SF = 4.1 ± 0.8 Myr; Figure 15) that predates the Tr 16 and 14 clusters. Ages varying from 1.1 to 4.3 Myr among the various subclusters in the South Pillars were reported by Getman et al. (2014). Individual subclusters display varying disk fractions as well, indicated qualitatively by the varying counts of YSOs versus diskless sources in each one (red versus goldenrod dots in Figure 3). Bochum 11 has few associated YSOs, Collinder 228 has relatively more, and an unnamed, obscured cluster to the immediate southwest of the Treasure Chest hosts many YSOs. The Treasure Chest itself is perhaps the youngest embedded cluster in the CNC (Smith et al. 2005), but its density and high MIR nebulosity precluded detection of its YSO population using Spitzer.

Zoom In Zoom Out Reset image size

Half of the CCCP stellar members lie outside of the three main spatial overdensities in a distributed population (Region D; Feigelson et al. 2011). Not all of this population is truly distributed; some localized overdensities have been grouped into small subclusters (Kuhn et al. 2014; Buckner et al. 2019). YSOs mixed within the distributed population and its most obvious smaller subclusters (Figure 3) reveal that region D traces a wide range of ages, likely spanning the full star-forming history of the CNC. This age range is reflected in the relatively broad peak of the t_⋆ distribution in Figure 16. Compared to the other subregions, we have decreased the histogram bin size and adjusted the break point, yielding τ_SF = 7.2 ± 2.3 Myr for Region D. We find that the duration of star formation in the CNC distributed population, while dominated by the most evolved stars in our CCCP IR-bright, diskless sample, is comparable to the age of the most evolved massive cluster, Tr 15.

Zoom In Zoom Out Reset image size

Our τ_SF results (Table 2) reveal that star formation began at different times at different locations within the CNC, producing a stellar population with a hierarchical, multiclustered structure. Various locations within the CNC continue to host ongoing star formation, as evidenced by protostellar populations (Povich et al. 2011a; Roccatagliata et al. 2013). Our age distributions (with the possible exception of Region A1) would be populated to the lowest measurable t_⋆ if YSOs with IR excess emission were included.

We can piece together a global star formation history of the CNC (Figure 17). The first stars formed 7–10 Myr ago as a GMC with initial mass likely exceeding 10⁶ M_⊙ (Grabelsky et al. 1988; Preibisch et al. 2012). Extrapolating from the remnant morphology apparent in IR dust maps, the structure of the original GMC was probably dominated by one large filament running from southeast to northwest with one or more intersecting, smaller filaments. The initial collapse phase produced, within a few megayears, Tr 15, multiple smaller clusters (including Bochum 11), and a distributed stellar population formed from lower-density regions that were never gravitationally bound.

Zoom In Zoom Out Reset image size

The most spectacular collapse phase commenced <3 Myr ago, as dense gas channeled along the filaments to their main intersection nodes in the center of the nebula rapidly collapsed to form the massive Tr 16¹¹ and 14 clusters (as well as Collinder 232¹² ) in rapid succession. Feedback from the radiation and stellar winds of the very massive stars formed in this event (including η Car in Tr 16 and the O3.5–4 V((f+)) stars in Tr 14; Walborn et al. 2002; Smith & Brooks 2007), and possibly the onset of supernovae (Townsley et al. 2011a, 2011b; Wang et al. 2011) from the most massive stars formed in the initial collapse phase sculpted the remaining molecular material, producing the global, bipolar superbubble structure of the Carina Nebula. The remaining molecular filaments were eroded, forming complex structures in the South Pillars and several pillars in the northern regions of the nebula (Hartigan et al. 2015). The remnants of the original filament–node structure are apparent in the large, V-shaped dust lane immediately south of Tr 16 and 14, and the Great Pillar farther south that points almost precisely toward the vertex of the V (Smith & Brooks 2008).

The older, distributed population is consistent with dynamical evolution of many smaller subclusters over several-megayear timescales, perhaps accelerated by the evaporation of dense gas by massive star feedback. The current, and likely final, phase of star formation in the CNC is underway within the clumpy molecular remnants, possibly triggered by massive star feedback and revealed by very young clusters of YSOs revealed by the retreat of the remaining pillars (Smith et al. 2010). In the near future, multiple supernova explosions from Tr 16 and 14 may remove the remaining dense gas, destroying the Carina Nebula and halting star formation for good.

Previous studies of the star formation history in the CNC using observations at visual wavelengths (e.g., DeGioia-Eastwood et al. 2001) were necessarily restricted to just the central, lightly obscured region containing Tr 14, 15, and 16. But one recent study leveraged CCCP X-ray and NIR data sets to measure stellar ages across a much wider field, probing the more reddened populations. Introducing a new technique called AgeJX, Getman et al. (2014) used a combination of CCCP X-ray data and deep Very Large Telescope/HAWK-I JHK_S photometry (Preibisch et al. 2011) to derive median isochronal ages for each of 19 small subclusters identified by Kuhn et al. (2014) plus the unclustered population. The Getman et al. (2014) study of the CNC was therefore restricted to the smaller HAWK-I field that covered the central 25% of the 1.4 deg² CCCP survey area. AgeJX is complementary to our methodology in many respects. It uses the same Chandra/ACIS X-ray point-source data, assuming an X-ray spectral shape characteristic of low-mass T Tauri stars to assign a mass to each star using an empirical L_X–M_⋆ relation (Preibisch et al. 2005). An extinction-corrected M_J is then used to place stars on the same Siess et al. (2000) pre-MS isochrones used in our naked SED models. AgeJX was restricted to ages <5 Myr and low-mass (M_⋆ < 1.2 M_⊙) stars, requiring deep J-band counterpart photometry for faint X-ray sources that may additionally lie behind large absorbing columns. The overlap between individual AgeJX sources and our CCCP IR-bright diskless sample is minimal.

The number of stars per subcluster used in the AgeJX analysis of the CNC ranged from 6 to 111, with typical clusters containing 20–40 stars suitable for AgeJX dating. Most of the Kuhn et al. (2014) subclusters are too small to contain a statistically robust sample of intermediate-mass stars, but each belongs to one of the larger-scale spatial overdensities (Feigelson et al. 2011) analyzed here. Specifically, region A1 (Tr 15 and environs) contains subclusters G, F, H, and I (AgeJX 2.8–4.8 Myr); Region A2 (Tr 14 and Collinder 232) contains subclusters A–D (AgeJX 1.5–2.8 Myr); Region B (Tr 16) contains subclusters E, J, K, and L (AgeJX 2.4–3.6 Myr); and Region C (South Pillars) contains subclusters M, O, P, Q, R, S, and T (AgeJX 1.1–4.3 Myr). The upper bounds of these AgeJX ranges agree very well with our τ_SF for the Tr 14 and South Pillars regions (Table 2). The oldest AgeJX measurement in Tr 16 is 3.6 ± 0.7 Myr for subcluster K, which is marginally greater than τ_SF. Subcluster K contains η Car and is therefore absent from our analysis, due to MIR saturation, but we find good agreement between AgeJX and τ_SF for the remaining subclusters in Tr 16. The truncation of AgeJX at 5 Myr will cause this method to underestimate the true ages of Tr 15 and the distributed population. In the latter case, AgeJX returns an age of 4.0 Myr based on 354 stars, compared to τ_SF = 7.2 ± 2.3 Myr from 1050 stars (Table 2). Part of this discrepancy may be due to the restriction of AgeJX to the central regions of the CCCP survey area, where even stars found to be unclustered may be younger on average than those in the farther reaches of the distributed population.

We conclude that our results are broadly consistent with past estimates of stellar ages in the CNC. Nearly all of the discrepancies previously reported in the literature can be attributed to difficulties in defining membership of individual clusters, a task made even more challenging by the ubiquitous presence of the large and systematically older, distributed young stellar population.

A number of systematics, in particular neglect of stellar accretion histories and treating binary systems as single stars, both of which we have done here, are known to contribute to apparent isochronal age spreads of ∼10 Myr on HRDs of very young (<20 Myr) star-forming regions (see Soderblom et al. 2014 for a thorough review). That said, we can be reasonably confident in relative measurements of isochronal ages among different stellar subpopulations, and the observed intrinsic age spread across the global CNC stellar population, with its numerous subclusters, is real. The accuracy of absolute isochronal ages, however, is a matter of ongoing investigation and debate. In many evolutionary models, including those used in this study, a single isochrone fails to connect intermediate-mass stars to subsolar-mass pre-MS stars in the same population, predicting younger ages for lower-mass stars (see, e.g., Hillenbrand & White 2004; David et al. 2019).

Relatively few pre-MS evolutionary models provide detailed calculations for the >2 M_⊙ stars that dominate our CCCP IR-bright sample (most exclude this mass range altogether to focus on low-mass pre-MS evolution). David et al. (2019) compare the results of placing the three most massive stellar components of eclipsing binary systems (1.4, 2.6, and 5.6 M_⊙, respectively) in the Upper Scorpius OB Association on isochrones from the Dartmouth (Dotter et al. 2008), MIST (Choi et al. 2016; Dotter 2016), and PARSEC (v1.0; Bressan et al. 2012) evolutionary models. These authors reported good agreement among all isochrones for the ages of these intermediate-mass stars, and while the older Siess et al. (2000) isochrones were not considered, the shapes and locations of their pre-MS isochrones are very similar to the three newer model sets analyzed by David et al. (2019) across the 2–6 M_⊙ mass range.

The pre-MS isochrones and evolutionary tracks adopted in our naked SED models (Bernasconi & Maeder 1996; Siess et al. 2000), in spite of their advanced age, still give reasonable results for intermediate-mass stars. The most promising new grid of evolutionary models appears to be the recent extension of the Geneva evolutionary tracks to the pre-MS phase by Haemmerlé et al. (2019). These models include a detailed treatment of accretion history and computation of the stellar birth line across the entire range of present-day stellar masses. They predict that stars of <3 M_⊙ are fully convective when accretion ends, which is precisely what we have observed with our large sample of X-ray-bright, diskless 2–3 M_⊙ IMPS in the CNC. The locations and shapes of the evolutionary tracks and isochrones for t_⋆ ≥ 1 Myr as well as the ZAMS arrival times reported by Haemmerlé et al. (2019) for 2–5 M_⊙ stars are very similar to those implemented in our models.

Some recent models for low-mass pre-MS evolution introduce the effects of magnetic fields (magnetic pressure support or large starspot covering fractions) that inhibit the contraction of fully convective stars on Hayashi tracks, increasing the isochronal ages of <1 M_⊙ pre-MS stars by factors of ∼2 (Feiden et al. 2015; Somers & Pinsonneault 2015; Feiden 2016; Jeffries et al. 2017). Intrinsic, magnetocoronal X-ray emission characteristic of the cool, convective IMPS in our sample (E. H. Nuñez et al. 2019, in preparation) raises the question of whether a similar correction for isochronal ages of IMPSs might be warranted. We suspect any such correction would be small, however, because a full factor of 2 increase in our reported ages would place the birth of the massive Tr 14 and 16 stellar clusters at 5–6 Myr ago, which exceeds the lifetimes of their most massive stars. This would, in turn, require the massive stars in Tr 16 and Tr 14 to form ≳2 Myr later than the intermediate-mass stars in the same (sub)clusters.

Uncertainties in pre-MS isochronal ages provide the most important systematic uncertainty in our results, but based on our review of the current literature, this systematic appears to be considerably smaller than the factor of ∼2 differences between the magnetic and nonmagnetic low-mass isochrones.

Nuñez et al. (2019, in preparation) present a spectral fitting analysis of the several hundred brightest CCCP X-ray sources (excluding OB stars) and found that those with IR counterparts classified as IMPS on the pHRD were systematically more luminous in X-rays compared to lower-mass T Tauri stars or intermediate-mass AB stars. IMPS extend the convective T Tauri star relation of L_X/L_bol ∼ 10^−3.6 (Preibisch et al. 2005) to higher X-ray luminosities (L_X ∼ 10³² erg s⁻¹). The strongest stellar coronal X-ray flares may thus be produced by IMPS. However, after IMPS complete their descent of the convective Hayashi tracks, they develop radiative cores that grow rapidly, sending stars horizontally across Henyey tracks toward the ZAMS (producing the above-mentioned R-C gap described by Mayne et al. 2007, which is apparent on all of our pHRDs but most pronounced for Tr 14 and 16; see Figures 13 and 14, respectively). The growth of the radiative core drastically alters the magnetic field structure of the star (Gregory et al. 2012). Unlike the case for low-mass stars, the IMPS dynamo-driven X-ray emission decays rapidly after the R-C transition (Gregory et al. 2016) and should disappear completely prior to arrival on the fully radiative ZAMS.

Because a greater fraction of intermediate-mass stars reach the ZAMS and become X-ray dark as a population ages, we expect that the magnitude of this deficit, (the "aridity" of the X-ray desert) will increase as a function of τ_SF.¹³ To explore this effect, for the mass distribution of each region or sample in Table 2, we define a parameter f_XIM, the fraction of X-ray-detected intermediate-mass stars with M_⋆ ≥ 1.8 M_⊙. We estimate the expected number N_IM of intermediate-mass stars in the parent population by integrating the Salpeter IMF (scaled to the value of the M_⋆ distribution at the cutoff mass M_C). After correcting for the varying mass incompleteness in each M_⋆ distribution in the range 1.8 M_⊙ ≤ M_⋆ < M_C using the same scaled Salpeter IMF, we integrate over this corrected distribution to find f_XIM.

As expected, we observe a trend of decreasing f_XIM with increasing τ_SF (Figure 18; error bars on f_XIM are estimated assuming Poisson counting statistics for each bin in the M_⋆ distribution). It would be tempting to use this qualitative trend to infer the real binary fraction among intermediate-mass stars, for example, but systematic effects preclude us from drawing such firm, quantitative conclusions. One complication is that N_IM provides only a lower limit on the true stellar population, because many intermediate-mass stars in the CNC still have disks and were hence not counted among our sample. Furthermore, we have no information about the frequency at which real binary companions are detected as X-ray point sources, the sensitivity of which varies across the large CCCP mosaic (Broos et al. 2011b). But the systematic most relevant to the current analysis involves the shape of the M_⋆ distributions themselves, given the challenges of precisely constraining stellar masses from SED modeling. For the Gaia-ESO spectroscopic comparison sample, we find a deeper deficit in intermediate-mass stars with the more stringent P(T_eff,S) constraints on the SED models than with our standard P(τ_d, τ_x) weighting (f_XIM = 0.78 versus 0.94; see Figure 11 and Table 2). This is purely an artifact of limited precision in measuring stellar masses via broadband photometry alone. A tentative correction for this effect, simply f_XIM − 0.16, is also illustrated in Figure 18. We caution that this correction is not strictly appropriate, given the additional selection criteria imposed by D17 and by us on the Gaia-ESO comparison sample that were not imposed on the full CCCP IR-bright sample, but it does extend the uncertainty on f_XIM to encompass a more plausible range of values.

Zoom In Zoom Out Reset image size