© ESO 2021

1. Introduction

The study of galaxy evolution will see tremendous breakthroughs in the next decade as a consequence of the huge effort by the scientific community using multi-object spectroscopic facilities. Surveys carried out with instruments such as 4MOST (de Jong et al. 2019), WEAVE (Dalton 2016), MOONS (Cirasuolo & MOONS Consortium 2016; Cirasuolo et al. 2020; Maiolino et al. 2020), and PFS (Takada et al. 2014) will enormously increase the number of available galaxy spectra, opening new windows into unexplored areas of research.

These spectra will tackle the crucial problem of the assembly of galaxies, shedding light on the physical processes shaping galaxies throughout their evolution. From this perspective, it is of paramount importance to refine and test reliable analysis tools that are able to infer the properties of a galaxy from medium- or high-resolution spectra. To achieve this goal a widely used technique is so-called spectral synthesis, which comes in two flavours: evolutionary and population synthesis methods (see the seminal works of Struck-Marcell & Tinsley 1978 and Faber 1972).

Evolutionary synthesis methods estimate the evolution of the spectral energy distribution (SED) of galaxies according to a specific set of star formation histories (SFH) and chemical evolution models (CEM). The observed spectrum is compared to different sets of models, and the best-fit provides the physical parameters of interest. The main limitation of such methods is in the underlying assumption of a well-defined star formation law, which even with more complex parametrisations (e.g., delayed exponential, double exponential), or SFHs derived from cosmological simulations (e.g., Pacifici et al. 2015), lacks the needed accuracy to reproduce the complexity of star formation processes within a galaxy (see a more complete discussion in Gomes & Papaderos 2017).

Population synthesis methods instead do not assume a fixed SFH and/or CEM, since by definition they attempt to approximate the SFH of a galaxy as a linear superposition of simple stellar populations (SSPs), each one corresponding to the SED of an instantaneously formed stellar population of a given age, metallicity, and initial mass function (IMF; Tinsley 1980; Bruzual 2010). This approach is complementary to evolutionary methods and has a number of advantages. The most important advantage is that this approach is free of prior assumptions about the SFH and CEM. This allows the investigation of galaxies with complex SFHs: sudden bursts due to a merger, active galactic nuclei (AGNs) triggering, or the abrupt decline of star formation due to ram pressure stripping (e.g., Cardoso et al. 2017; Pappalardo et al. 2010).

The main limits of spectral synthesis codes are the following: (a) degeneracies inherent to theoretical stellar physics; (b) different mathematical recipes of the available codes, which have their own intrinsic limitations; (c) the non-perfect coverage of the adopted SSP library, resulting in increasing uncertainties in modelling techniques; and, (d) the neglect of important physical ingredients to the SED (e.g., nebular continuum emission), as pointed out in Gomes & Papaderos (2017).

Different approaches have been proposed to overcome this limitation. An interesting example is to join semi-analytical libraries of stars to reproduce the spectra in missing regions of the parameter space (Coelho et al. 2020). This sparse coverage can introduce large errors; depending on the application, semi-empirical or fully theoretical models can be effective. With the development of more advanced computational techniques, analysis methods have become more sophisticated, and new approaches such as genetic algorithms (Gomes & Papaderos 2017), principal component analysis (PCA, McGurk et al. 2010), mixed PCA-neural networks analysis (Alsing et al. 2020), or Monte Carlo methods (Leja et al. 2017) have been implemented. Despite the significant improvements provided by these approaches, the SFH and CEM remain complicated elements to disentangle, and different methods have been refined to break the well-known problem of age-metallicity-extinction degeneracy.

For example in the STEllar Content and Kinematics via Maximum A Posteriori likelihood tool (STECKMAP, Ocvirk et al. 2006a,b), the matrix solution is regularised through the introduction of a penalty function, while in Iyer et al. (2019) the star formation is modelled by adding a series of Gaussian bursts, reproducing various star formation episodes.

A different approach is taken in STARLIGHT (Cid Fernandes et al. 2005), where an observed spectrum is fit with a linear combination of SSPs of several ages and metallicities, assuming a given extinction law and a Gaussian distribution for stellar velocities. The results are very accurate, recovering reliable stellar ages and metallicities but, as noted in Cid Fernandes et al. (2005), these analyses are valid “as long as one does not attempt a very detailed description of the star formation and chemical histories”. For this reason subsequent versions of the code have tackled different aspects of the issues mentioned in their paper, including photometric data from ultraviolet (UV, see e.g., López Fernández et al. 2016). This approach is not free of caveats, one of which is that aperture effects and differences in the point spread function (PSF) of GALEX and optical spectrographs need to be accounted for.

The problem of determining the SFH and CEM in a galaxy from spectra remains hard to solve, and the main outcome from these previous investigations is that only for high-resolution spectra with sufficient signal-to-noise ratio (S/N ∼ 10) it is possible to achieve reliable results. This is because it is only possible to disentangle two bursts with 0.8 dex difference in age at high S/N, as shown in Ocvirk et al. (2006a). However, this effect depends on the ages considered and could also be a consequence of the regularisation itself and thus specific to STECKMAP. This is because the penalty function tends to avoid solutions with high derivative, excluding short bursts and favouring a more continuous SFH. Beyond the quality of the data, we should also consider the quality of the models used to fit the data, which has a crucial impact on the final outcome of the analysis (Koleva et al. 2008).

In this context, the Fitting Analysis using Differential evolution Optimization tool (FADO, Gomes & Papaderos 2017) represents a conceptually novel approach because it is the first code to reproduce the nebular continuum in star-forming galaxies. This approach is of paramount importance to quantify the role of warm ionised gas in star-forming regions, particularly in starburst systems. Dealing self-consistently with nebular and stellar emission allows FADO to reliably reproduce emission fluxes and equivalent widths (EW) during starburst phases of galaxies and to identify the stellar continuum in UV emitting regions, thereby removing the bias introduced by neglecting the nebular gas (Gomes & Papaderos 2017). A correct modelling of this component will become more relevant in the near future because upcoming surveys will provide us with many more galaxies at higher redshift, where we will analyse spectra that contain or are dominated by nebular emission. This is because the rest-frame blue part of the spectrum of a galaxy will be redshifted into optical spectral window and also because star-forming galaxies become more common as we move to higher redshifts.

In light of these future surveys, it is necessary to understand and quantify the reliability of the results obtained with such methods by determining their ability to recover fundamental evolutionary parameters. Motivated by these problems, and having in mind the plethora of multi-object spectrograph instruments in the deployment phase, in this paper we aim to investigate the limits and reliability of different spectral synthesis methods in the estimation of two specific quantities being representatives of the evolutionary status of a galaxy: the mean stellar age and mean stellar metallicity. These two parameters are fundamental for determining the assembly history of a galaxy, by providing key insights into its SFH and stellar mass growth. The latter will be the subject of a forthcoming article of this series.

The main question we want to address in this work is at which S/N values is it still possible to determine, from a galaxy spectrum, its mean stellar age and metallicity and how this depends on the tool used to model the spectra.

To address this question we built a set of simulated spectra reproducing the evolution of a galaxy considering two limiting cases: a continuous star formation and a single initial burst with an exponentially declining star formation law. We degraded the synthetic spectra built from these two SFHs to different S/N, simulating the different conditions encountered within the different surveys mentioned at the beginning of Sect. 1.

The obtained mock spectra were then analysed with three widely used spectral synthesis codes: FADO, STECKMAP, and STARLIGHT assuming similar initial set-ups and spectral bases. The main goals are as follows: (a) to determine at which S/N the results for the parameters considered are still reliable; (b) to identify which stellar ages the fitting tools are more sensitive to; and (c) to quantify the effect of nebular emission component on the physical quantities inferred for the stellar component in population synthesis codes.

The paper is organised as follows: Sections 2 and 3 describe the construction of synthetic spectra and the spectral fitting tools used for the analysis, respectively. Sections 4 and 5 report and discuss the results, while Sect. 6 presents the main conclusion of this work.

2. Synthetic spectra: Dancing REBETIKO

The synthetic spectra analysed in this paper were built with the evolutionary synthesis code Reckoning galaxy Emission By means of Evolutionary Tasks with Input Key Observables (REBETIKO), applying the following equation:

$$\begin{array}{cc}& F(\lambda ,t)={\displaystyle {\int }_0^t}\mathrm{\Phi }(t-t^{\prime })F_{\lambda }^{\mathrm{SSP}}(t^{\prime },Z)\mathrm{d}{t}^{\prime }\hfill \\ \hfill & +\frac{{\gamma }_{\mathrm{eff}}(T_e)}{{\alpha }_B(T_e)}{\displaystyle {\int }_0^t}\mathrm{\Phi }(t-t^{\prime })q^{\mathrm{SSP}}(t^{\prime },Z)\mathrm{d}{t}^{\prime }.\hfill \end{array}$$(1)

Equation (1) reproduces the evolution of the SED of a galaxy as a function of time. For each t the flux F(λ, t) at a specific wavelength λ depends on two terms: The first term represents the emission due to the stellar component, which is the integral of the different SSPs weighted for the star formation rate (SFR; Φ). The second term represents the flux contribution due to the nebular continuum, which depends on the number of hydrogen ionising photons q^SSP weighted for the star formation (analogously to the stellar component). The only difference here is the addition of an extra term quantifying the efficiency of the process, which is estimated by the ratio of the photons continuous emission γ_eff to the photons recombination rate, α_B. These values depend on the gas temperature and the star formation-weighted nebular emission since both processes drive the number of photons emitted (Schaerer & de Barros 2009). We assume all ionising photons produced by the stars are absorbed by the interstellar medium (ISM), and the nebular continuum is computed assuming case B recombination under typical physical conditions of HII regions, that is an electron density of n_e = 100 cm⁻³ and a temperature of T_e = 10⁴ K (Osterbrock & Ferland 2006; Osterbrock 1989).

In the construction of the models, we did not take into account potential misalignments of stellar and gas kinematics. Decoupling of gas and stars in a galaxy can be indicative of gas accretion (Sarzi et al. 2007; Chen et al. 2019), AGN or starburst-driven feedback (Florian et al. 2020; Cicone et al. 2016), and/or gas loss due to environmental effects (Duckworth et al. 2020). The misalignment affects mainly early-type galaxies and only a small fraction of star-forming discs (Jin et al. 2016). For such cases, the integral field unit analysis of Martinsson et al. 2013 found a misalignment on the order of 2 km s⁻¹. Our mock galaxies are built assuming a closed-box model neglecting external processes, such as accretion, ram pressure stripping, and merging. Finally, emission lines are masked out in fits with STECKMAP or STARLIGHT, thus a kinematical decoupling between gas and stars does not notably influence these fits. Nebular continuum emission (excluding its Balmer and Paschen jump) is featureless, so a kinematical decoupling does not play a role. As for FADO, it measures the Hα and Hβ (within ±400 km s⁻¹ from the stellar velocity) and uses their luminosity to self-consistently infer the SFH, so kinematical decoupling is not an issue.

One of the goals of this paper is to quantify the effect of the inclusion of nebular continuum in the determination of evolutionary parameters inferred by the different codes and describing the evolution of a galaxy. To accomplish such an objective the strategy is to estimate the reliability of the results obtained with different analysis methods as a function of S/N. For these reasons we built our synthetic spectra assuming two SFHs. First, the Cont model simulates a galaxy undergoing a continuous accretion of gas, which is therefore able to sustain star formation at a constant rate. Second, the TAU1 model simulates a galaxy undergoing an initial instantaneous burst, on a timescale of ∼10⁸ Myr, followed by an exponentially declining SFR (Φ = e^−t/τ, with τ = 1 Myr).

From the mathematical point of view extracting at the same time age and metallicity from a spectrum is an ‘ill-posed’ problem because of the well-known age-metallicity degeneracy. The introduction of a complex chemical evolutionary history for our models would complicate the interpretation of the results, and for this reason, we assumed for our models a constant solar metallicity with Z_⊙ = 0.02. We note that Bruzual & Charlot (2003) made use of the STELIB stellar library (Le Borgne et al. 2003), which is biased towards solar metallicity. The spectra of the simulated galaxies cover 3400–8900 Å and have a spectral resolution of 2.3 Å and a Δλ = 1 Å. The wavelengths coverage emulates the average range observed by the Sloan Digital Sky Survey (SDSS, Ahumada et al. 2020; Doi et al. 2010), including the Balmer and Paschen jumps, which are relevant features at young ages (<10⁷ yr) and are impossible to fit with purely stellar templates (Cardoso et al. 2019). The evolution of a galaxy is divided into 716 time-steps spanning 1 <  t <  1500 Myr, assuming as a spectral basis the 221 ages in Bruzual & Charlot (2003), a Chabrier (2003) IMF, and Padova 1994 evolutionary tracks (Alongi et al. 1993; Bressan et al. 1993; Fagotto et al. 1994a,b; Girardi et al. 1996).

The models assume no extinction (V-band extinction of A_V = 0 mag) and no extra line broadening from stellar kinematics to reduce the parameters space studied. The latter is a necessary step in spectral fitting methods, usually estimated convolving the models with a proper kinematical kernel, including physical and instrumental effects. This procedure allows us to improve the constraints on the galaxy kinematics, spot an anomalous PSF, and determine the redshift with higher accuracy when unknown. The REBETIKO spectra for the CONT and TAU1 models are shown in Fig. 1 (adapted from Cardoso et al. 2019).

Fig. 1. Synthetic spectra normalised at λ = 4020 Å for TAU1 (instantaneous burst, top panel) and CONT (continuous, bottom panel) SFH models. Figures adapted from Figs. 1 and. 2 of Cardoso et al. 2019.

The evolution of the synthetic SEDs simulated in REBETIKO is characterised by two parameters: the mean stellar age and the mean stellar metallicity. Each of these parameters is weighted by the light (mass) contribution to the total galaxy light (mass), considering logarithmic time-steps, as in Cid Fernandes et al. 2005:

$$\begin{array}{cc}& {⟨log(t_{\ast })⟩}_L={\displaystyle \sum _{j=1}^{N_{\ast }}}{x}_j·log(t_j),\hfill \\ \hfill & {⟨log(t_{\ast })⟩}_M={\displaystyle \sum _{j=1}^{N_{\ast }}}{\mathrm{\mu }}_j·log(t_j).\hfill \end{array}$$(2)

Analogously for the metallicity,

$$\begin{array}{cc}& {⟨log(Z_{\ast })⟩}_L={\displaystyle \sum _{j=1}^{N_{\ast }}}{x}_j·log(Z_j),\hfill \\ \hfill & {⟨log(Z_{\ast })⟩}_M={\displaystyle \sum _{j=1}^{N_{\ast }}}{\mathrm{\mu }}_j·log(Z_j).\hfill \end{array}$$(3)

In Eqs. (2) and (3), t_j (Z_j) represents the age (metallicity) of the jth SSP element, where x_j and μ_j are their associated light and mass fractions at 4020 Å, respectively. Specifically, μ_j denotes the existent stellar mass M_* corrected for the fraction returned to the ISM during stellar evolution. These parameters characterise the stellar population mixture of a galaxy, tracing the main contributors to the observed light and mass. To analyse more quantitatively the evolution of these parameters we defined

$$\begin{array}{cc}& {\mathrm{\Delta }}_L={⟨log(t_{\ast })⟩}_L-log(t)\hfill \\ \hfill & {\mathrm{\Delta }}_M={⟨log(t_{\ast })⟩}_M-log(t),\hfill \end{array}$$(4)

which is the ratio of the mean light (mass) weighted stellar age to the galaxy age. The parameters in Eq. (4) quantify the uncertainties related to the use of average quantities, as in Eq. (2), to estimate the age of a galaxy. Figure 2 shows the evolution of such parameters for the CONT model.

Fig. 2. Differences between luminosity-weighted (solid black line) and mass-weighted (dashed red line) mean stellar ages and the galaxy ages (ΔL and ΔM in Eq. (4)) as a function of galaxy age for a continuous SFH (CONT). The horizontal dotted line reports 0 dex variations. On top of this line are reported the average ⟨Δ⟩L (black text) and ⟨Δ⟩M (red text) for young (7 <  log(t/yr)), intermediate(7 <  log(t/yr) < 9), and old (log(t/yr) > 9) ages, separated by vertical dotted lines.

Variations of light-weighted quantities with respect to the galaxy age are mainly due to the dominant luminosity contribution that stellar populations younger than 10⁸ yr give to the total galaxy luminosity. This is a well-known selection effect (Martinsson et al. 2013; Gomes & Papaderos 2016), in which young populations are more luminous than older populations, biasing the determination of the galaxy age. The new generation of stars systematically exceeds the luminosity of older populations, increasing over time the difference between the galaxy age and the luminosity-weighted mean stellar age (Conroy 2013), as shown in in Fig. 2.

Mass- and light-weighted quantities underestimate the galaxy ages in the first 10⁷ yr when star formation begins and the mass and luminosity of the stellar populations are growing. After 10⁷ yr, the mass-weighted stellar ages evolve steadily, showing a constant offset with respect to the galaxy age of −0.45 dex. To better quantify this trend we divided our simulated galaxy into three different time-steps representing its main evolutionary phases: young (7 <  log(t/yr), intermediate (7 <  log(t/yr) < 9), and old (log(t/yr) > 9) ages (see also Westera et al. 2004).

Figure 2 shows that after 10⁷ yr the mass-weighted stellar ages consistently underestimate the galaxy ages by −0.45 dex (a factor of 3), while for luminosity-weighted quantities the underestimation grows constantly with time, reaching for ages above 1 Gyr an average value of −1.58 dex (a factor of 40). This implies that for galaxies undergoing recent bursts of star formation, because of mergers or self-induced starburst phases, the average stellar population estimated with luminosity-weighted quantities could underestimate the real galaxy age by more than one order of magnitude.

A further consideration regards ⟨log(t_*)⟩_M, which better reproduces the overall galaxy ages, and besides this, it also has a more physical meaning than luminosity-weighted quantities. As shown in the following sections, despite this theoretical advantage, the determination of ⟨log(t_*)⟩_M remains very challenging because of the large variations of mass-to-light ratios as stellar populations evolve (Leitherer et al. 1999; Cid Fernandes et al. 2005; Gomes & Papaderos 2017).

To better characterise REBETIKO galaxies, we use the Hα line emission line EW, EW(Hα), which gives an indication about the contribution of the nebular emission to the total observed continuum. This contribution, produced by ionised gas, can be significant in specific phases of the evolution of a galaxy, for example during starburst phases, producing an overall increase in the EW of the emitted Hα line. Figure 3 reports the evolution of EW(Hα) as a function of galaxy age for CONT and TAU1 SFHs. We are interested in two specific phases of the evolution of galaxies, indicated by the vertical dotted lines in the figure where the galaxy experiences a peak in the production of ionising radiation resulting from the starburst phase.

Fig. 3. Evolution of EW(Hα) as a function of galaxy age for the CONT (solid black line) and TAU1 (dashed red line) models. The solid black and dashed vertical red lines show the ages at which, for each SFH model, the nebular emission contribution is dominant or negligible.

The Lyman continuum production rate, hence the EW(Hα) steeply decrease at these ages, albeit remaining high (200–1000 Å for both star formation scenarios). This phase is identified in Fig. 3 by the regions EW(Hα) > 1000 Å, corresponding to galaxy ages of 5.6 Myr and 3.8 Myr for the CONT and TAU1 models, respectively. After the initial burst, the nebular emission decreases rapidly, depending on the SFH. This is an intermediate phase where nebular emission continues providing an appreciable fraction of the total emission to the case of a continuous SFH; but the intensity of ionising radiation slowly decreases, producing, in turn, weaker emission lines. There is still a rapid jump around 100 Myr, owing to the onset of the post-AGB phase, which produces further gas photo-ionisation (Binette et al. 1994; Stasińska et al. 2008; Cid Fernandes et al. 2011; Gomes et al. 2016). But overall the decline of the nebular emission finally becomes negligible for EW(Hα) < 60 Å, which is the typical value of late-type galaxies in nearby Universe (see e.g., Breda & Papaderos 2018). These phases are fairly different for the two SFHs chosen: in TAU1 this drop in the nebular continuum is immediately after the burst, at 8.8 Myr, while in CONT model it occurs at later stages, around 3.7 Gyr. This implies for such models a presence of nebular contamination for most of the evolution of the galaxies (Leitherer et al. 1999). These timescales must be taken into account when comparing the different results of our analysis in the following sections.

3. Methods: FADO, STARLIGHT, and STECKMAP

We add to the input REBETIKO spectra random Gaussian noise with a varying standard deviation, scaled to a fixed value of the S/N. This procedure assures that each spectrum has a constant S/N throughout the entire range of wavelengths considered in the analysis. In this way we obtained a set of spectra with constant S/N values of 3, 5, 10, 20, 50, and 100, reproducing the range covered by SDSS DR7 (e.g., Fig. 2 in Brinchmann et al. 2004)

The spectral basis for the analysis is built from the BaseL of Bruzual & Charlot 2003, considering 100 SSPs for 25 ages (t = 0.001, 0.00209, 0.00316, 0.00501, 0.00661, 0.00871, 0.01, 0.01445, 0.02512, 0.04, 0.055, 0.10152, 0.1609, 0.28612, 0.5088, 0.90479, 1.27805, 1.434, 2.5, 4.25, 6.25, 7.5, 10, 13, and 15 Gyr) and four metallicities (Z = 0.004, 0.008, 0.02 and 0.05).

We selected three widely used spectral fitting tools for the analysis, which are representative of different approaches to extract evolutionary parameters from a galaxy spectrum: FADO (Gomes & Papaderos 2017), which includes in its analysis nebular emission; the purely stellar code STARLIGHT (Cid Fernandes et al. 2005); and STECKMAP (Ocvirk et al. 2006a,b), which in many aspects can be considered a ‘hybrid’ approach (see Sect. 3.3). In the following section we describe the most important aspects of each code, which are summarised in Table 1.

3.1. FADO

The FADO (Gomes & Papaderos 2017) population synthesis code allows for simultaneous modelling of stellar and nebular emission in a self-consistent framework. It is currently the only publicly available tool permitting determination of SFH and CEM that consistently accounts for the nebular characteristics (Hα, Hβ luminosity and EW). This approach reduces the biases in the determination of mass and light-weighted stellar ages, as shown in Cardoso et al. (2019), and improves the accuracy of the SFH estimation. The code employs a differential genetic optimisation approach called differential evolution optimisation (DEO, Storn & Price 1996, 1997), which assures convergence at an affordable expense of computational time.

The nebular continuum is estimated on the fly, together with electron density and temperature. This separation between stellar and nebular continuum is also maintained for the determination of the intrinsic extinction, which is followed by an automated spectroscopic classification based on diagnostic emission-line ratios. Important features of the code are the following: (a) the inclusion of the nebular continuum to the fit, which can have a significant effect on stellar mass determination (Gomes & Papaderos 2017; Cardoso et al. 2019; Yuan et al. 2019; b) the automatic masking of emission lines and spurious features (e.g., residuals from the sky line subtraction); (c) a pre-defined extinction law, which is set as input parameter; (d) a Gaussian kernel considered for broadening effects, both instrumental and physical.

3.2. STARLIGHT

The STARLIGHT code (Cid Fernandes et al. 2005) fits the observed spectrum with a linear combination of SSPs, modelling the extinction with a parametric law and assuming a Gaussian distribution for the line-of-sight stellar motions. Stellar continuum is normalised to an anchor point (usually set at 4020 Å) varying its shape through a user-defined extinction law.

The minimisation is done through simulated annealing plus the Metropolis scheme, masking regions around emission lines, and bad pixels (MacKay 2003). The Metropolis exploration of the parameter space identifies regions of convergence within the solutions space and after an iteration, the step-size of each parameter is changed accordingly. At the same time, the number of steps of the parameters varies by an amount that scales with the step-size variation. The overall effect of such iterations is a gradual focus of the solutions towards the most likely region in the parameter space, thereby avoiding the problem of local minima.

3.3. STECKMAP

The STECKMAP (Ocvirk et al. 2006a,b) is a tool for interpreting galaxy spectra in terms of combinations of different stellar populations. The high level of degeneracy for the solutions is tackled by regularising the problem through different penalty functions in a Bayesian framework. The main advantages of this approach are the following: (a) the method is non-parametric and does not assume any specific shape for the SFH and CEM; (b) stellar continuum is normalised to have an average of one, and a user defined extinction law is applied on the fly during the minimisation steps; (c) the degeneracy problem is partially solved through a proper regularisation, calibrated with high-resolution spectra; and (d) the minimal request for the solutions is the continuity and an overall low gradient within the solutions space.

At the beginning of this section, we referred to STECKMAP as an example of a hybrid method. This definition originates from the idea that STECKMAP cannot be fully considered as a pure population synthesis code since the matrix of the solution is regularised through an ad hoc penalty function, which forces the solution towards a smooth shape, suppressing large gradients. The a priori exclusion of specific solutions is per se incompatible with the concept of population synthesis, which relies totally on the algorithm chosen for the fitting procedure. To be exact, all the methods considered in this work have some sort of a priori assumption in the sense that all the algorithms exclude negative contributions. This can be considered a ‘pseudo’ penalisation, but it is more related to the mathematical formalism of the problem than the physical motivation behind it. For this reason, we consider STECKMAP to be a hybrid method because it approaches the problem with a minimisation procedure similar to many spectral synthesis methods, but reserving the right to remove valid solutions on the basis of an underlying assumption about the continuity of the solution.

4. Results

The synthetic spectra built with REBETIKO in Sect. 2 were analysed with the tools described in Sect. 3, for six S/N values in the range 3 ≤ S/N ≤ 100. For each SFH we estimated the mean stellar age and metallicity, weighted by mass or light, as defined in Eqs. (2) and (3). Table 2 lists the median of the logarithmic difference between the values recovered by each tool and the REBETIKO input parameter. The results for CONT and TAU1 models are discussed separately in the following subsections.

4.1. Continuous star formation: CONT models

Figure 4 shows the difference between ⟨log(t_*)⟩_L, ⟨log(t_*)⟩_M, ⟨log(Z_*)⟩_L, and ⟨log(Z_*)⟩_M for CONT models with respect to the input REBETIKO value, whose median differences are reported in Table 2. All tools show large differences of the fit at S/N = 3; FADO and STARLIGHT have a median Δ ⟨log(t_*)⟩_L ∼ 0.1 dex, while STECKMAP have a higher Δ ⟨log(t_*)⟩_L of 0.2. All tools have some problems at ages log(t/yr) < 8, where the results have a large spread as a consequence of the instability of the solutions. The variation in the solutions at those ages and S/N values have their origin in the large error bars of spectra, which enormously increases the number of models able to fit the data. This in turn produces almost flat χ² maps, where spectral fitting algorithms could easily remain stuck in a local minimum. From Table 2 we see that FADO and STARLIGHT underestimate ⟨log(t_*)⟩_L by ∼0.1 dex, and in FADO this value is driven mainly by the large errors at log(t/yr) < 8, where results vary by ∼0.3–0.4 dex (∼50%). For older ages, the results become more stable and have a mild tendency to underestimate ⟨log(t_*)⟩_L. For S/N = 3 STECKMAP shows the highest discrepancies, have a median Δ ⟨log(t_*)⟩_L ∼ 0.2 dex, and large differences in the results.

Fig. 4. Results of CONT model analysis. Luminosity- and mass-weighted mean stellar age (first and second row), and metallicity (third and fourth row) residuals for FADO (blue line), STARLIGHT (red line), STECKMAP (green line) as a function of the galaxy age. Luminosity-weighted quantities for FADO and STARLIGHT are obtained considering the best-fitting populations vector at 4020 Å. On top of each panel the S/N of the input REBETIKO spectrum are reported.

At S/N = 5 the solutions span a similar range in terms of Δ ⟨log(t_*)⟩_L, varying between 0.1 and 0.2 dex. The uncertainties of the fits at galaxy ages log(t/yr) < 7.5 are still present; FADO and STARLIGHT show some spread. However, the median Δ ⟨log(t_*)⟩_L for FADO stays within the 0.1 dex because of the higher reliability of the results at older galaxy ages. This trend is confirmed for spectra with S/N >  20, where the uncertainties at log(t/yr) < 8 disappear, decreasing the Δ ⟨log(t_*)⟩_L ∼ 0.04 dex and keeping results differences under the 10%.

For S/N >  20 the median differences between the input REBETIKO ⟨log(t_*)⟩_L and the best fits from FADO are of the order of 0.03 dex (∼7%), a factor 3 and 4 lower than the 0.08 dex (∼20%), and 0.12 dex (∼30%) obtained from STARLIGHT and STECKMAP, respectively. In contrast to FADO, STARLIGHT shows a systematic underestimation of the input ⟨log(t_*)⟩_L by 0.08 dex, 20%, for S/N >  10. The fits show small variations, and the Δ ⟨log(t_*)⟩_L at ages above log(t/yr)∼8 decreases gradually, in correspondence with the predominance of the post-AGB stars. After this phase, we observe higher stability of the fit at log(t/yr) > 9. An interesting aspect to note are the systematic uncertainties in determining the mean stellar ages at log(t/yr)∼7, that is epochs with relevant nebular contribution to the overall continuum (see Fig. 3).

The STECKMAP tool shows median differences of the order of 0.15 dex (∼40%) at all S/N, as shown in Table 2. On average the results are very stable between 7 <  log(t/yr) < 9 yr, with underestimation of ⟨log(t_*)⟩_L ∼0.08 dex at log(t/yr)≤8, similarly to that seen for STARLIGHT, and an overestimation of ∼20% for ages log(t/yr) > 9, owing to an effect already mentioned in Cardoso et al. (2019), and Ocvirk (2010). Both of these papers point towards the sparse coverage of the spectral library used for the analysis at ages 6 <  log(t/yr) < 7, which prevents a proper fit for such epochs and produces, in turn, a decline of the luminosity-weighted stellar age distribution (SAD) at these time-steps. This lack of SSPs at those ages is compensated by an increased fraction of old stellar populations which, having a higher mass-to-light ratio, increases the value of the total mean ⟨log(t_*)⟩_L. This is seen clearly in Fig. 5, in which we show the luminosity-weighted SAD as a function of galaxy age for the three different models, at 2.4, 8.8, and 13.2 Gyr. The vertical lines report the value of ⟨log(t_*)⟩_L recovered by the different tools and the input value from REBETIKO. As we can see for each of the models shown, there is a prominent contribution in light from the first time bin at 1 Myr, followed by lower contributions from SSPs with ages between 6 < log(t/yr) < 7. Then the shape of the SAD returns to the expected trend, with continuous star formation, but to compensate the missing light at 6 <  log(t/yr) < 7, the fitting algorithm during the minimisation increases the overall contribution of older SSPs. The final result of this process is the overestimation of ⟨log(t_*)⟩_L for those galaxy ages (see Fig. 5).

Fig. 5. Luminosity-weighted SAD vs. SSPs age for CONT models and different REBETIKO models: 2.4 (top), 8.8 (middle), and 13.2 (bottom) Gyr. The vertical lines represent the galaxy age from the model (magenta), the input REBETIKO value of ⟨log(t*)⟩L (black), and the values recovered with STECKMAP (green), FADO (blue), and STARLIGHT (red). The top left of each panel indicate the χ2 of the STECKMAP best-fit model.

The third row panels of Fig. 4 show the difference of luminosity-weighted mean stellar metallicity ⟨log(Z_*)⟩_L at different S/N with respect to the REBETIKO input parameter. Broadly speaking, all tools converge towards the input solar metallicity, once input spectra have S/N above 10 and galaxy models have ages above log(t/yr) > 8. The results of FADO show small differences at S/N >  5 on the order of 0.1 dex; there is a mild underestimation of the input value. On the other hand, STARLIGHT shows some discrepancies, Δ ⟨log(Z_*)⟩_L ∼0.2 dex, for ages below log(t/yr) < 8, with a similar pattern at S/N = 50 and 100. The STECKMAP analysis tool shows the largest scatter in the results at galaxy ages log(t/yr) < 8; there is a mild overestimation of the metallicity at higher S/N and old ages, in any case, lower than 0.01 dex.

Mass-weighted quantities show fits that are more noisy with respect to light-weighted fits, as shown in Fig. 4 and Table 2. The FADO uncertainties remain around ∼0.09 dex, 25%, at S/N ≥ 10, which is slightly larger than the 0.04 dex seen for ⟨log(t_*)⟩_L. An interesting feature to note is the oscillation of the results for ages above 4 Gyr, which corresponds to the epochs where the contribution of nebular emission to the total continuum becomes negligible (see Fig. 3).

The STARLIGHT results show median discrepancies between 0.12–0.24 dex and there is no real trend with increasing S/N. These values are driven by the large uncertainties shown at all S/N for ages younger than log(t/yr) < 8 because of the combined effect of rapid variations of mass-to-light ratios at such ages, together with a large contribution of ionised gas to the observed continuum. This trend is confirmed in STECKMAP, which shows similar Δ ⟨log(t_*)⟩_M for log(t/yr) < 8 at all S/N. With respect to STARLIGHT, these discrepancies decrease strongly with S/N, and uncertainties decrease from 0.23 to 0.1 dex for S/N = 10 and S/N = 50; this feature is discussed in detail in Sect. 5. We confirm an overestimation of the ⟨log(t_*)⟩_M for ages above 3 Gyr, as seen already in Cardoso et al. (2019). These variations are of the order of 0.2–0.3 dex for STARLIGHT, but decrease at 0.1 dex in FADO. For such ages, STECKMAP shows more consistent results, where Δ⟨log(t_*)⟩_M <  0.1 dex.

Also mass-weighted mean stellar metallicities show for all codes more unstable results with respect to light-weighted quantities. In particular, all the tools have difficulties recovering the metallicity evolution at S/N <  20. In FADO we see more stable fits at higher S/N, with mild differences at ages above log(t/yr) > 9, a feature seen also in Cardoso et al. (2019). The STARLIGHT best fits systematically underestimate the metallicity by ∼0.01 dex, in particular at young ages, and this trend rises to 0.6 dex for log(t/yr) < 8, increasing with the S/N. In this context, STECKMAP shows more stable results, with a mild underestimation of the metallicity at older ages, similarly to the light-weighted quantities. However, for this tool, we have to consider that the luminosity-weighted quantities are converted into mass-weighted quantities considering a conversion factor that takes into account the metallicity evolution obtained from the fit, normalised by construction. This implies that the differences in the metallicity are of the same order of magnitude and the effect of the mass-to-light ratio overall is negligible (see Eq. (3) in Ocvirk et al. 2006b).

4.2. Instantaneous burst: TAU1 models

Figure 6 and Table 2 report the same analysis as the previous subsection, considering SFH for the TAU1 models. We have already seen in Sect. 2 that the nebular contribution for this SFH models declines very rapidly, and for this reason, we expect higher reliabilities from pure stellar synthesis codes. This is confirmed in Fig. 6, which shows an overall consistency of the results; the median Δ ⟨log(t_*)⟩_L is approximately 2-3% for FADO and STARLIGHT already at S/N = 5. Slight oscillations are seen in STECKMAP at ages ∼100–400 Myr, as a consequence of the post-AGB stars phase, which adds further degeneracy to the fit procedure, but overall the results are in good agreement with the input REBETIKO parameters. A feature in common between FADO and STARLIGHT is a small bump at ages between 3–5 Gyr, presumably due to degeneracies between old SSPs in the adopted SSP library (see Cardoso et al. 2019). The results of STECKMAP are dependent on the input S/N; a reasonable accuracy is reached at S/N >  10. For lower S/Ns the fits are more unstable and we observe an overestimation of ⟨log(t_*)⟩_L at log(t/yr) > 7.5, whose origins are clarified in Sect. 5.

Fig. 6. Same caption of Fig. 4 for TAU1 SFH models.

The results for mean metallicities are more homogeneous than CONT models between the different tools. The solutions are more unstable at log(t/yr)∼7; FADO and STARLIGHT show a peak of ∼0.2 dex at ∼200 Myr in ⟨log(Z_*)⟩_L and ⟨log(Z_*)⟩_M resulting from the poor coverage of the spectral basis at those ages (see Cardoso et al. 2019). Curiously, STECKMAP has similar problems when fitting mean metallicity at ages around log(t/yr)∼8, but the sign of the deviation from the true value goes in the opposite direction with respect to FADO and STARLIGHT, where STECKMAP underestimates the REBETIKO values for low S/N and the other two codes showing an overestimation. In contrast to FADO and STARLIGHT, STECKMAP fits show high instabilities at older ages above log(t/yr) > 9.

Analogous to the results obtained for CONT models, mass-weighted quantities show a larger scatter, as shown in Fig. 6 and Table 2. Overall, the results of FADO fits have smaller differences, showing similar Δ ⟨log(t_*)⟩_L and Δ ⟨log(t_*)⟩_M, ∼0.1−0.2 dex, while for STARLIGHT and STECKMAP the average uncertainties increase by a factor of ∼3. Both methods show large uncertainties at ages below log(t/yr) < 8, where the nebular contribution is still significant. For STECKMAP the scatter decreases at S/N >  10, while for FADO and STARLIGHT the fits are more stable for S/N ≥ 5.

5. Discussion

As we have seen in Sect. 4, the reliability of the results for each code is strongly dependent not only on the input SFH but also on other factors: the S/N; the contribution of the nebular continuum, particularly for pure stellar synthesis codes; and the age coverage of the spectral basis. Regarding stellar ages, the accuracy of the estimated ⟨log(t_*)⟩_L and ⟨log(t_*)⟩_M changes drastically according to the input SFH. For CONT models, Table 2 shows that for FADO the uncertainties increase by a factor between 3–5 from the use of light- and mass-weighted mean stellar ages. For TAU1 models there are no relevant differences between light- and mass-weighted stellar ages. In STARLIGHT the uncertainties are about a factor 2–4, for both CONT and TAU1 models. The average accuracy of 0.15 dex for ⟨log(t_*)⟩_M seen for STARLIGHT confirms the results of a previous analysis of Cid Fernandes et al. (2005), and are typical of purely stellar codes (Ocvirk et al. 2006a,b; Tojeiro et al. 2007; Koleva et al. 2009; Wilkinson et al. 2017; Cardoso et al. 2019).

Larger discrepancies are associated with evolutionary phases where the nebular contribution to the continuum is not negligible, which is particularly evident in CONT spectra. This confirms the importance of considering nebular emission for starburst systems, where this component can reach up to 50% of the total optical and near-infrared emission (Krueger et al. 1995; Izotov et al. 1997; Papaderos et al. 1998; Leitherer et al. 1999; Schaerer & de Barros 2009).

As discussed in Sect. 2, the parameters defined in Eq. (2) to quantify the assembly history of a galaxy must be used with caution. By definition, these are averaged quantities that miss important details of galaxy evolution, such as starburst episodes or sudden shutdown of the star formation, and they should be considered as the first moment of the distribution (Cid Fernandes et al. 2005). These details, however, can be investigated more thoroughly by analysing the luminosity-weighted SAD and the SFH recovered from the fitting process. Such analysis allows for a more complete insight into the fitting procedure, providing hints about the observed discrepancies and helping to interpret the results.

For the mean stellar metallicities defined in Eq. (3) the same point is valid, even if there is a further warning to consider. The synthetic spectra built in Sect. 2 assumed a constant solar metallicity: since the SFHs chosen to model our galaxies are different, the constant metallicity imposed for CONT and TAU1 models is a simplification. For CONT models, a constant SFR implies that the quantity of processed material within the galaxy increases with time, raising, in turn, the value of the average metallicity. To keep constant Z, a mechanism responsible for metals dilution must be introduced, such as a constant cosmic gas accretion from the circum-galactic medium. On the other hand, TAU1 models require the opposite effect: a closed-box model, in which the galaxy does not produce further stellar populations, and the constant metallicity is guaranteed by the absence of any replenishment or inflow of low-metallicity gas from the circum-galactic medium (see Tumlinson et al. 2017, for a more complete discussion) and the absence of new star formation episodes, responsible for ISM metal enrichment. Some additional cautionary remarks are relevant, for example the assumptions of (a) no velocity broadening, (b) constant spectral resolution, (c) constant noise across the considered spectral range, (d) constant electron density and temperature, (e) no Lyman continuum photon leakage throughout the evolution of a galaxy, and (f) no shocks are strong simplifications. Therefore the presented tests do not mean that application of the three codes to real spectra will closely match what is found in this work. However, it is reasonable to assume that in broad terms the biases found in this study, in particular the overestimation of ⟨log(t_*)⟩_M in phases with strong nebular emission, should affect studies of real galaxy spectra.

5.1. Uncertainties introduced by a different spectral basis or a change in the wavelength coverage

An important and complex aspect to consider in our analysis is the impact on the results of a different spectral basis or a change in the wavelength coverage. These issues are beyond the goal of this paper, and it is part of ongoing works within our team (Gomes et al., in prep.). However, in this section, a qualitative analysis gives some hints on possible systematics, reproducing in our tests existing instruments.

A different SSP basis would certainly affect our results, reflecting at a first-order the systematic differences in the stellar basis used to build the SSPs, as shown for example in Koleva et al. (2008). The use of different stellar libraries has a strong effect on the near-infrared regime, while the effect in the optical is negligible, as seen for example in Baldwin et al. (2018; see their Fig. 7). In their work Baldwin et al. (2018) show that the largest effect on the derived SFHs is related to the stellar spectral library chosen for the analysis, impacting more the results than the method considered. To better quantify this effect we compare in Fig. 7 the residuals between the mean stellar ages estimated with STECKMAP and the REBETIKO input obtained with BaseL as a function of the age of the galaxy for the TAU1 models. The results were obtained using as a spectral basis the standard BaseL, defined in Sect. 3 (yellow line), and MILES (Vazdekis et al. 2010, green line).

Fig. 7. Luminosity- and mass-weighted mean stellar age (1st and 2nd row) residuals for STECKMAP, obtained with BaseL (yellow line, Bruzual & Charlot 2003) and MILES (green line, Vazdekis et al. 2010) for an instantaneous burst SFH (TAU1 models). The S/N of the input REBETIKO spectra are indicated on the top of each panel.

The main differences in the light-weighted quantities are seen for ages below 10⁸ yr, where the age coverage of MILES is sparser than BaseL, resulting in larger uncertainties. For older stellar ages MILES fits overestimate the input value and have differences below 0.5 dexes. For mass-weighted stellar ages, the trends reported are qualitatively similar, and the overestimation starts at ∼10⁸ yr for both spectral bases. The overall differences are reported in Table 3, showing a general consistency of the results obtained with the two SSPs basis.

The problem of the spectral basis considered in a fit has been thoroughly addressed in Cardoso et al. (2019), which showed how larger base libraries do not necessarily improve the results. The important part in building a good base is highly dependent on the type of sample in question (e.g., spectral coverage; passive or star-forming galaxies). For instance, a good base could require sufficient young stellar populations <100 Myr in age to properly cover the star-forming regimes, while at the same time not having too many old stellar populations >1–10 Gyr in age because of the expected nature of the objects analysed.

Concerning the wavelength coverage issue, Figs. 8 and 9 show the difference between the estimated stellar properties (mean ages and metallicities) as a function of the age of the galaxy model for continuous (Fig. 8) and instantaneous (Fig. 9) SFHs. The results shown have been obtained by fitting the synthetic spectra from REBETIKO in two spectral ranges: one simulating that of MUSE instrument (Bacon et al. 2010, ∼4800–9300 Å, including Hα+Hβ as well as the Paschen jump), and another in the range 3500–7000 Å, simulating the combined V1200+V500 data from the Potsdam Multi-Aperture Spectrophotometer (PMAS, Kelz et al. 2006), used in the CALIFA survey (Sánchez et al. 2012). The blue and red lines indicate the results obtained with MUSE coverage with FADO and STARLIGHT, while the cyan and orange lines show the best fits for the CALIFA coverage. The plots show results for S/N = 5, 20, and 100, but the trends seen are present in the other S/N values previously discussed, only varying in the amount of noise in the biases.

Fig. 8. Luminosity- and mass-weighted mean stellar age (first and second row), and metallicity (third and fouth row) residuals for FADO (blue and cyan lines) and STARLIGHT (red and orange lines), obtained with MUSE (4800–9300 Å) and CALIFA (3500–7000 Å) spectral coverage (see labels) for a continuous SFH (CONT models). The S/N of the input REBETIKO spectra are indicated on the top of each panel.

Fig. 9. Same caption of Fig. 8 for TAU1 SFH models.

As a general comment, FADO fits show overall a better agreement with the input spectra compared to STARLIGHTin both coverage set-ups. However, each set-up yields new trends in both codes, justifying a more detailed interpretation in a future paper, as already mentioned at the beginning of this section. A bluer coverage in FADO, such as that in CALIFA, yields better light-weighted age estimates than its red counterpart (MUSE) for both SFHs. Interestingly, both FADO and STARLIGHT show an overestimation of stellar ages when MUSE coverage is considered for a continuous SFH with increasing galaxy age, linked to degeneracies between the old SSPs in the base (Fig. 8). This is rather intuitive since the continua of old SSPs (>1 Gyr) do not strongly differ in 6000–8900 Å. A bluer coverage in STARLIGHT also yields better results than a redder coverage.

The results flip when looking at the light-weighted metallicity, and for such cases, a redder coverage gives better results than a bluer coverage in FADO and STARLIGHT for an instantaneous burst SFH (Fig. 9). The same holds for a continuous SFH in STARLIGHT. When looking at the mass-weighted age and metallicity, FADO yields better results when assuming a red coverage; uncertainties in blue are probably linked to a poor age coverage of the adopted base library. The STARLIGHT tool also yields better mass-weighted age estimates when red coverage is available, but this effect is less prominent for an instantaneous burst. Finally, the blue coverage leads to a systematic metallicity overestimation in both codes (especially for older galaxy models), whereas a red coverage leads to a systematic metallicity underestimation in STARLIGHT for decreasing model age.

In summary, FADO seems to be more efficient than STARLIGHT in recovering the mean age and metallicity during strong star-forming evolutionary stages (t >  0.1 Gyr) for MUSE and CALIFA coverages, regardless of S/N and SFH. These results show that the impact of the spectral coverage on this type of analysis is not negligible and this is something we will explore in a future paper (Gomes et al., in prep.). For example, an important caveat to this trend is the systematic overestimation of the light-weighted mean stellar age by FADO in the MUSE set-up after t >  0.1 Gyr and for a continuous SFH, which is likely due to the absence of the Balmer jump at 3646 Å needed to assure self-consistency between the stellar and nebular spectral contributions. This absence does not affect the instantaneous burst SFH (TAU1 models) because after t >  0.1 Gyr there is only a residual nebular spectral contribution coming from post-AGB stars.

5.2. Star formation history from purely stellar codes

A systematic feature seen in STARLIGHT and STECKMAP is the overestimation of ⟨log(t_*)⟩_M for log(t/yr) < 8 (see Figs. 4 and 6). This feature affects both CONT and TAU1 models, and there is a tendency to decrease with increasing S/N. The large errors, up to 3 dex, are associated with galaxy ages for which the contribution of nebular emission to the overall continuum is relevant (Fig. 3) and points towards the origin of these discrepancies.

Figure 10 shows an example of a best-fit spectrum obtained with FADO for a metal-poor blue compact galaxy (BCD CGCG 007-025). The contribution to the overall spectrum due to the nebular continuum is shown as a red line below the spectrum. The peak of this contribution is at wavelengths below 4000 Å, coinciding with the Balmer jump at ∼3646 Å. The contribution of the nebular continuum increases with increasing wavelength, becoming dominant in the region around the Paschen jump at ∼8207 Å. In this part of the spectrum, both contributions are relatively small because of the decreased stellar continuum, which is typical of star-forming galaxies. In the redder part of the spectrum the nebular continuum is dominant, and this has profound implications for pure spectral synthesis codes: as discussed in Izotov et al. (2011) and Papaderos & Östlin (2012) the reddish nebular continuum of starburst galaxies monotonically increases with wavelength and forces purely stellar SED fitting codes to artificially raise the fraction of old stars.

Fig. 10. FADO spectral fit example, from Gomes & Papaderos (2017), for the metal-poor blue compact galaxies (BCD CGCG 007-025, yellow line) observed in the SDSS (see also Guseva et al. 2007). Panel a: best-fitting synthetic SED with stellar and nebular continuum emission, shown in blue, grey and red, respectively. Panel b: residuals between fit and observed spectrum; the shaded area and the dashed curve delineate the ±1σ and ±3σ error spectrum, respectively.

For luminosity-weighted quantities this effect is less relevant, even in the presence of a strong nebular component. Figure 11 shows the luminosity-weighted SAD and the SFH obtained for CONT and TAU1 models with STECKMAP at galaxy ages of 30 Myr. The CONT models show a small bump of old stellar populations around log(t/yr) > 9 (top left panel); but this bump, because of the high mass-to-light ratios of these stellar populations, does not affect the estimation of ⟨log(t_*)⟩_L, which is recovered with an error of ∼0.1 dex. A similar effect is seen for TAU1 models (bottom left panel), where STECKMAP fits report a very small contribution from a stellar population with ages above 10 Gyr; this is barely visible from the figure.

Fig. 11. Luminosity-weighted SAD (left panels) and SFR (right panels) vs. the SSPs ages obtained with STECKMAP for CONT (top) and TAU1 (bottom) models with S/N = 10 at a galaxy age of 30 Myr. The vertical lines are shown as in Fig. 5.

When we consider mass-weighted quantities, the small bumps seen in the SAD of Fig. 11 become the dominant stellar populations (right panels), producing an overall overestimation of ⟨log(t_*)⟩_M, in agreement with the results of (Cardoso et al. 2019, see their Fig. B.1). These differences, as mentioned at the beginning of Sect. 5, depend strongly on the contribution of the nebular emission to the overall continuum. For this reason in the CONT and TAU1 models at 110 Myr, where the nebular emission is lower, the fits consistently recover the main contributors in both ⟨log(t_*)⟩_L and ⟨log(t_*)⟩_M, as shown in Fig. 12. In other words where the nebular contribution is relevant there is a systematic overestimation of ⟨log(t_*)⟩_M, which has led to an overestimation of the stellar masses (see also Gomes & Papaderos 2017 and Cardoso et al. 2019).

Fig. 12. Same caption of Fig. 11 for a galaxy age of 110 Myr.

5.3. The Balmer jump effect

Table 2 and Figs. 4, and 6 show that the median Δ ⟨log(t_*)⟩_M measured with STECKMAP span a relatively extended range of differences from REBETIKO input parameters, for both CONT and TAU1 models: from ∼0.4 dex at S/N = 3 up to 0.01 dex at S/N = 100. A relevant gap is seen between S/N = 50 and S/N = 10, where for example in CONT models Δ⟨log(t_*)⟩_M decreases from 0.23 (∼70%) up to 0.1 (∼25%) dex. We explore the reasons for this decline.

From the fourth row panels of Figs. 4 and 6 it is clear that the improvement for STECKMAP at higher S/N is mainly driven by better fits at log(t/yr) < 8. At those ages, STARLIGHT also shows a large bump of up to 3 dex, similar to that seen with STECKMAP. In this context the FADO results are quite different, showing a small bump at log(t/yr) < 7 for S/N = 5, and overall a good agreement for higher S/N. This indicates that the contamination of nebular emission to the overall continuum is responsible for this scatter. However, since STECKMAP fits improve with increasing S/N even for models in which the nebular contribution decreases relatively quickly, this effect must be related to specific features in the spectra involving other physical mechanisms.

To investigate the origin of such S/N-linked overestimation of ⟨log(t_*)⟩_M in detail, we analyse the SFH recovered from the STECKMAP fits at different ages that are representative of galaxy evolutionary phases with prominent, intermediate, or negligible nebular contributions. Figure 13 shows the SFH recovered with STECKMAP for CONT and TAU1 models at galaxy ages of 4.4, 180, and 320 Myr and S/N 10 and 50, respectively.

Fig. 13. Star formation rate vs. the SSPs ages obtained with STECKMAP for CONT (top panels) and TAU1 (bottom panels) models with S/N = 10 (left column), and 50 (right column) at a galaxy age of 4.4 (top row), 180 (middle row), and 320 (bottom row) Myr. The vertical lines are indicated as in Fig. 5.

For CONT models at S/N = 10 STECKMAP fits have large discrepancies at 4.4 and 180 Myr. The best-fit spectra for the ages of these galaxies assume a large contribution from old stellar populations, with low mass-to-light ratios and a huge weight in term of total star formation. The fit improves substantially at 320 Myr, where the contribution of the different SSPs is well modelled, as Fig. 4 confirms. When increasing the S/N to 50 the results at 320 Myr remains almost unchanged, but the ⟨log(t_*)⟩_M estimated at 4.4 and 180 Myr improves notably.

For TAU1 models we observe a similar trend, but since the nebular contribution decreases more quickly for such SFHs, this effect is more evident only at 4.4 Myr, and ⟨log(t_*)⟩_M is positively recovered at 180 and 320 Myr even for S/N = 10. For this reason for the rest of the analysis, we focus on the CONT models.

Increasing the S/N, STECKMAP improves the quality of the fits, reproducing more the input spectrum more accurately. Figure 14 shows the best-fit spectra obtained with STECKMAP for the ages investigated in Fig. 13, 4.4, 180, and 320 Myr, and focussing on the region of the spectrum between 3300–4500 Å. This was done to underline an important effect visible regarding discontinuities, in particular around the Balmer jump, at 3646 Å.

Fig. 14. STECKMAP best-fit spectrum for CONT models at 4.4, 180, and 320 Myr and S/N = 10 (top panels), and 50 (bottom panels). The black and red lines show the input and the best-fit spectrum, respectively. The blue, magenta, and green lines show the mask, the residual and the errors, respectively.

As seen in Fig. 10, the nebular emission, according to the evolutionary stage of the galaxy, can contribute conspicuously to the shape and amplitude of the Balmer jump. The shape and the strengths of the Balmer jump are set by the evolution of the number of ionising photons, which evolves according to the stellar ages of the SSPs contributing to the observed spectrum. This is not only a theoretical expectation but it has been observed in galaxies with strong starburst activity (e.g., Guseva et al. 2007).

Figure 14 shows that one of the causes of the high ⟨log(t_*)⟩_M recovered with purely stellar codes is related to the difficulty of such methods in reproducing the Balmer jump (see Fig. 1). To quantify this effect better, we defined a new parameter for each REBETIKO synthetic spectrum built for CONT models called the Balmer jump strength (BJS), which is written as

$$\begin{array}{c}\hfill \mathrm{BJS}=\frac{{\int }_{3580\mathrm{\AA }}^{3600\mathrm{\AA }}F(\lambda ,t)\mathrm{d}\lambda }{{\int }_{3660\mathrm{\AA }}^{3680\mathrm{\AA }}F(\lambda ,t)\mathrm{d}\lambda }.\end{array}$$(5)

The evolution of the BJS can be considered as a diagnostic of the ‘height’ of the Balmer jump, and varies according to the evolutionary phase of the galaxy, as shown for different S/N and CONT models in Fig. 15. The effect of S/N is rather small, as expected, and its values span a range between 0.8 < BJS < 1.8, with an average of ⟨BJS⟩ = 1.06 ± 0.11. Within these REBETIKO synthetic spectra we defined three subgroups (indicated by the horizontal lines in Fig. 15), according to the BJS value as follows:

• High BJS: 1.2 < BJS < 2

• Intermediate BJS: 1.08 < BJS < 1.12

• Low BJS: 0.92 < BJS < 0.98.

Fig. 15. Evolution of the BJS (see Eq. (5)) as a function of galaxy ages for CONT models at S/N = 5 (black), 10 (yellow), 50 (magenta), and 250 (cyan). Horizontal lines show the subgroup identifying ages with weak (0.92 < BJS < 0.98, red line), intermediate (1.08 < BJS < 1.12, green line) and strong (1.2 < BJS < 2, blue line) BJS.

These subsamples define families of spectra with different strengths in the Balmer jump, representatives of different evolutionary phases of a galaxy. The next step is to investigate the correlation of the Δ ⟨log(t_*)⟩_M analysed in Sect. 4 with the BJS.

Figures 4 and 16 show the residual of ⟨log(t_*)⟩_M as a function of galaxy age, from analysis restricted to the STECKMAP results and flagging the synthetic spectra with high (blue), intermediate (green), and low (red) BJS. According to the selection criteria defined above, we can identify three different cases:

• Case 1: The spectrum has a high BJS because of the relative nebular continuum contribution to the overall SED. For these cases, purely stellar codes have immense difficulties in determining mean stellar ages independently of the quality of data, as shown by blue circles in Fig. 16. We see from the two lower panels of Figs. 4 and 6 that the metallicity inferred by STECKMAP is strongly sub-solar (off-scale in some panels). A reasonable explanation for this is that purely stellar codes, such as STECKMAP and STARLIGHT, have no other option for fitting the sudden increase in flux at 3650 Å than to assume a population of very hot, thus low-metallicity, massive stars.

• Case 2: The spectra has an intermediate BJS, with a moderate contribution from nebular emission. In general, while the galaxy evolves, this component decreases and the fits converge to the input ⟨log(t_*)⟩_M. However, for specific galaxy ages, purely stellar codes are not able to discriminate the correct stellar ages, unless the Balmer jump in these phases is small enough so that the uncertainties are dominated by the S/N and/or the intrinsic degeneracies between SSP, instead of the nebular continuum, as happens in the cases represented by the green triangles in Fig. 16. The SSPs at these ages change in short timescales, as shown in Fig. 3. The EW of Hα in these spectra shows rapid jumps around 100 Myr because of the onset of the post-AGB phase, which produces further gas photo-ionisation (Binette et al. 1994; Stasińska et al. 2008; Cid Fernandes et al. 2011; Gomes et al. 2016). The Δ ⟨log(t_*)⟩_M at ages higher than log(t/yr)∼8 decreases gradually, depending on the importance of the post-AGB stars, adding further degeneracy to the fit procedure.

• Case 3: The spectra has a low BJS, implying low contamination from nebular emission of the overall continuum or, alternatively, for evolved galaxies in which old stellar populations dominate the spectrum. For these evolutionary phases, purely stellar codes recover the ⟨log(t_*)⟩_M with an accuracy better than 0.4 dex at S/N ≥ 5, as shown with red stars in Fig. 16.

Fig. 16. Δ⟨log(t*)⟩M obtained from STECKMAP vs. galaxy age at different S/N, reported on the top right of each panel. The blue circle, green triangles, and red stars point out the synthetic REBETIKO spectrum with high, intermediate, and low BJS, respectively.

These considerations imply that when a galaxy enters into an evolutionary phase in which the observed spectrum has a strong BJS, the effect of neglecting the nebular emission in the fitting process has a strong impact on the estimation of the SFH at such ages. But even in the presence of mild contamination by nebula continuum, there is still a possibility to misinterpret the data, as a consequence of its poor quality. This can produce an overestimation of the stellar ages and a simultaneous underestimation of stellar metallicity, since the intermediate Balmer jump height cannot be fit by these purely stellar codes, introducing a higher contribution of old, low mass-to-light ratio SSPs to compensate for this effect. In our mock experiment, we identified this transition phase in the range 10 <  S/N <  50, which is also the transition where purely stellar codes showed the largest differences in the results, as a consequence of unstable fits.

6. Conclusions

In this paper we investigated the limits and reliability of three spectral synthesis codes for the estimation of two specific parameters representative of the evolution of a galaxy: the mean stellar age and the mean stellar metallicity. The question addressed in the introduction was at which S/N values is it still possible to determine from a galaxy spectrum its star formation and chemical evolution history and how does this depend on the tool used to model the spectra.

To address this question we built a set of synthetic spectra with both stellar and nebular emission, thereby reproducing the evolution of a galaxy with two different SFHs: a continuous star formation (CONT) and a single initial burst with an exponentially declining star formation law (TAU1). The synthetic spectra were modified for a range of S/N ratios and were analysed with the three spectral synthesis codes: FADO, STECKMAP, and STARLIGHT. The main results from our analysis are as follows:

• At S/N = 3 all tools show large differences in the results; FADO and STARLIGHT have a median Δ ⟨log(t_*)⟩_L ∼0.1 dex and STECKMAP has a higher Δ ⟨log(t_*)⟩_L of 0.2. The codes have a large spread of ⟨log(t_*)⟩_L at ages log(t/yr) < 8, as a consequence of the instability of the solutions.

• For S/N >  20 the median differences between the input REBETIKO ⟨log(t_*)⟩_L and the best fits from FADO are of the order of 0.03 dex (∼7%), which is a factor 3 and 4 lower than the 0.08 dex (∼20%) and 0.11 dex (∼30%) obtained from STARLIGHT and STECKMAP, respectively. In contrast to FADO, STARLIGHT shows a systematic underestimation of the input ⟨log(t_*)⟩_L = 0.08 dex, 20%, for S/N >  10.

• Mass-weighted quantities show results that are less reliable than light-weighted results (Fig. 4), with differences around 0.08 dex, 20%, at S/N ≥ 10 for FADO. The STARLIGHT results show median discrepancies between 0.13–0.19 dex; there is no real trend with increasing S/N.

• Since the nebular contribution for TAU1 models declines very rapidly, the results for this SFH shows an overall consistency, with median Δ ⟨log(t_*)⟩_L around 5% for FADO and STARLIGHT already at S/N = 5 (Fig. 6).

The impact of a different spectral basis and/or wavelength coverage on the result is very complex and is postponed to a subsequent paper of this series (Gomes et al., in prep.). Our preliminary analysis has shown that a different spectral basis affects our results mainly as a consequence of the systematic differences in the stellar basis used to build the SSPs. Comparing the standard BaseL basis defined in Sect. 3 with MILES (Vazdekis et al. 2010), we see that the main differences in both light- and mass-weighted quantities arise for ages below 10⁸ yr, where the age coverage of MILES is sparser than BaseL. This confirms the importance of building an appropriate spectral basis that is calibrated for the type of galaxies investigated.

The variation of wavelength coverage affects mean ages and metallicities in opposite ways: To estimate mean stellar ages a bluer coverage yields better estimates than its red counterpart, where FADO shows overall a better agreement compared to STARLIGHT. For the mean metallicity, the trend goes in the opposite direction, and a redder coverage gives better results than a bluer coverage with both FADO and STARLIGHT.

Larger discrepancies in the results are associated with evolutionary phases in which the nebular contribution to the overall SED is significant, which is particularly evident in CONT spectra. Although less significant for its light-weighted counterparts, the presence of a strong nebular component leads to significant biases in the estimated mass-weighted quantities, mainly owing to a severe overestimation of old SSPs in pure-stellar fitting approaches (Fig. 11). This implies that if the contribution from nebular continuum is relevant, there is a systematic overestimation of ⟨log(t_*)⟩_M, which has led to an overestimation of the stellar masses up to 50% for t >  300 Myr (see also Gomes & Papaderos 2017; Cardoso et al. 2019).

Detailed investigations of the best-fit spectrum in galaxies with overestimated ⟨log(t_*)⟩_M allowed the determination of the origin of such discrepancies: the difficulties of such methods to reproduce the Balmer jump at 3646 Å. Quantifying this effect through the definition of the BJS parameter (Eq. (5)), we identified three possible cases, discriminated by the presence of high, intermediate, or weak BJS (Fig. 16):

• Case 1: For galaxies with high BJS purely stellar codes have major difficulties in determining mean stellar ages, irrespective of the data quality.

• Case 2: For galaxies with intermediate BJS pure stellar codes are not able to discriminate the correct stellar ages, unless the spectrum has a sufficiently high S/N, from our test at least S/N >  10.

• Case 3: For galaxies with weak BJS purely stellar codes recover the ⟨log(t_*)⟩_M with an accuracy < 0.4 dex, at S/N = 5, thereby providing consistent results.

These considerations imply that when a galaxy enters into an evolutionary phase in which the observed spectrum has a high BJS, the effect of neglecting the nebular emission in the fitting process has a strong impact on the estimation of the SFH of such galaxy. But even in the presence of a mild contribution from nebular continuum, there is still the possibility to misinterpret the data as a consequence of the poor quality of the observations. This can produce an overestimation of the stellar ages, for example Fig. 4, since intermediate Balmer jumps cannot be fit by a purely stellar code, which models the SFH of the spectrum adding an higher contribution from old SSPs with low mass-to-light ratios (Cardoso et al. 2019).

Our work underlines the importance of using appropriate codes when considering the spectra from galaxies with high SFRs. The introduction of a self-consistent treatment of nebular emission in this context assumes a fundamental role, being the only viable route towards a reliable determination of the assembly history of galaxies in starburst phase, which are more common at higher redshifts.

Acknowledgments

The authors would like to thank the anonymous referee for the constructive feedback on this manuscript, which helped to improve the content and the presentation of the article. This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through the research grants PTDC/FIS-AST/29245/2017, UID/FIS/04434/2019, UIDB/04434/2020 and UIDP/04434/2020. C. P. acknowledges support from DL 57/2016 (P2460) from the ‘Departamento de Física, Faculdade de Ciências da Universidade de Lisboa’ L.S.M.C. acknowledges support by the project “Enabling Green E-science for the SKA Research Infrastructure (ENGAGE SKA)” (reference POCI-01-0145-FEDER-022217), funded by COMPETE2020 and FCT. P.P. acknowledges support by the project “Identifying the Earliest Supermassive Black Holes with ALMA (IdEaS with ALMA)” (PTDC/FIS-AST/29245/2017). T.S. acknowledges support from DL 57/2016/CP1364/CT0009 from the ‘Centro de Astrofísica da Universidade do Porto, Portugal’. R.C. acknowledges support from the Fundação para a Ciência e a Tecnologia (FCT) through the Fellowship PD/BD/150455/2019 (PhD:SPACE Doctoral Network PD/00040/2012) and POCH/FSE (EC). A.P.A. acknowledge support by FCT under Project CRISP PTDC/FIS-AST-31546/2017 and UIDB/00099/2020. We would like to thank D. Munro for freely distributing his Yorick programming language (available at http://www.maumae.net/yorick/doc/index.html). C. P. acknowledges the Luisinho family, for the support offered during pandemic times in their house, and the countryside of Olho Marinho (Obidos), where most of this manuscript has been thought and written.

References

All Tables

All Figures

	Fig. 1. Synthetic spectra normalised at λ = 4020 Å for TAU1 (instantaneous burst, top panel) and CONT (continuous, bottom panel) SFH models. Figures adapted from Figs. 1 and. 2 of Cardoso et al. 2019.
In the text

Fig. 2. Differences between luminosity-weighted (solid black line) and mass-weighted (dashed red line) mean stellar ages and the galaxy ages (ΔL and ΔM in Eq. (4)) as a function of galaxy age for a continuous SFH (CONT). The horizontal dotted line reports 0 dex variations. On top of this line are reported the average ⟨Δ⟩L (black text) and ⟨Δ⟩M (red text) for young (7 <  log(t/yr)), intermediate(7 <  log(t/yr) < 9), and old (log(t/yr) > 9) ages, separated by vertical dotted lines.

In the text

	Fig. 3. Evolution of EW(Hα) as a function of galaxy age for the CONT (solid black line) and TAU1 (dashed red line) models. The solid black and dashed vertical red lines show the ages at which, for each SFH model, the nebular emission contribution is dominant or negligible.
In the text

Fig. 4. Results of CONT model analysis. Luminosity- and mass-weighted mean stellar age (first and second row), and metallicity (third and fourth row) residuals for FADO (blue line), STARLIGHT (red line), STECKMAP (green line) as a function of the galaxy age. Luminosity-weighted quantities for FADO and STARLIGHT are obtained considering the best-fitting populations vector at 4020 Å. On top of each panel the S/N of the input REBETIKO spectrum are reported.

In the text

Fig. 5. Luminosity-weighted SAD vs. SSPs age for CONT models and different REBETIKO models: 2.4 (top), 8.8 (middle), and 13.2 (bottom) Gyr. The vertical lines represent the galaxy age from the model (magenta), the input REBETIKO value of ⟨log(t*)⟩L (black), and the values recovered with STECKMAP (green), FADO (blue), and STARLIGHT (red). The top left of each panel indicate the χ2 of the STECKMAP best-fit model.

In the text

	Fig. 6. Same caption of Fig. 4 for TAU1 SFH models.
In the text

	Fig. 7. Luminosity- and mass-weighted mean stellar age (1st and 2nd row) residuals for STECKMAP, obtained with BaseL (yellow line, Bruzual & Charlot 2003) and MILES (green line, Vazdekis et al. 2010) for an instantaneous burst SFH (TAU1 models). The S/N of the input REBETIKO spectra are indicated on the top of each panel.
In the text

	Fig. 8. Luminosity- and mass-weighted mean stellar age (first and second row), and metallicity (third and fouth row) residuals for FADO (blue and cyan lines) and STARLIGHT (red and orange lines), obtained with MUSE (4800–9300 Å) and CALIFA (3500–7000 Å) spectral coverage (see labels) for a continuous SFH (CONT models). The S/N of the input REBETIKO spectra are indicated on the top of each panel.
In the text

	Fig. 9. Same caption of Fig. 8 for TAU1 SFH models.
In the text

Fig. 10. FADO spectral fit example, from Gomes & Papaderos (2017), for the metal-poor blue compact galaxies (BCD CGCG 007-025, yellow line) observed in the SDSS (see also Guseva et al. 2007). Panel a: best-fitting synthetic SED with stellar and nebular continuum emission, shown in blue, grey and red, respectively. Panel b: residuals between fit and observed spectrum; the shaded area and the dashed curve delineate the ±1σ and ±3σ error spectrum, respectively.

In the text

	Fig. 11. Luminosity-weighted SAD (left panels) and SFR (right panels) vs. the SSPs ages obtained with STECKMAP for CONT (top) and TAU1 (bottom) models with S/N = 10 at a galaxy age of 30 Myr. The vertical lines are shown as in Fig. 5.
In the text

	Fig. 12. Same caption of Fig. 11 for a galaxy age of 110 Myr.
In the text

	Fig. 13. Star formation rate vs. the SSPs ages obtained with STECKMAP for CONT (top panels) and TAU1 (bottom panels) models with S/N = 10 (left column), and 50 (right column) at a galaxy age of 4.4 (top row), 180 (middle row), and 320 (bottom row) Myr. The vertical lines are indicated as in Fig. 5.
In the text

	Fig. 14. STECKMAP best-fit spectrum for CONT models at 4.4, 180, and 320 Myr and S/N = 10 (top panels), and 50 (bottom panels). The black and red lines show the input and the best-fit spectrum, respectively. The blue, magenta, and green lines show the mask, the residual and the errors, respectively.
In the text

	Fig. 15. Evolution of the BJS (see Eq. (5)) as a function of galaxy ages for CONT models at S/N = 5 (black), 10 (yellow), 50 (magenta), and 250 (cyan). Horizontal lines show the subgroup identifying ages with weak (0.92 < BJS < 0.98, red line), intermediate (1.08 < BJS < 1.12, green line) and strong (1.2 < BJS < 2, blue line) BJS.
In the text

	Fig. 16. Δ⟨log(t*)⟩M obtained from STECKMAP vs. galaxy age at different S/N, reported on the top right of each panel. The blue circle, green triangles, and red stars point out the synthetic REBETIKO spectrum with high, intermediate, and low BJS, respectively.
In the text