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

TL;DR

This paper introduces a new infrared SED modeling method for young stellar populations in star-forming regions with high extinction, applied to the Carina Nebula, revealing varied star formation durations and recent activity peaks.

Contribution

The paper develops a likelihood-based SED modeling approach for heavily obscured young stars, validated on the Carina Nebula, providing insights into star formation history and public data tools.

Findings

Star formation in Carina began ~10 Myr ago and continues today.

Star formation rate peaked less than 3 Myr ago with cluster formation.

Significant variation in star formation duration among different regions.

Abstract

We present a novel infrared spectral energy distribution (SED) modeling methodology that uses likelihood-based weighting of the model fitting results to construct probabilistic H-R diagrams (pHRD) for X-ray identified, intermediate-mass (2-8 $M_{\odot}$), pre-main sequence young stellar populations. This methodology is designed specifically for application to young stellar populations suffering strong, differential extinction ($\Delta A_V > 10$ mag), typical of Galactic massive star-forming regions. We pilot this technique in the Carina Nebula Complex (CNC) by modeling the 1-8 $\mu$m SEDs of 2269 likely stellar members that exhibit no excess emission from circumstellar dust disks at 4.5 $\mu$m or shorter wavelengths. A subset of ${\sim}100$ intermediate-mass stars in the lightly-obscured Trumpler 14 and 16 clusters have available spectroscopic $T_{\rm eff}$, measured from the Gaia-ESO survey. We correctly identify the stellar temperature in 70% of cases, and the aggregate pHRD for all sources returns the same peak in the stellar age distribution as obtained using the spectroscopic $T_{\rm eff}$. The SED model parameter distributions of stellar mass and evolutionary age reveal significant variation in the duration of star formation among four large-scale stellar overdensities within the CNC and a large distributed stellar population. Star formation began ${\sim}10$ Myr ago and continues to the present day, with the star formation rate peaking ${<}3$ Myr ago when the massive Trumpler 14 and 16 clusters formed. We make public the set of 100,000 SED models generated from standard pre-main sequence evolutionary tracks and our custom software package for generating pHRDs and mass-age distributions from the SED fitting results.

Tables

Table 1. Table 1: Summary of T eff subscript 𝑇 eff T_{\rm eff} accuracy grades among the Gaia-ESO comparison sample

Grade			Definition
A			Tightly-constrained $T_{\mathrm{eff}}$ parameter
B			Well-constrained $T_{\mathrm{eff}}$ parameter
C			Poorly-constrained $T_{\mathrm{eff}}$ parameter
D			Poorly-constrained $T_{\mathrm{eff}}$ parameter
F
[Early]	14		Spectra and SED modeling agree on B-type star (AB grade equivalent).

Table 2. Table 2: Duration of Star Formation in the Various CCCP Sub-Populations

Region/sample	$N$	$N_{t\star }$	${\tau }_{\mathrm{SF}}$ (Myr)	$M_C$ ($M_{☉}$)	$N_C$	$N_{\mathrm{IM}}$	$f_{\mathrm{XIM}}$
Gaia-ESO ($T_{\mathrm{eff},S}$)	aaThis sample includes 14 early-type stars with no $T_{\mathrm{eff},S}$ reported by D17, which we modeled using only our $W_i\left({\chi }^2\right)$ weighting function.		$2.7\pm 0.8$	2.3
Gaia-ESO			$2.7\pm 0.8$
Tr 15	253	124	$6.5\pm 1.3$	1.8	203	335	0.60
Tr 14/Cr 232	356	194	$2.6\pm 0.5$	3.0	124	353	0.85
Tr 16	290bbThe Tr 16 region sample size and inferred stellar population has been significantly reduced by the photometric zone-of-avoidance around $\eta$ Car.	144	$2.6\pm 0.5$	3.0	85bbThe Tr 16 region sample size and inferred stellar population has been significantly reduced by the photometric zone-of-avoidance around $\eta$ Car.	280bbThe Tr 16 region sample size and inferred stellar population has been significantly reduced by the photometric zone-of-avoidance around $\eta$ Car.	0.80
South Pillars	319	156	$4.1\pm 0.8$	2.0	191	297	0.81
Distributed	1050	550	$7.2\pm 2.3$	1.8	760	1184	0.64
All	2269	1168	$3.6\pm 0.6$	1.8	1651	2118	0.77

Equations

$$\frac{{\chi }_i^2-{\chi }_0^2}{N_{\mathrm{data}}}\le 1.$$

$$P_i({\tau }_d,{\tau }_X)=P_n{W}_i({\chi }^2)W_i({\tau }_d)W_i({\tau }_X)$$

$$W_i({\tau }_X)=\{\begin{array}{ll}{\left(\frac{t_{\star ,i}}{{\tau }_X}\right)}^{-{\beta }_X},& t_{\star ,i}>{\tau }_X\\ 1,& t_{\star ,i}\le {\tau }_X.\end{array}$$

$$\frac{{\tau }_X(M_{\star })}{\mathrm{yr}}=\{\begin{array}{ll}10^6{\left(\frac{1.494}{M_{\star }}\right)}^{2.364},& 1<M_{\star }/M_{\odot }\le 15\\ 2.583\times 10^6,& M_{\star }/M_{\odot }\le 1.\end{array}$$

$${\beta }_X=\{\begin{array}{ll}1.19,& 2<M_{\star }/M_{\odot }\le 15\\ 0.86,& 1.5<M_{\star }/M_{\odot }\le 2\\ 0.75,& M_{\star }/M_{\odot }\le \mathrm{1.5.}\end{array}$$

$$P_i(T_{\mathrm{eff},\mathrm{S}})=P_n{W}_i({\chi }^2)W_i(T_{\mathrm{eff},\mathrm{S}}),$$

$$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}}).$$

Full text

methods: data analysis — stars: formation — stars: pre-main-sequence — X-rays: stars — open clusters and associations, individual (Carina Nebula)

\turnoffeditone\turnoffedittwo

1 Introduction

The timescale over which giant molecular clouds (GMCs) collapse and produce new stars places fundamental constraints on theories of star formation (McKee & Ostriker, 2007; Krumholz et al., 2012; Burkhart, 2018) calibrations of star formation rates (SFRs; Chomiuk & Povich, 2011; Kennicutt & Evans, 2012), and the dynamical evolution of H II  regions (Weaver et al., 1977; Koo & McKee, 1992; Arthur et al., 2011; Zamora-Avilés et al., 2019). Massive GMCs produce the O- and early B-type (OB) stars that dominate feedback on the galactic scale, hence measuring the duration of star formation in these structures calibrates the tachometers for the engines driving galaxy evolution (Vogelsberger et al., 2014; Hopkins et al., 2014, 2018).

Our knowledge of constituent stellar populations in the most massive Galactic GMCs has been limited by the combination of relatively large heliocentric distances (${\gtrsim}2$ kpc), relatively large angular extent (typically tens of arcmin to ${>}1^{\circ}$ in diameter) high differential extinction ($\Delta A_{V}>10$ mag), and overwhelming contamination by unassociated field stars in the Galactic plane (Povich et al., 2011a, 2017). Large archival datasets from the Chandra X-ray Observatory and Spitzer Space Telescope have finally flung open an observational window to identify and resolve the young stellar populations associated with the most massive Galactic GMCs (Benjamin et al., 2003; Townsley et al., 2011a; Feigelson et al., 2013; Povich et al., 2013; Wright et al., 2014; Townsley et al., 2014, 2018).

To constrain the duration of star formation, inter-mediate-mass (2–8 $M_{\sun}$) stars are potentially the most valuable, yet still under-utilized, segment of the stellar mass distribution. OB stars are relatively rare and reach the zero-age main sequence (ZAMS) while still embedded in their dusty, natal cores (Zinnecker & Yorke, 2007). The main-sequence turnoff for OB stars gives the age of individual massive clusters (e.g., Figer et al., 1999; Fall et al., 2005), but cannot probe timescales shorter than the ${\sim}3$–20 Myr MS lifetimes of the most massive stars formed in a given cluster and breaks down if massive and low-mass stars are not coeval (Massey & Hunter, 1998). Isochronal ages for low-mass, pre–main-sequence (pre-MS) stars are frequently employed, but these stars evolve relatively slowly along the fully-convective Hayashi tracks, requiring very precise placement on the Hertzsprung-Russell (HR) diagram to establish precise ages. Distant GMCs produce additional challenges, as low-mass T Tauri stars are faint and difficult to detect, and their photospheres are frequently veiled by circumstellar dust in a disk and/or contaminated by accretion signatures such as strong emission lines (Soderblom et al., 2014). By contrast, intermediate-mass, pre-MS stars (IMPS111This class of stars, the cooler progenitors of Herbig Ae/Be stars, have also been called intermediate-mass T Tauri stars (IMTTS; Calvet et al., 2004) and G-type T Tauri stars (GTTS Herbst & Shevchenko, 1999). We introduce the new acronym IMPS because it is more physically descriptive of cool, partially- or fully-convective stars (Palla & Stahler, 1991). We hope to avoid further overworking the T Tauri moniker with additional modifiers, considering that two members of the eponymous triple system, T Tauri Sa and T Tauri N, each have masses of ${\sim}2$ $M_{\sun}$(Köhler et al., 2016), so they are themselves IMPS.) evolve more rapidly than low-mass T Tauri stars. Some fraction of IMPS shed their circumstellar disks in 1–2 Myr (Hillenbrand et al., 1992; Hernández et al., 2007; Povich et al., 2016) to reveal cool, convective photospheres before transitioning to the blue ZAMS in ${<}10$ Myr (Bernasconi & Maeder, 1996; Siess et al., 2000; Haemmerlé et al., 2019). These stars have active dynamos powering magneto-coronal X-ray emission (Gregory et al., 2012, 2016), hence they are identifiable as Chandra point-sources even when the lack of an inner dust disk means there is no measurable mid-infrared (MIR) excess emission above the stellar photosphere. Although they may still possess debris disks, we will henceforth use the shorthand “diskless” to refer to X-ray selected young stars with MIR colors consistent with bare photospheres (for $\lambda\leq 4.5~{}\mu$m).

The wide-field, homogeneous X-ray and infrared datasets produced for the Chandra Carina Complex Project (CCCP; Townsley et al., 2011a) provide an ideal test case for our new methodology. The Carina Nebula Complex (CNC) contains a very large population of ${\gtrsim}5\times 10^{4}$ stars arranged hierarchically among three major, massive clusters (Trumpler 14, 15, and 16), two dozen smaller clusters (including Collinder 228 and 232, Bochum 10 and 11, and the Treasure Chest), and as many as half of the young stars comprise a distributed population (Feigelson et al., 2011; Kuhn et al., 2014). This spatial complexity reflects a complicated star-forming history over at least 10 Myr (DeGioia-Eastwood et al., 2001). The ages of individual subclusters within the CNC span at least several Myr (Feinstein et al., 1980; Carraro, 2002; Tapia et al., 2003; Smith et al., 2005; Ascenso et al., 2007; Hur et al., 2012; Getman et al., 2014) and numerous regions of currently-active, possibly feedback-driven star formation permeate the remaining molecular cloud material (Smith et al., 2010; Povich et al., 2011b; Roccatagliata et al., 2013). The CNC therefore boasts one of the largest, relatively nearby (${<}3$ kpc from the Sun) populations of IMPS, Herbig Ae/Be stars, and intermediate-mass ZAMS stars formed from a single GMC. Our methodology must succeed in reproducing the full range of ages represented by this population while handling the significant differential extinction (negligible $A_{V}$ in certain places, $A_{V}>15$ mag in others, with a reddening law that is known to vary with location and with increasing extinction; see Povich et al. 2011a and references therein).

2 Data, Models, and Source Sample Construction

2.1 Source Datasets

This pilot study uses previously-published CCCP datasets, specifically a large subset of the 4,664 Chandra X-ray point sources that were positionally matched to mid-infrared (MIR) sources from the Spitzer Space Telescope Vela-Carina Survey Point Source Archive (Broos et al., 2011a; Povich et al., 2011b).222Vela-Carina survey mosaics and point-source lists were produced using the GLIMPSE pipeline; for details of the processing and data products, go to http://irsa.ipac.caltech.edu/data/SPITZER/GLIMPSE/doc/glimpse1_dataprod_v2.0.pdf. The Archive provides broadband photometry data at 3.6, 4.5, 5.8, and 8.0 $\mu$m, plus $JHK_{S}$ near-IR (NIR) photometry from the 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]\lesssim 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. (2011b) identified 3,444 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 3,657 sources form the basis of our source sample, as they were possible X-ray detected members of the Carina Nebula Complex but did not exhibit significant 4.5 $\mu$m excess emission above a stellar photosphere.

2.2 Refinement of CCCP Membership List Using Gaia DR2 Parallaxes

Broos et al. (2011b) classified 75% of X-ray point sources in CCCP as probable members of the Carina Nebula young stellar population . 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., 2011a). The nearest Gaia DR2 source falling within 1.2″ of the 2MASS source was declared matched. We then analyzed the parallaxes of 4053 Gaia DR2 matches with $G>8$ mag and $\texttt{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\sigma$ outliers. Our median parallax $\varpi_{0}=0.40\pm 0.04$ mas gives a distance of $2.50^{+0.28}_{-0.23}$ 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 zeropoint. This moves the Carina nebula to a larger distance than the 2.3 kpc assumed in the original CCCP studies, however 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\sigma$ 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. (2011b). 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.

2.3 The SED Models

For this and future papers in this series we employ a set of $10^{5}$ “naked” pre-MS spectral energy distribution (SED) models that was first described by Povich et al. (2016, hereafter P16). These models are publicly-available.333https://doi.org/10.5281/zenodo.2647586 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_{\rm eff}\geq 4,000$ K. All R17 model sets describe the central sources solely in terms of stellar radius $R_{\rm eff}$ and temperature $T_{\rm eff}$, with no reference to a particular set or 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_{\star}$ and age $t_{\star}$ were sampled uniformly in logarithmic space (Robitaille et al., 2006, P16) and then converted to ($T_{\rm eff}$,$R_{\rm 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_{\rm eff}$ and $T_{\rm eff}$.

The regions of the traditional Hertzsprung-Russell diagram (HRD; $L_{\rm bol}$–$T_{\rm eff}$ space) and of $M_{\star}$–$t_{\star}$ 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_{\star}}$ maps to the familiar, highly non-uniform 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 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 1a). These models do not include post–main-sequence evolution, creating an unpopulated region in the upper-right of the mass–age parameter space (Figure 1b). 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.

Figure 2 illustrates the time-evolution of the model SEDs for a 1, 2, and 3 $M_{\sun}$ star over the first 10 Myr. The 1 $M_{\sun}$ 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_{\rm eff}$. The intermediate-mass model SEDs, by constrast, 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_{\sun}$ evolutionary track and between 1 and 3 Myr for 3 $M_{\sun}$). 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_{\rm eff}$, while Spitzer/IRAC photometry at ${\geq}3.6$ $\mu$m places strong constraints on the Rayleigh-Jeans tail of the SED and hence $R_{\rm eff}$ (or $L_{\rm bol}$) for a specified $T_{\rm eff}$.

2.4 Fitting Models to the Data

We fit the SED models to available broadband photometry data using the publicly-available Python port (R17)444http://doi.org/10.5281/zenodo.235786 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 found consistent with stellar photospheres by Povich et al. (2011b), 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–4.5 $\mu$m. Because the SED declines toward longer IR wavelengths for these diskless sources and the Carina Nebula itself produces very bright nebular IR background, detection in the IRAC [5.8] and [8.0] bands was less frequent. For the ${\sim}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_{\rm eff}/d$ of the model SEDs by restricting the distance to $2.3~{}{\rm kpc}\leq d\leq 2.8$ kpc.

Following our well-established practice, we use the (unreduced) $\chi^{2}$ goodness-of-fit parameter normalized to the number of photometry data points used in the fit ($N_{\rm data}$) to identify sources that can be successfully modeled with our chosen model set (e.g. Povich et al., 2011b, 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^{2}_{0}/N_{\rm data}\leq 1$. We explain our rationale for this choice of cutoff value in Section 3.1 below.

2.5 Likelihood-Based Weighting of SED Model Fits

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

[TABLE]

This is a more strict cutoff than was used by P16, but still typically generates sets of hundreds or even thousands of well-fit model parameters 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_{\rm 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.555The software routines and step-by-step procedures used to carry out our advanced analysis of the SED fitting results from this point onwards are publicly-available as an IDL library at https://doi.org/10.5281/zenodo.3234101. We plan to produce a python port in the future. 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

[TABLE]

(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}(\chi^{2})$ (Equation 4 of P16). The likelihood that a star of a particular mass and age is diskless is $W_{i}(\tau_{d})$, computed for each $(M_{\star},t_{\star})$ 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}(\tau_{d})$ weighting term is only appropriate for sources where we can be confident that the stellar photosphere is dominates the detected MIR emission at 4.5 $\mu$m.

The final weighting term, $W_{i}(\tau_{X})$, is the likelihood that the star is X-ray detected based on the model $(M_{\star},t_{\star})$ 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_{\sun}$, 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:

[TABLE]

We have introduced two new mass-dependent parameters, $\tau_{X}$ and $\beta_{X}$. The critical age at which a radiative core first develops is

[TABLE]

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

[TABLE]

For low-mass T Tauri stars, $\beta_{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}(\tau_{X})$ constraints on the SED model fitting results.

2.6 The CCCP IR-Bright, Diskless Source Sample

After selecting only probable members from Broos et al. (2011b), 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 sub-clustering 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 sub-regions, A1 and A2 (separated by the dashed line in Figure 3). This yielded ${\gtrsim}300$ sources per sub-region (1050 sources in the distributed population), sufficiently large samples for robust statistical analysis of the SED fitting results in each.

3 The Gaia-ESO Spectroscopic Comparison Sample

Damiani et al. (2017, hereafter D17) analyzed a large set of spectra, obtained as part of the Gaia-ESO survey (Gilmore et al., 2012), for 1085 stars in a relatively small, lightly-obscured field containing the massive Trumpler (Tr) 14 and 16 clusters (white box in Figure 3). This stellar sample was selected using the optical photometric catalog of Hur et al. (2012). D17 used a combination of CCCP X-ray detection, radial velocity, and Lithium equivalent width (Li EW) $>150$ mÅ to identify 286 likely Carina members with spectral types A or later. D17 published the spectroscopically-determined and uncertainties for all of these stars. Of these, 125 were also found within our IR-bright, diskless sample of CCCP sources (white diamonds in Figure 4).

D17 also identified 110 early-type members for which no $T_{\rm eff,S}$ measurements were possible. Two of these had spectra consistent with obscured O-type stars (previously reported as candidate X-ray emitting O stars by Povich et al. 2011a) while the rest were likely B-type stars. A large majority of the B-type stars were not detected as CCCP sources, presumably because they are X-ray quiet, and only 14 are found in our sample (orange diamonds in Figure 4).666Because D17 did not tabulate CCCP associations for the small minority of their early-type stars that had them, we instead cross-matched their list to our sources using the 2MASS source IDs.

3.1 A Grading Scheme for SED Fitting Accuracy

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 $\chi_{0}^{2}/N_{\rm data}\leq 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_{\rm eff}$ model parameter distributions for each source in the Gaia-ESO comparison sample and compared them to the spectroscopically measured 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 the true $T_{\rm eff,S}$ of the star, with and decreasing precision from A to C grades . DF grades, by contrast, reveal strong or irreconcilable disagreement between the SED models and spectroscopy. D17 could not report spectroscopic $T_{\rm 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.

These $T_{\rm 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^{2}_{0}/N_{\rm data}\leq 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^{2}_{0}\leq N_{\rm data}$ as the robust cutoff for reliable SED fits.

3.2 Discrepancies Between SED Fitting and Spectroscopy

4 Results: Probabilistic H-R Diagrams

Thus far we have examined only the SED fitting results for individual sources. However, we are primarily interested in the aggregate properties of stellar populations. We obtain these by summing the normalized SED model parameter distributions of all constituent sources within the various subsets of interest.

4.1 Temperature and Age Distributions for the Spectroscopic Comparison Sample

The aggregate $T_{\rm 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_{\rm eff}/{\rm[K]}}=3.7$.

We also construct joint probability distributions of $T_{\rm eff}$ and $L_{\rm bol}$ and plot them directly on the traditional H–R diagram as two-dimensional histograms. Summing the probability distributions for all of the individual stars in the spectroscopic comparison sample, we obtain the probabilistic H–R diagram (pHRD; Figure 9). pHRDs give similar information to a traditional color-magnitude diagram, and indeed they appear qualitatively similar to the distribution of points on the $V_{0}$ versus $(V-I)_{0}$ CMD presented by D17 for the entire Gaia-ESO sample (their Fig. 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.

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 (; Section 2.5) and (bottom panel) a Gaussian weighting factor using and its uncertainty reported by D17. The weighting function in the latter case is

[TABLE]

which typically (and unsurprisingly) yields much tighter constraints on all model parameters for individual sources. Two other features of the pHRD produced using our likelihood weighting (Figure 9, top panel) are missing from the $T_{\rm eff,S}$-constrained pHRD (bottom panel): (1) a residual locus of cool ($\log{(T_{\rm eff}/{\rm[K]})}<3.6$), low-mass ($M_{\star}<1$ $M_{\sun}$) pre-MS stars and (2) the intermediate-mass (3–7 $M_{\sun}$) ZAMS.

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_{\star}$ and $M_{\star}$ (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}(\tau_{d},\tau_{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 binsize 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_{\star}$ and $t_{\star}$ were based on the total number of sources included in each distribution.

We create and analyze age and mass distributions (except for the $P_{i}(T_{\rm eff,S})$-constrained models) using a two-step, iterative process. In the first iteration, we identify the modal bin where $M_{\star}=M_{C1}$ in the mass distribution produced with our standard $P_{i}(\tau_{d},\tau_{X})$ likelihood weighting. We include only those models with in the $t_{\star}$ distributions (Figure 10). We adopted this cutoff mass to avoid biasing the age distributions . Younger, low-mass pre-MS stars are larger and and cooler, hence they may be detected and included in our sample while older stars of the same mass fall below our IR photometric sensitivity limits.

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

Both $t_{\star}$ distributions for the comparison sample give the same value for $t_{\rm SF}=2.7\pm 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.

mass distributions for the spectroscopic comparison sample are plotted in Figure 11.

Sample bias introduced by our X-ray selection criteria explains this observed deficit.

introducing a new Gaussian weighting parameter $W_{i}(\tau_{\rm SF})$. The complete weighting function for this second iteration is hence

[TABLE]

For determining 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 to employ a mass–luminosity relation appropriate to the age of the population (Offner et al., 2014). Here we have adopted $\tau_{\rm SF}$ from the first iteration as the best isochrone for converting from $L_{\rm bol}$ to $M_{\star}$ for our ensemble population.

4.2 Age and Mass Distributions for the CCCP Sub-Regions

Having validated our methodology on the spectroscopic comparison sample, we applied the same analysis to each of our five CCCP spatial sub-regions (Figure 3). In Figures 12–16 we present the pHRDs, age, and mass distributions for the sources in each sub-region. 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 $\tau_{\rm SF}$ and break mass $M_{C}$, respectively (as in Figures 10 and 11). These and other important quantities describing each sub-region 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 $\tau_{\rm SF}=6.5\pm 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 IMPS, as most stars with $M_{\star}>2$ $M_{\sun}$have already reached the ZAMS (Figure 12).

We find $\tau_{\rm SF}=2.6\pm 0.5$ Myr for sub-regions A2 and B, containing Tr 14 and Tr16, respectively, each presenting a very luminous pre-MS populated by 2–3 $M_{\sun}$ 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., 2011b) and X-ray bright IMPS (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_{\star}$ distributions for both sub-regions A2 and B break from the Salpeter slope below a relatively high $M_{C}=3.0$ $M_{\sun}$, indicating severe photometric incompleteness caused by the combination of crowding (Tr 14), proximity to the very bright IR source $\eta$ Car (Tr 16), and bright, diffuse background IR nebulosity (both).

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 giant molecular clouds (Smith et al., 2010; Povich et al., 2011b; Roccatagliata et al., 2013). The X-ray bright, diskless population, however, reveals that the evident embedded population intermingles with a significantly more evolved population ($\tau_{\rm SF}=4.1\pm 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.

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_{\star}$ distribution in Figure 16. Compared to the other sub-regions, we have decreased the histogram binsize and adjusted the break point, yielding $\tau_{\rm SF}=7.2\pm 2.3$ Myr for Region D. We find that the duration of start 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.

5 Discussion

5.1 The Complicated Star Formation History of the Carina Nebula Complex

Our $\tau_{\rm 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, multi-clustered structure. Various locations within the CNC continue to host ongoing star formation, as evidenced by protostellar populations (Povich et al., 2011b; Roccatagliata et al., 2013). Our age distributions (with the possible exception of Region A1) would be populated to the lowest measurable $t_{\star}$ 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^{6}$ $M_{\sun}$  (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 Myr, Tr 15, multiple smaller clusters (including Bochum 11), and a distributed stellar population formed from lower-density regions that were never gravitationally bound.

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 16777Tr 16 is in reality a collection of smaller subclusters (Feigelson et al., 2011; Kuhn et al., 2014; Buckner et al., 2019) that are approximately coeval (Getman et al., 2014). and 14 clusters (as well as Collinder 232888Previous studies at visual wavelengths (Tapia et al., 2003; Hur et al., 2012) questioned the existence of the Collinder 232 cluster, probably because parts of it are heavily obscured by a dust pillar. Collinder 232 presents distinct overdensities in both X-ray sources and intermediate-mass YSOs (see Figure 4).) in rapid succession. Feedback from the radiation and stellar winds of the very massive stars formed in this event (including $\eta$ 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 (Wang et al., 2011; Townsley et al., 2011a, b) 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 further 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 sub-clusters over several Myr 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 datasets 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 VLT/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 deg2 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_{\star}$ 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_{\star}<1.2$ $M_{\sun}$) 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 $\tau_{\rm SF}$ for the Tr 14 and South Pillars regions (Table 2). The oldest AgeJX measurement in Tr 16 is $3.6\pm 0.7$ Myr for subcluster K, which is marginally greater than $\tau_{\rm SF}$. Subcluster K contains $\eta$ Car and is therefore absent from our analysis due to MIR saturation, but we find good agreement between AgeJX and $\tau_{\rm 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 $\tau_{\rm SF}=7.2\pm 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.

5.2 Impact of the Choice of Pre-MS Evolutionary Models

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 ${\sim}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_{\sun}$ 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_{\sun}$, 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 intermedate-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_{\sun}$ mass range.

The pre-MS isochrones and evolutionary 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 birthline across the entire range of present-day stellar masses. They predict that stars of ${<}3$ $M_{\sun}$ 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_{\sun}$ IMPS in the CNC. The locations and shapes of the evolutionary tracks and isochrones for $t_{\star}\geq 1$ Myr as well as the ZAMS arrival times reported by Haemmerlé et al. (2019) for 2–5 $M_{\sun}$ stars are very similar 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_{\sun}$ pre-MS stars by factors of ${\sim}2$ (Somers & Pinsonneault, 2015; Feiden et al., 2015; Feiden, 2016; Jeffries et al., 2017). Intrinsic, magneto-coronal X-ray emission characteristic of the cool, convective IMPS in our sample (Nuñez et al. in preparation) raises the question of whether a similar correction for isochronal ages of IMPS 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 ${\gtrsim}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 factors of ${\sim}2$ differences between the magnetic and non-magnetic low-mass isochrones.

5.3 IMPS and the “X-ray Desert”

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_{\rm bol}\sim 10^{-3.6}$ (Preibisch et al., 2005) to higher X-ray luminosities ($L_{X}\sim 10^{32}$ erg s*-1*). 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 $\tau_{\rm SF}$.999For the youngest stellar populations we expect a mass distribution consistent with Salpeter for X-ray selected, diskless IMPS, as observed by P16 for the IR dark cloud M17 SWex. To explore this effect, for the mass distribution of each region/sample in Table 2 we define a parameter $f_{\rm XIM}$, the fraction of X-ray detected intermediate-mass stars with $M_{\star}\geq 1.8$ $M_{\sun}$. We estimate the expected number $N_{\rm IM}$ of intermediate-mass stars in the parent population by integrating the Salpeter IMF (scaled to the value of the $M_{\star}$ distribution at the cutoff mass $M_{C}$). After correcting for the varying mass-incompleteness in each $M_{\star}$ distribution in the range 1.8 $M_{\sun}$ $\leq M_{\star}<M_{C}$ using the same scaled Salpeter IMF, we integrate over this corrected distribution to find $f_{\rm XIM}$.

As expected, we observe a trend of decreasing $f_{\rm XIM}$ with increasing $\tau_{\rm SF}$ (Figure 18; error bars on $f_{\rm XIM}$ are estimated assuming Poisson counting statistics for each bin in the $M_{\star}$ 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_{\rm 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., 2011a). But the systematic most relevant to the current analysis involves the shape of the $M_{\star}$ 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_{\rm eff,S})$ constraints on the SED models than with our standard $P(\tau_{d},\tau_{x})$ weighting ($f_{\rm XIM}=$ versus ; 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_{\rm XIM}-$, 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_{\rm XIM}$ to encompass a more plausible range of values.

6 Conclusions

We have presented a pilot study in which we have developed a novel methodology for constraining the duration of star formation and applied it to multiple young stellar populations observed in the Carina Nebula Complex (CNC). We employ the results of IR SED fitting to individual X-ray-selected, diskless intermedate-mass stars to synthesize models of the ensemble population and probe distributions of stellar masses and ages. Because IMPS rapidly evolve along radiative tracks to become late B- and A-type stars on the ZAMS (including Herbig Ae/Be stars), they serve as the most sensitive available chronometers for the first ${\sim}10$ Myr in the evolution of massive stellar clusters and associations.

Our SED modeling methodology uses all available photometric data simultaneously, incorporating external age constraints based on the lack of mid-IR excess emission and the presence of X-ray emission that are communicated to the model via weighting functions. We treat each star as a cloud of probability in SED model parameter space, analyzing the $\log{L_{\rm bol}}$–$\log{T_{\rm eff}}$ plane (pHRD) coupled with mass and age distributions. This avoids the pitfalls of interpreting individual photometry data points with gaussian error bars, which can be large and are inherently asymmetric with respect to non-linear isochrones and evolutionary tracks. We validate our technique using a comparison sample of stars with spectroscopic $T_{\rm eff,S}$ measurements or broad constraints from the Gaia-ESO survey (D17), and find that we can successfully identify the correct isochronal age from the weighted SED models fit to broadband IR photometry alone. This validation exercise provides a road map for what we would consider a “gold standard” for observational studies of young stellar populations in massive Galactic star-forming regions suffering from severe differential extinction: visual/NIR spectroscopy to measure $T_{\rm eff,S}$ for individual stellar sources combined with IR SED modeling analysis to tightly constrain $L_{\rm bol}$ and reddening.

Our study encompasses the entire 1.4 deg2 field of the CCCP, including the substantial fraction of distributed and/or significantly reddened stars in the CNC, providing a homogeneous analysis of wide-field X-ray and IR imaging data. We find that star formation commenced throughout the CNC ${\sim}10$ Myr in the past (as proposed by DeGioia-Eastwood et al., 2001), after which the SFR rapidly accelerated to produce the massive Tr 15 cluster ${\sim}6.5$ Myr ago and peaked 2–3 Myr ago with the birth of the very massive Tr 16 and 14 clusters, which contain some of the most massive stars known in the Galaxy. Our results generally agree with isochronal ages reported previously for various constituent CNC (sub)clusters (e.g., Feinstein et al., 1980; Carraro, 2002; Tapia et al., 2003; Ascenso et al., 2007; Hur et al., 2012; Getman et al., 2014).

Our X-ray selection criterion imposes unique selection biases on the resultant stellar samples. In particular, the powerful coronal X-ray emission produced by convective IMPS disappears once these objects reach the ZAMS along radiative tracks (Gregory et al., 2016, Nuñez et al. in preparation). We observe the effects of this selection bias as a deficit of intermediate-mass stars in our modeled mass distributions with respect to a standard Salpeter IMF slope, the magnitude of which increases with age of the stellar population. Late B- and A-type stars (including Herbig Ae/Be stars) near the ZAMS can also be included in our samples if they have lower-mass, convective stellar companions producing detectable X-ray emission.

Subsequent papers in this series will implement our methodology for a uniform, comparative study of the duration of star-formation in more than two dozen other Galactic massive star-forming regions for which comparable X-ray and IR photometric data exist. The results of this extended study will enable improved measurements of star-formation rates, star-forming efficiencies, and the dynamical evolution of giant H II regions.

We thank M. A. Kuhn and L. A. Hillenbrand for numerous discussions and suggestions that substantially improved this work. This research was supported by the NSF through grant CAREER-1454333 (PI M. S. Povich). JTM and EHN acknowledge support from the Cal-Bridge program through NSF awards DUE-1356133 and AST-1559559. The scientific results are based in part on observations made by the Chandra X-ray Observatory and published previously in cited articles. This work is based in part on archival data obtained with the Spitzer Space Telescope, which is operated by the Jet Propulsion Laboratory, California Institute of Technology under a contract with NASA. This publication makes use of data products from the Two Micron All-Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by NASA and the NSF. CXO (ACIS), Spitzer (IRAC), CTIO:2MASS