2019ApJ...881...37P

We present a novel infrared spectral energy distribution (SED) modeling methodology that uses likelihood-based weighting of the model fitting results to construct probabilistic Hertzsprung-Russell diagrams (pHRD) for X-ray-identified, intermediate-mass (2-8 M_☉), pre-main-sequence young stellar populations. This methodology is designed specifically for application to young stellar populations suffering strong, differential extinction (Δ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 µm SEDs of 2269 likely stellar members that exhibit no excess emission from circumstellar dust disks at 4.5 µm or shorter wavelengths. A subset of ∼100 intermediate-mass stars in the lightly obscured Trumpler 14 and 16 clusters have available spectroscopic T_eff, measured from the Gaia-ESO survey. We correctly identify the stellar temperature in 85% of cases, and the aggregate pHRD for all sources returns the same peak in the stellar age distribution as obtained using the spectroscopic T_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 ∼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.