11institutetext: Instituto de Astrofísica e Ciências do Espaço, Universidade
de Lisboa - OAL, Tapada da Ajuda, PT1349-018 Lisboa, Portugal 22institutetext:
Departamento de Física, Faculdade de Ciências da Universidade de Lisboa,
Edifício C8, Campo Grande, PT1749-016 Lisboa, Portugal 33institutetext:
Department of Mathematics, University of Évora, R. Romão Ramalho 59, 7000
Évora, Portugal 44institutetext: Instituto de Astrofísica e Ciências do
Espaço, Universidade do Porto - CAUP, Rua das Estrelas, PT4150-762 Porto,
Portugal 55institutetext: CENTRA/COSTAR, Instituto Superior Técnico,
Universidade de Lisboa, Av. Rovisco Pais 1, 1049-001 Lisboa, Portugal

# Self-consistent population spectral synthesis with Fado:

II. Star formation history of galaxies in spectral synthesis methods

Ciro Pappalardo 1122 Leandro S.M. Cardoso 3344 Jean Michel Gomes 44
Polychronis Papaderos 112244 José Afonso 1122 Iris Breda 11 Andrew Humphrey 44
Tom Scott 44 Stergios Amarantidis 1122 Israel Matute 1122 Rodrigo Carvajal
1122 Silvio Lorenzoni 11 Patricio Lagos 44 Ana Paulino-Afonso 55 Henrique
Miranda 22

###### Abstract

Context. The field of galaxy evolution will make a great leap forward in the
next decade as a consequence of the huge effort by the scientific community in
multi-object spectroscopic facilities. Various future surveys will enormously
increase the number of available galaxy spectra, providing new insights into
unexplored areas of research. To maximise the impact of such incoming data,
the analysis methods must also step up, extracting reliable information from
the available spectra. It is therefore urgent to refine and test reliable
analysis tools that are able to infer the properties of a galaxy from medium-
or high-resolution spectra.

Aims. In this paper we aim to investigate the limits and the reliability of
different spectral synthesis methods in the estimation of the mean stellar age
and metallicity. These two quantities are fundamental to determine the
assembly history of a galaxy by providing key insights into its star formation
history. The main question this work aims to address is which signal-to-noise
ratios (S/N) are needed to reliably determine the mean stellar age and
metallicity from a galaxy spectrum and how this depends on the tool used to
model the spectra.

Methods. To address this question we built a set of realistic simulated
spectra containing stellar and nebular emission, reproducing the evolution of
a galaxy in two limiting cases: a constant star formation rate and an
exponentially declining star formation with a single initial burst. We
degraded the synthetic spectra built from these two star formation histories
(SFHs) to different S/N and analysed with three widely used spectral synthesis
codes, namely Fado, Steckmap, and Starlight, assuming similar fitting set-ups
and the same base of spectral templates.

Results. For S/N $\leq 5$ all the three tools show a large diversity in the
results. The Fado and Starlight tools find median differences in the light-
weighted mean stellar age of $\sim$0.1 dex, while Steckmap shows a higher
value of $\sim$0.2 dex. For S/N $>50$ the median differences in Fado are
$\sim$0.03 dex ($\sim 7$%), a factor 3 and 4 lower than the 0.08 dex
($\sim$20%) and 0.11 dex ($\sim$30%) obtained from Starlight and Steckmap,
respectively. Detailed investigations of the best-fit spectrum for galaxies
with overestimated mass-weighted quantities point towards the inability of
purely stellar models to fit the observed spectral energy distribution around
the Balmer jump.

Conclusions. Our results imply that when a galaxy enters a phase of high
specific star formation rate (sSFR) the neglect of the nebular continuum
emission in the fitting process has a strong impact on the estimation of its
SFH when purely stellar fitting codes are used, even in presence of high S/N
spectra. The median value of these differences are of the order of 7% (Fado),
20% (Starlight), and 30% (Steckmap) for light-weighted quantities, and 20%
(Fado), 60% (Starlight), and 20% (Steckmap) for mass-weighted quantities. More
specifically, for a continuous SFH both Steckmap and Starlight overestimate
the stellar age by $>$2 dex within the first $\sim$100 Myr even for high S/N
spectra. This bias, which stems from the neglect of nebular continuum
emission, obviously implies a severe overestimation of the mass-to-light ratio
and stellar mass. 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. Our work underlines once
more the importance of a self-consistent treatment of nebular emission, as
implemented in FADO, which, according to our analysis, is the only viable
route towards a reliable determination of the assembly of any high-sSFR galaxy
at high and low redshift.

###### Key Words.:

Star formation; Galaxies Evolution; Spectral Synthesis Codes

##  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
A. 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
$\sim$ 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 Sec. 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 Section 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{split}F(\lambda,t)=\int_{0}^{t}\Phi(t-t^{\prime})&F^{SSP}_{\lambda}(t^{\prime},Z)dt^{\prime}\\\
&+\frac{\gamma_{eff}(T_{e})}{\alpha_{B}(T_{e})}\int_{0}^{t}\Phi(t-t^{\prime})q^{SSP}(t^{\prime},Z)dt^{\prime}\end{split}.$
(1)

Eq. 1 reproduces the evolution of the SED of a galaxy as a function of time.
For each $t$ the flux $F(\lambda,t)$ at a specific wavelength $\lambda$
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; $\Phi$). 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 $\gamma_{eff}$ to
the photons recombination rate, $\alpha_{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-3 and a temperature of $T_{e}=10^{4}$ 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-1. 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$\alpha$ and H$\beta$ (within $\pm$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 $\sim 10^{8}$ Myr, followed by an exponentially declining SFR
($\Phi$ = e-t/τ, with $\tau$ = 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_{\odot}$ = 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 $\Delta\lambda$ = 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^{7}$
yr) and are impossible to fit with purely stellar templates (Cardoso, Gomes &
Papaderos 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 AV = 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, Gomes & Papaderos
2019).

Figure 1: Synthetic spectra normalised at $\lambda$ = 4020 Å for TAU1
(instantaneous burst, top panel) and CONT (continuous, bottom panel) SFH
models. Figures adapted from Fig. 1 and Fig. 2 of Cardoso, Gomes & Papaderos
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{split}\left<log(t_{\ast})\right>_{L}=\sum_{j=1}^{N_{\ast}}x_{j}\cdot
log(t_{j}),\\\
\left<log(t_{\ast})\right>_{M}=\sum_{j=1}^{N_{\ast}}\mu_{j}\cdot
log(t_{j})\end{split}.$ (2)

Analogously for the metallicity,

$\begin{split}\left<log(Z_{\ast})\right>_{L}=\sum_{j=1}^{N_{\ast}}x_{j}\cdot
log(Z_{j}),\\\
\left<log(Z_{\ast})\right>_{M}=\sum_{j=1}^{N_{\ast}}\mu_{j}\cdot
log(Z_{j}).\end{split}$ (3)

In Eq. 2 and 3, $t_{j}$ ($Z_{j}$) represents the age (metallicity) of the j-th
SSP element, where $x_{j}$ and $\mu_{j}$ are their associated light and mass
fractions at 4020 Å, respectively. Specifically, $\mu_{j}$ denotes the
existent stellar mass $M_{\ast}$ 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{split}\Delta_{L}=\left<log(t_{\ast})\right>_{L}-log(t)\\\
\Delta_{M}=\left<log(t_{\ast})\right>_{M}-log(t),\end{split}$ (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. Fig. 2 shows
the evolution of such parameters for the CONT model.

Figure 2: Differences between luminosity-weighted (solid black line) and mass-
weighted (dashed red line) mean stellar ages and the galaxy ages ($\Delta_{L}$
and $\Delta_{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 $\left<\Delta\right>_{L}$ (black text) and
$\left<\Delta\right>_{M}$ (red text) for young ($7<\log(t/\mathrm{yr})$),
intermediate($7<\log(t/\mathrm{yr})<9$), and old ($\log(t/\mathrm{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 108 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
107 yr when star formation begins and the mass and luminosity of the stellar
populations are growing. After 107 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/\mathrm{yr})$, intermediate ($7<\log(t/\mathrm{yr})<9$), and old
($\log(t/\mathrm{yr})>9$) ages (see also Westera et al. 2004).

Fig. 2 shows that after 107 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 $\left<log(t_{\ast})\right>_{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
$\left<log(t_{\ast})\right>_{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$\alpha$ line emission
line EW, EW(H$\alpha$), 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$\alpha$ line. Fig. 3 reports the
evolution of EW(H$\alpha$) 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.

The Lyman continuum production rate, hence the EW(H$\alpha$) 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$\alpha$)
$>1000\ \AA$, 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$\alpha$) $<60\ \AA$, 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.

Figure 3: Evolution of EW(H$\alpha$) 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.

##  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 Sec. 3.3). In the following section we
describe the most important aspects of each code, which are summarised in
Table 1.

| Fado | Starlight | Steckmap
---|---|---|---
Algorithm used | DEO111Differential evolution optimisation (Storn & Price 1996, 1997) | SA+MS222Simulated annealing plus Metropolis scheme (MacKay 2003) | MAP333Maximum a posteriori method
Nebular emission | Yes | No | No
Automatic masking | Yes | No | No
Predefined ext. law | Yes | Yes | Yes
Gaussian kernel | Yes | Yes | Yes
Table 1: Main features of the spectral synthesis codes used in this paper.

###  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$\alpha$, H$\beta$ luminosity and EW). This approach reduces
the biases in the determination of mass and light-weighted stellar ages, as
shown in Cardoso, Gomes & Papaderos (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, Gomes &
Papaderos 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.

Figure 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.

##  4 Results

The synthetic spectra built with Rebetiko in Section 2 were analysed with the
tools described in Section 3, for six S/N values in the range $3\leq$ S/N
$\leq$ 100\. For each SFH we estimated the mean stellar age and metallicity,
weighted by mass or light, as defined in Eq. 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.

tool | S/N | SFH | $\Delta t_{L}$ | $\Delta t_{M}$ | $\Delta Z_{L}$ | $\Delta Z_{M}$ | SFH | $\Delta t_{L}$ | $\Delta t_{M}$ | $\Delta Z_{L}$ | $\Delta Z_{M}$
---|---|---|---|---|---|---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12)
Fado | 3 | CONT | 0.09 | 0.16 | 0.112 | 0.093 | TAU1 | 0.03 | 0.05 | 0.041 | 0.032
| 5 |  | 0.10 | 0.13 | 0.103 | 0.058 |  | 0.02 | 0.03 | 0.031 | 0.020
| 10 |  | 0.08 | 0.12 | 0.084 | 0.047 |  | 0.02 | 0.02 | 0.021 | 0.013
| 20 |  | 0.06 | 0.11 | 0.072 | 0.035 |  | 0.01 | 0.02 | 0.017 | 0.007
| 50 |  | 0.04 | 0.09 | 0.049 | 0.020 |  | 0.01 | 0.02 | 0.008 | 0.005
| 100 |  | 0.03 | 0.07 | 0.037 | 0.013 |  | 0.01 | 0.01 | 0.007 | 0.004
Starlight | 3 |  | 0.11 | 0.15 | 0.111 | 0.096 |  | 0.03 | 0.06 | 0.036 | 0.028
| 5 |  | 0.10 | 0.12 | 0.110 | 0.079 |  | 0.02 | 0.04 | 0.022 | 0.018
| 10 |  | 0.08 | 0.08 | 0.093 | 0.059 |  | 0.01 | 0.03 | 0.014 | 0.010
| 20 |  | 0.08 | 0.14 | 0.082 | 0.043 |  | 0.01 | 0.03 | 0.009 | 0.008
| 50 |  | 0.08 | 0.20 | 0.068 | 0.044 |  | 0.01 | 0.03 | 0.007 | 0.006
| 100 |  | 0.08 | 0.24 | 0.070 | 0.042 |  | 0.01 | 0.03 | 0.007 | 0.005
Steckmap | 3 |  | 0.21 | 0.37 | 0.081 | 0.081 |  | 0.09 | 0.26 | 0.026 | 0.026
| 5 |  | 0.19 | 0.33 | 0.055 | 0.055 |  | 0.08 | 0.20 | 0.021 | 0.021
| 10 |  | 0.16 | 0.23 | 0.032 | 0.032 |  | 0.05 | 0.14 | 0.007 | 0.007
| 20 |  | 0.15 | 0.19 | 0.026 | 0.026 |  | 0.04 | 0.11 | 0.003 | 0.003
| 50 |  | 0.13 | 0.10 | 0.025 | 0.025 |  | 0.02 | 0.06 | $<10^{-4}$ | $<10^{-4}$
| 100 |  | 0.11 | 0.07 | 0.023 | 0.024 |  | 0.01 | 0.04 | $<10^{-4}$ | $<10^{-4}$
Table 2: Median of the logarithmic difference between the values of
$\left<log(t_{\ast})\right>_{L}$, $\left<log(t_{\ast})\right>_{M}$,
$\left<log(Z_{\ast})\right>_{L}$, and $\left<log(Z_{\ast})\right>_{M}$ in dex
recovered by each tool and the Rebetiko input parameter. Col. 1: Analysis
tool; col. 2: Signal-to-noise ratio of the input Rebetiko spectra; cols. 3 and
8: SFH of the models; cols. 4-7: median
$\Delta$$\left<log(t_{\ast})\right>_{L}$,
$\Delta$$\left<log(t_{\ast})\right>_{M}$,
$\Delta$$\left<log(Z_{\ast})\right>_{L}$, and
$\Delta$$\left<log(Z_{\ast})\right>_{M}$ obtained from CONT models; cols.
9-12: median $\Delta$$\left<log(t_{\ast})\right>_{L}$,
$\Delta$$\left<log(t_{\ast})\right>_{M}$,
$\Delta$$\left<log(Z_{\ast})\right>_{L}$, and
$\Delta$$\left<log(Z_{\ast})\right>_{M}$ obtained from TAU1 models.

###  4.1 Continuous star formation: CONT models

Fig. 4 shows the difference between $\left<log(t_{\ast})\right>_{L}$,
$\left<log(t_{\ast})\right>_{M}$, $\left<log(Z_{\ast})\right>_{L}$, and
$\left<log(Z_{\ast})\right>_{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
$\Delta$$\left<log(t_{\ast})\right>_{L}$$\sim$0.1 dex, while STECKMAP have a
higher $\Delta$$\left<log(t_{\ast})\right>_{L}$of 0.2. All tools have some
problems at ages $\log(t/\mathrm{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 $\chi^{2}$ maps, where
spectral fitting algorithms could easily remain stuck in a local minimum. From
Table 2 we see that Fado and Starlight underestimate
$\left<log(t_{\ast})\right>_{L}$ by $\sim$ 0.1 dex, and in Fado this value is
driven mainly by the large errors at $\log(t/\mathrm{yr})<8$, where results
vary by $\sim$ 0.3-0.4 dex ($\sim$50%). For older ages, the results become
more stable and have a mild tendency to underestimate
$\left<log(t_{\ast})\right>_{L}$. For S/N = 3 Steckmap shows the highest
discrepancies, have a median $\Delta$$\left<log(t_{\ast})\right>_{L}$$\sim$
0.2 dex, and large differences in the results.

At S/N = 5 the solutions span a similar range in terms of
$\Delta$$\left<log(t_{\ast})\right>_{L}$, varying between 0.1 and 0.2 dex. The
uncertainties of the fits at galaxy ages $\log(t/\mathrm{yr})<7.5$ are still
present; Fado and Starlight show some spread. However, the median
$\Delta$$\left<log(t_{\ast})\right>_{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/\mathrm{yr})<8$ disappear, decreasing the
$\Delta$$\left<log(t_{\ast})\right>_{L}$$\sim$0.04 dex and keeping results
differences under the 10%.

For S/N $>20$ the median differences between the input Rebetiko
$\left<log(t_{\ast})\right>_{L}$ and the best fits from Fado are of the order
of 0.03 dex ($\sim 7$%), a factor 3 and 4 lower than the 0.08 dex ($\sim$20%),
and 0.12 dex ($\sim$30%) obtained from Starlight and Steckmap, respectively.
In contrast to Fado, Starlight shows a systematic underestimation of the input
$\left<log(t_{\ast})\right>_{L}$ by 0.08 dex, 20%, for S/N $>$ 10\. The fits
show small variations, and the $\Delta$$\left<log(t_{\ast})\right>_{L}$ at
ages above $\log(t/\mathrm{yr})\sim 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/\mathrm{yr})>9$. An interesting aspect
to note are the systematic uncertainties in determining the mean stellar ages
at $\log(t/\mathrm{yr})\sim 7$, that is epochs with relevant nebular
contribution to the overall continuum (see Fig. 3).

Figure 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 $\left<log(t_{\ast})\right>_{L}$ (black), and the values recovered
with Steckmap (green), Fado (blue), and Starlight (red). The top left of each
panel indicate the $\chi^{2}$ of the Steckmap best-fit model.

The Steckmap tool shows median differences of the order of 0.15 dex
($\sim$40%) at all S/N, as shown in Table 2. On average the results are very
stable between $7<\log(t/\mathrm{yr})<9$ yr, with underestimation of
$\left<log(t_{\ast})\right>_{L}$$\sim$0.08 dex at $\log(t/\mathrm{yr})\leq 8$,
similarly to that seen for Starlight, and an overestimation of $\sim$20% for
ages $\log(t/\mathrm{yr})>9$, owing to an effect already mentioned in Cardoso,
Gomes & Papaderos (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/\mathrm{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
$\left<log(t_{\ast})\right>_{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 $\left<log(t_{\ast})\right>_{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/\mathrm{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/\mathrm{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 $\left<log(t_{\ast})\right>_{L}$ for those
galaxy ages (see Fig. 5).

The third row panels of Fig. 4 show the difference of luminosity-weighted mean
stellar metallicity $\left<log(Z_{\ast})\right>_{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/\mathrm{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, $\Delta$$\left<log(Z_{\ast})\right>_{L}$$\sim 0.2$ dex, for
ages below $\log(t/\mathrm{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/\mathrm{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 $\sim$0.09 dex, 25%, at S/N $\geq 10$, which is slightly larger than
the 0.04 dex seen for $\left<log(t_{\ast})\right>_{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/\mathrm{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 $\Delta$$\left<log(t_{\ast})\right>_{M}$ for
$\log(t/\mathrm{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
Section 5. We confirm an overestimation of the
$\left<log(t_{\ast})\right>_{M}$ for ages above 3 Gyr, as seen already in
Cardoso, Gomes & Papaderos (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
$\Delta$$\left<log(t_{\ast})\right>_{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/\mathrm{yr})>9$, a feature seen also in Cardoso, Gomes &
Papaderos (2019). The Starlight best fits systematically underestimate the
metallicity by $\sim$0.01 dex, in particular at young ages, and this trend
rises to 0.6 dex for $\log(t/\mathrm{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: Same caption of Fig. 4 for TAU1 SFH 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 Sec. 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 $\Delta$$\left<log(t_{\ast})\right>_{L}$ is approximately 2-3% for Fado
and Starlight already at S/N = 5. Slight oscillations are seen in Steckmap at
ages $\sim$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, Gomes &
Papaderos 2019). The results of Steckmap are dependent on the input S/N; a
reasonable accuracy is reached at S/N $>$ 10\. For lower S/N ratios the fits
are more unstable and we observe an overestimation of
$\left<log(t_{\ast})\right>_{L}$ at $\log(t/\mathrm{yr})>7.5$, whose origins
are clarified in Section 5.

The results for mean metallicities are more homogeneous than CONT models
between the different tools. The solutions are more unstable at
$\log(t/\mathrm{yr})\sim 7$; Fado and Starlight show a peak of $\sim$0.2 dex
at $\sim$200 Myr in $\left<log(Z_{\ast})\right>_{L}$ and
$\left<log(Z_{\ast})\right>_{M}$ resulting from the poor coverage of the
spectral basis at those ages (see Cardoso, Gomes & Papaderos 2019). Curiously,
Steckmap has similar problems when fitting mean metallicity at ages around
$\log(t/\mathrm{yr})\sim 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/\mathrm{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
$\Delta$$\left<log(t_{\ast})\right>_{L}$ and
$\Delta$$\left<log(t_{\ast})\right>_{M}$, $\sim 0.1-0.2$ dex, while for
Starlight and Steckmap the average uncertainties increase by a factor of
$\sim$3\. Both methods show large uncertainties at ages below
$\log(t/\mathrm{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 $\geq$ 5.

##  5 Discussion

As we have seen in Sec. 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 $\left<log(t_{\ast})\right>_{L}$ and
$\left<log(t_{\ast})\right>_{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
$\left<log(t_{\ast})\right>_{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, Gomes & Papaderos 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 Section 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
Section 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
$\left<log(t_{\ast})\right>_{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 Sec. 3 (yellow line), and MILES (Vazdekis et al. 2010, green
line).

Figure 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
108 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 $\sim 10^{8}$ 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.

TAU1 Models  tool S/N SSP $\Delta t_{L}$ $\Delta t_{M}$ SSP $\Delta t_{L}$
$\Delta t_{M}$ (1) (2) (3) (4) (5) (6) (7) (8) Steckmap 5 BL 0.08 0.2 ML 0.09
0.23 20 0.04 0.11 0.04 0.19 100 0.01 0.04 0.12 0.24

Table 3: Median of the logarithmic difference between the values of
$\left<log(t_{\ast})\right>_{L}$and $\left<log(t_{\ast})\right>_{M}$ in dex
recovered by Steckmap and the Rebetiko input parameter for TAU1 models. Col.
1: Analysis tool; col. 2: Signal-to-noise ratio of the input Rebetiko spectra;
cols. 3, and 6: spectral basis used: BaseL as defined in Sec. 3 (BL), and
MILES (ML, Vazdekis et al. 2010); cols. 4-5-7-8: median
$\Delta$$\left<log(t_{\ast})\right>_{L}$,
$\Delta$$\left<log(t_{\ast})\right>_{M}$ for each SSP.

The problem of the spectral basis considered in a fit has been thoroughly
addressed in Cardoso, Gomes & Papaderos (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.

Figure 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. Figure 9: Same caption of Fig. 8 for TAU1 SFH models.

Concerning the wavelength coverage issue, Figures 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, $\sim$4800-9300 Å, including
H$\alpha$+H$\beta$ 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.

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
$\left<log(t_{\ast})\right>_{M}$ for $\log(t/\mathrm{yr})<8$ (see Fig. 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: 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 $\pm 1\sigma$ and $\pm 3\sigma$
error spectrum, respectively.

Fig. 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 $\sim$3646 Å. The contribution of the
nebular continuum increases with increasing wavelength, becoming dominant in
the region around the Paschen jump at $\sim$ 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.

SFH = CONT

SFH = TAU1

Figure 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.

For luminosity-weighted quantities this effect is less relevant, even in the
presence of a strong nebular component. Fig. 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/\mathrm{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 $\left<log(t_{\ast})\right>_{L}$, which is recovered with an
error of $\sim$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.

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 $\left<log(t_{\ast})\right>_{M}$, in agreement with
the results of (Cardoso, Gomes & Papaderos 2019, see their Fig. B.1). These
differences, as mentioned at the beginning of Section 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
$\left<log(t_{\ast})\right>_{L}$ and $\left<log(t_{\ast})\right>_{M}$, as
shown in Fig. 12. In other words where the nebular contribution is relevant
there is a systematic overestimation of $\left<log(t_{\ast})\right>_{M}$,
which has led to an overestimation of the stellar masses (see also Gomes &
Papaderos 2017 and Cardoso, Gomes & Papaderos 2019).

SFH = CONT

SFH = TAU1

Figure 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
$\Delta$$\left<log(t_{\ast})\right>_{M}$ measured with Steckmap span a
relatively extended range of differences from Rebetiko input parameters, for
both CONT and TAU1 models: from $\sim$ 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 $\Delta$$\left<log(t_{\ast})\right>_{M}$ decreases from
0.23 ($\sim$70%) up to 0.1 ($\sim$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/\mathrm{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/\mathrm{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.

SFH = CONT

SFH = TAU1

Figure 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.

To investigate the origin of such S/N-linked overestimation of
$\left<log(t_{\ast})\right>_{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. Fig. 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.

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 $\left<log(t_{\ast})\right>_{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 $\left<log(t_{\ast})\right>_{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.

S/N = 10

S/N = 50

Figure 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.

Increasing the S/N, Steckmap improves the quality of the fits, reproducing
more the input spectrum more accurately. Fig. 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 Å.

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

Fig. 14 shows that one of the causes of the high
$\left<log(t_{\ast})\right>_{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

$\mathrm{BJS}=\frac{\int_{3580\AA}^{3600\AA}F(\lambda,t)d\lambda}{\int_{3660\AA}^{3680\AA}F(\lambda,t)d\lambda}.$
(5)

Figure 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.

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 $\left<BJS\right>$ = 1.06$\pm$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 .

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
$\Delta$$\left<log(t_{\ast})\right>_{M}$ analysed in Sec. 4 with the BJS.

Figure 16: $\Delta$$\left<log(t_{\ast})\right>_{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.

Fig. 4 and Fig. 16 show the residual of $\left<log(t_{\ast})\right>_{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 $\left<log(t_{\ast})\right>_{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$\alpha$ 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 $\Delta$$\left<log(t_{\ast})\right>_{M}$ at ages higher than
$\log(t/\mathrm{yr})\sim 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 $\left<log(t_{\ast})\right>_{M}$ with
an accuracy better than 0.4 dex at S/N $\geq$ 5, as shown with red stars in
Fig. 16.

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 $\Delta$$\left<log(t_{\ast})\right>_{L}$$\sim$0.1 dex and
STECKMAP has a higher $\Delta$$\left<log(t_{\ast})\right>_{L}$of 0.2. The
codes have a large spread of $\left<log(t_{\ast})\right>_{L}$ at ages
$\log(t/\mathrm{yr})<8$, as a consequence of the instability of the solutions.

  * •

For S/N $>20$ the median differences between the input Rebetiko
$\left<log(t_{\ast})\right>_{L}$ and the best fits from Fado are of the order
of 0.03 dex ($\sim 7$%), which is a factor 3 and 4 lower than the 0.08 dex
($\sim$20%) and 0.11 dex ($\sim$30%) obtained from Starlight and Steckmap,
respectively. In contrast to Fado, Starlight shows a systematic
underestimation of the input $\left<log(t_{\ast})\right>_{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 $\geq
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
$\Delta$$\left<log(t_{\ast})\right>_{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 Sec. 3 with MILES (Vazdekis et al. 2010), we
see that the main differences in both light- and mass-weighted quantities
arise for ages below 108 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 $\left<log(t_{\ast})\right>_{M}$,
which has led to an overestimation of the stellar masses up to 50% for $t>300$
Myr (see also Gomes & Papaderos 2017 and Cardoso, Gomes & Papaderos 2019).

Detailed investigations of the best-fit spectrum in galaxies with
overestimated $\left<log(t_{\ast})\right>_{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
$\left<log(t_{\ast})\right>_{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, Gomes &
Papaderos 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.

###### Acknowledgements.

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

  * Ahumada et al. (2020) Ahumada, R., Allende Prieto, C., Almeida, A., et al. 2020, ApJS, 249, 3
  * Alsing et al. (2020) Alsing, J., Peiris, H., Leja, J., et al. 2020, ApJS, 249, 5. doi:10.3847/1538-4365/ab917f
  * Alongi et al. (1993) Alongi, M., Bertelli, G., Bressan, A., et al. 1993, A&AS, 97, 851
  * Bacon et al. (2010) Bacon, R., Accardo, M., Adjali, L., et al. 2010, Proc. SPIE, 7735, 773508. doi:10.1117/12.856027
  * Baldwin et al. (2018) Baldwin, C., McDermid, R. M., Kuntschner, H., et al. 2018, MNRAS, 473, 4698. doi:10.1093/mnras/stx2502
  * Binette et al. (1994) Binette, L., Magris, C. G., Stasińska, G., et al. 1994, A&A, 292, 13
  * Breda & Papaderos (2018) Breda, I., & Papaderos, P. 2018, A&A, 614, A48
  * Bressan et al. (1993) Bressan, A., Fagotto, F., Bertelli, G., et al. 1993, A&AS, 100, 647
  * Brinchmann et al. (2004) Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151. doi:10.1111/j.1365-2966.2004.07881.x
  * Bruzual & Charlot (2003) Bruzual, G., & Charlot, S. 2003, MNRAS, 344, 1000
  * Bruzual A. (2010) Bruzual A., G. 2010, Philosophical Transactions of the Royal Society of London Series A, 368, 783
  * Cardoso, Gomes & Papaderos (2019) Cardoso, L. S. M., Gomes, J. M., & Papaderos, P. 2019, A&A, 622, A56
  * Cardoso et al. (2017) Cardoso, L. S. M., Gomes, J. M., & Papaderos, P. 2017, A&A, 604, A99
  * Chabrier (2003) Chabrier, G. 2003, ApJ, 586, L133
  * Chen et al. (2019) Chen, Y.-Y., Chen, Y.-M., Gu, Q.-S., et al. 2019, Research in Astronomy and Astrophysics, 19, 081. doi:10.1088/1674-4527/19/6/81
  * Cicone et al. (2016) Cicone, C., Maiolino, R., & Marconi, A. 2016, A&A, 588, A41. doi:10.1051/0004-6361/201424514
  * Cid Fernandes et al. (2005) Cid Fernandes, R., Mateus, A., Sodré, L., et al. 2005, MNRAS, 358, 363
  * Cid Fernandes et al. (2011) Cid Fernandes, R., Stasińska, G., Mateus, A., et al. 2011, MNRAS, 413, 1687
  * Cirasuolo & MOONS Consortium (2016) Cirasuolo, M., & MOONS Consortium 2016, Multi-object Spectroscopy in the Next Decade: Big Questions, Large Surveys, and Wide Fields, 109
  * Cirasuolo et al. (2020) Cirasuolo, M., Fairley, A., Rees, P., et al. 2020, The Messenger, 180, 10
  * Coelho et al. (2020) Coelho, P. R. T., Bruzual, G., & Charlot, S. 2020, MNRAS, 491, 2025
  * Conroy (2013) Conroy, C. 2013, ARA&A, 51, 393
  * Dalton (2016) Dalton, G. 2016, Multi-object Spectroscopy in the Next Decade: Big Questions, Large Surveys, and Wide Fields, 97
  * Doi et al. (2010) Doi, M., Tanaka, M., Fukugita, M., et al. 2010, AJ, 139, 1628
  * de Jong et al. (2019) de Jong, R. S., Agertz, O., Berbel, A. A., et al. 2019, The Messenger, 175, 3
  * Duckworth et al. (2020) Duckworth, C., Tojeiro, R., & Kraljic, K. 2020, MNRAS, 492, 1869. doi:10.1093/mnras/stz3575
  * Faber (1972) Faber, S. M. 1972, A&A, 20, 361
  * Fagotto et al. (1994a) Fagotto, F., Bressan, A., Bertelli, G., et al. 1994, A&AS, 104, 365
  * Fagotto et al. (1994b) Fagotto, F., Bressan, A., Bertelli, G., et al. 1994, A&AS, 105, 29
  * Florian et al. (2020) Florian, J., Ziegler, B., Hirschmann, M., et al. 2020, A&A, 635, A41. doi:10.1051/0004-6361/201936441
  * Girardi et al. (1996) Girardi, L., Bressan, A., Chiosi, C., et al. 1996, A&AS, 117, 113
  * Gomes & Papaderos (2016) Gomes, J. M. & Papaderos, P. 2016, A&A, 594, A49. doi:10.1051/0004-6361/201628316
  * Gomes et al. (2016) Gomes, J. M., Papaderos, P., Kehrig, C., et al. 2016, A&A, 588, A68
  * Gomes & Papaderos (2017) Gomes, J. M., & Papaderos, P. 2017, A&A, 603, A63
  * Guseva et al. (2007) Guseva, N. G., Izotov, Y. I., Papaderos, P., et al. 2007, A&A, 464, 885
  * Heavens et al. (2004) Heavens, A., Panter, B., Jimenez, R., et al. 2004, Nature, 428, 625
  * Iyer et al. (2019) Iyer, K. G., Gawiser, E., Faber, S. M., et al. 2019, ApJ, 879, 116
  * Izotov et al. (1997) Izotov, Y. I., Lipovetsky, V. A., Chaffee, F. H., et al. 1997, ApJ, 476, 698
  * Izotov et al. (2011) Izotov, Y. I., Guseva, N. G., & Thuan, T. X. 2011, ApJ, 728, 161
  * Jin et al. (2016) Jin, Y., Chen, Y., Shi, Y., et al. 2016, MNRAS, 463, 913. doi:10.1093/mnras/stw2055
  * Kelz et al. (2006) Kelz, A., Verheijen, M. A. W., Roth, M. M., et al. 2006, PASP, 118, 129. doi:10.1086/497455
  * Koleva et al. (2008) Koleva, M., Prugniel, P., Ocvirk, P., et al. 2008, MNRAS, 385, 1998
  * Koleva et al. (2009) Koleva, M., de Rijcke, S., Prugniel, P., et al. 2009, MNRAS, 396, 2133
  * Krueger et al. (1995) Krueger, H., Fritze-v. Alvensleben, U., & Loose, H.-H. 1995, A&A, 303, 41
  * Le Borgne et al. (2003) Le Borgne, J.-F., Bruzual, G., Pelló, R., et al. 2003, A&A, 402, 433. doi:10.1051/0004-6361:20030243
  * Leitherer et al. (1999) Leitherer, C., Schaerer, D., Goldader, J. D., et al. 1999, ApJS, 123, 3
  * Leja et al. (2017) Leja, J., Johnson, B. D., Conroy, C., et al. 2017, ApJ, 837, 170
  * López Fernández et al. (2016) López Fernández, R., Cid Fernandes, R., González Delgado, R. M., et al. 2016, MNRAS, 458, 184
  * Maiolino et al. (2020) Maiolino, R., Cirasuolo, M., Afonso, J., et al. 2020, The Messenger, 180, 24
  * MacKay (2003) MacKay D. J. C., 2003, Information Theory, Inference and Learning Algorithms. Cambridge Univ. Press, Cambridge
  * Martinsson et al. (2013) Martinsson, T. P. K., Verheijen, M. A. W., Westfall, K. B., et al. 2013, A&A, 557, A130. doi:10.1051/0004-6361/201220515
  * McGurk et al. (2010) McGurk, R. C., Kimball, A. E., & Ivezić, Ž. 2010, AJ, 139, 1261
  * Ocvirk et al. (2006a) Ocvirk, P., Pichon, C., Lançon, A., et al. 2006, MNRAS, 365, 46
  * Ocvirk et al. (2006b) Ocvirk, P., Pichon, C., Lançon, A., et al. 2006, MNRAS, 365, 74
  * Ocvirk (2010) Ocvirk, P. 2010, ApJ, 709, 88
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of Gaseous Nebulae and Active Galactic Nuclei
  * Osterbrock & Ferland (2006) Osterbrock, D. E., & Ferland, G. J. 2006, Astrophysics of gaseous nebulae and active galactic nuclei
  * Pacifici et al. (2015) Pacifici, C., da Cunha, E., Charlot, S., et al. 2015, MNRAS, 447, 786
  * Panter et al. (2003) Panter, B., Heavens, A. F., & Jimenez, R. 2003, MNRAS, 343, 1145
  * Panter et al. (2007) Panter, B., Jimenez, R., Heavens, A. F., et al. 2007, MNRAS, 378, 1550
  * Papaderos et al. (1998) Papaderos, P., Izotov, Y. I., Fricke, K. J., et al. 1998, A&A, 338, 43
  * Papaderos & Östlin (2012) Papaderos, P., & Östlin, G. 2012, A&A, 537, A126
  * Pappalardo et al. (2010) Pappalardo, C., Lançon, A., Vollmer, B., et al. 2010, A&A, 514, A33
  * Sánchez et al. (2012) Sánchez, S. F., Kennicutt, R. C., Gil de Paz, A., et al. 2012, A&A, 538, A8. doi:10.1051/0004-6361/201117353
  * Sarzi et al. (2007) Sarzi, M., Bacon, R., Cappellari, M., et al. 2007, New A Rev., 51, 18. doi:10.1016/j.newar.2006.11.023
  * Schaerer & de Barros (2009) Schaerer, D., & de Barros, S. 2009, A&A, 502, 423
  * Stasińska et al. (2008) Stasińska, G., Vale Asari, N., Cid Fernandes, R., et al. 2008, MNRAS, 391, L29
  * Storn & Price (1996) Storn, R., & Price, K. V. , K. V. 1996, Minimizing the real function of the ICEC’96 contest by differential evolution, IEEE Conf. on Evolutionary Computation, 842, 844
  * Storn & Price (1997) Storn, R., & Price, K. V. 1997, J. Global Optimization, 11, 34
  * Struck-Marcell & Tinsley (1978) Struck-Marcell, C., & Tinsley, B. M. 1978, ApJ, 221, 562
  * Takada et al. (2014) Takada, M., Ellis, R. S., Chiba, M., et al. 2014, PASJ, 66, R1
  * Tinsley (1980) Tinsley, B. M. 1980, Fund. Cosmic Phys., 5, 287
  * Tojeiro et al. (2007) Tojeiro, R., Heavens, A. F., Jimenez, R., et al. 2007, MNRAS, 381, 1252
  * Tumlinson et al. (2017) Tumlinson, J., Peeples, M. S., & Werk, J. K. 2017, ARA&A, 55, 389
  * Vazdekis et al. (2010) Vazdekis, A., Sánchez-Blázquez, P., Falcón-Barroso, J., et al. 2010, MNRAS, 404, 1639. doi:10.1111/j.1365-2966.2010.16407.x
  * Yuan et al. (2019) Yuan, F.-T., Burgarella, D., Corre, D., et al. 2019, A&A, 631, A123
  * Westera et al. (2004) Westera, P., Cuisinier, F., Telles, E., et al. 2004, A&A, 423, 133
  * Wilkinson et al. (2017) Wilkinson, D. M., Maraston, C., Goddard, D., et al. 2017, MNRAS, 472, 4297