Cloudy-Maraston: Integrating nebular continuum and line emission with the Maraston stellar population synthesis models

The James Webb Space Telescope has ushered in an era of abundant high-redshift observations of young stellar populations characterized by strong emission lines, motivating us to integrate nebular emission into the new Maraston stellar population model which incorporates the latest Geneva stellar evolutionary tracks for massive stars with rotation. We use the photoionization code Cloudy to obtain the emergent nebular continuum and line emission for a range of modelling parameters, then compare our results to observations on various emission line diagnostic diagrams. We carry out a detailed comparison with several other models in the literature assuming different input physics, including modified prescriptions for stellar evolution and the inclusion of binary stars, and find close agreement in the H$\rm β$, H$\rm α$, [N II]$λ6583$, and [S II]$λ6731$ luminosities between the models. However, we find significant differences in lines with high ionization energies, such as He II$λ$1640 and [O III]$λ5007$, due to large variations in the hard ionizing photon production rates. The models differ by a maximum of $\hat{Q}_{\rm [O III]λ5007} = \rm 6 \times 10^9 \; s^{-1} \, M_{\odot}^{-1}$, where these differences are mostly caused by the assumed stellar rotation and effective temperatures for the Wolf Rayet phase. Interestingly, rotation and uncorrected effective temperatures in our single star population models alone generate [O III] ionizing photon production rates higher than models including binary stars with ages between 1 to 8 Myr. These differences highlight the dependence of derived properties from SED fitting on the assumed model, as well as the sensitivity of predictions from cosmological simulations.