Self-consistent population spectral synthesis with FADO I. The importance of nebular emission in modelling star-forming galaxies

Context. Spectral population synthesis (PS) is a fundamental tool in extragalactic research that aims to decipher the assembly history of galaxies from their spectral energy distribution (SED). Whereas this technique has led to key insights into galaxy evolution in recent decades, star formation histories (SFHs) inferred therefrom have been plagued by considerable uncertainties stemming from inherent degeneracies and the fact that until recently all PS codes were restricted to purely stellar fits, neglecting the essential contribution of nebular emission (ne). With the advent of FADO (Fitting Analysis using Differential evolution Optimisation), the now possible self-consistent modelling of stellar and ne opens new routes to the exploration of galaxy SFHs. Aims. The main goal of this study is to quantitatively explore the accuracy to which FADO can recover physical and evolutionary properties of galaxies and compare its output with that from purely stellar PS codes. Methods. FADO and Starlight were applied to synthetic SEDs that track the spectral evolution of stars and gas in extinction-free mock galaxies of solar metallicity that form their stellar mass (M-star) according to different parametric SFHs. Spectral fits were computed for two different set-ups that approximate the spectral range of SDSS and CALIFA (V500) data, using up to seven libraries of simple stellar population spectra in the 0.005-2.5 Z(circle dot) metallicity range. Results. Our analysis indicates that FADO can recover the key physical and evolutionary properties of galaxies, such as M-star and mass- and light-weighted mean age and metallicity, with an accuracy better than 0.2 dex. This is the case even in phases of strongly elevated specific star formation rate (sSFR) and thus with considerable ne contamination (EW(H alpha) > 10(3) angstrom). Likewise, population vectors from FADO adequately recover the mass fraction of stars younger than 10 Myr and older than 1 Gyr (M-star(<10 Myr) / M-star(total) and M-star (> 1 Gyr) / M-star(total), respectively) and reproduce with a high fidelity the observed H alpha luminosity. As for STARLIGHT, our analysis documents a moderately good agreement with theoretical values only for evolutionary phases for which ne drops to low levels (EW(H alpha) <= 60 angstrom) which, depending on the assumed SFH, correspond to an age between similar to 0.1 Gyr and 2-4 Gyr. However, fits with STARLIGHT during phases of high sSFR severely overestimate both M-star and the mass-weighted stellar age, whereas strongly underestimate the light-weighted age and metallicity. Furthermore, our analysis suggests a subtle tendency of STARLIGHT to favour a bi-modal SFH, as well a slightly overestimated M-star(< 10 Myr) / M-star(total), regardless of galaxy age. Whereas the amplitude of these biases can be reduced, depending on the specifics of the fitting procedure (e.g. accuracy and completeness of flagging emission lines, omission of the Balmer and Paschen jump from the fit), they persist even in the idealised case of a line-free SED comprising only stellar and nebular continuum emission. Conclusions. The insights from this study suggest that the neglect of nebular continuum emission in STARLIGHT and similar purely stellar PS codes could systematically impact M-star and SFH estimates for star-forming galaxies. We argue that these biases can be relevant in the study of a range of topics in extragalactic research, including the redshift-dependent slope of the star formation (SF) main sequence, the SF frosting hypothesis, and the regulatory role of supermassive black holes on the global SFH of galaxies.