Revisiting stellar properties of star-forming galaxies with stellar and nebular spectral modelling

Spectral synthesis is a powerful tool for interpreting the physical properties of galaxies by decomposing their spectral energy distributions into the main luminosity contributors (e.g. stellar populations or ionised gas). However, the impact nebular emission has on the inferred properties of star-forming (SF) galaxies has been largely overlooked over the years. The objective of this work is to estimate the relations between stellar properties of SF galaxies from SDSS DR7 by simultaneously fitting the stellar and nebular continua with FADO and comparing them to the results derived using STARLIGHT, a representative of purely stellar population synthesis codes. Differences between codes regarding average mass, mean age and mean metallicity values can go as high as $\sim$0.06 dex for the overall population of galaxies and $\sim$0.12 dex for SF galaxies (galaxies with EW(H$\alpha$)>3 \AA), with the most prominent difference between both codes in the light-weighted mean stellar age. A closer look into the average light- and mass-weighted star formation histories of intensively SF galaxies (EW(H$\alpha$)>75 \AA) suggests that STARLIGHT is underestimating the average light-weighted age of intensively SF galaxies by up to $\sim$0.17 dex and overestimating the light-weighted metallicity by up to $\sim$0.13 dex compared to FADO (or vice versa). The comparison between the average stellar properties of passive, SF and intensively SF galaxy samples also reveals that differences between codes increase with increasing EW(H$\alpha$) and decreasing total stellar mass. This work finds indirect evidence that a purely stellar population synthesis approach negatively impacts the inferred stellar properties of galaxies with relatively high star formation rates. In turn, this can bias interpretations of fundamental relations such as the mass-age or mass-metallicity.


Introduction
Understanding the complex physical processes behind Galaxy formation and evolution can be a daunting endeavour, especially with the increasing technical difficulties when looking further back in time. However, several keystone results have been obtained over the past decades: (a) high-mass galaxies assemble most of their stellar content early in their lives, whereas low-mass galaxies display relatively high specific star-formation rates (sSFRs) at a late cosmic epoch (e.g. Brinchmann et al. 2004;Noeske et al. 2007a), a phenomenon known as galaxy downsizing (e.g. Cowie et al. 1996;Heavens et al. 2004; that is usually referred to as the SF main sequence of galaxies (e.g. Brinchmann et al. 2004;Noeske et al. 2007a).
Notwithstanding these insights, the overall picture of galaxy formation and evolution is far from complete. For instance, a persistent issue in interpreting the physical characteristics of SF galaxies from spectral synthesis has been the relative contribution of nebular emission and its potential impact on the estimated stellar and overall galaxy properties. Indeed, several studies have shown that nebular emission (continuum and lines) can account up to ∼60% of the monochromatic luminosity at ∼5000 Å in AGN and starburst galaxies (e.g. Krüger, Fritze-v. Alvensleben & Loose 1995;Papaderos et al. 1998;Zackrisson et al. 2001;Zackrisson, Bergvall & Leitet 2008;Schaerer & de Barros 2009;Papaderos & Östlin 2012) and in optical broadband photometry, specially at low metallicities (Anders & Fritze-v. Alvensleben 2003).
For instance, Krüger, Fritze-v. Alvensleben & Loose (1995) noted that strong star-formation activity can lead nebular continuum to contribute ∼30-70% to the total optical and nearinfrared emission, whereas emission lines alone account for up to ∼45% in optical broadband photometry. From a different perspective, Reines et al. (2010) carried out spectral modelling of two young massive star clusters in NGC 4449 and showed that both the nebular continuum and line emission have a major impact on the estimated magnitudes and colours of young clusters (≲ 5 Myr), thus also affecting their age, mass and extinction estimates. Moreover, Izotov, Guseva & Thuan et al. (2011) studied a sample of 803 SF luminous compact galaxies (i.e. 'green peas') from the SDSS (York et al. 2000) and warned that a purely stellar spectral modelling of such objects can lead to the overestimation of the relative contribution of old stellar populations by as much as a factor of four. As the authors noted, there are two reasons for this overestimation: (a) nebular continuum increases the galaxy luminosity and (b) the nebular continuum is flatter than the stellar, thus making the overall continuum redder that it would be if it were purely stellar. Papaderos & Östlin (2012) also showed that nebular emission can introduce strong photometric biases of 0.4-1 mag in galaxies with high specific SFRs, whereas Pacifici et al. (2015) found that that SFRs can be overestimated by up ∼0.12 dex if nebular emission is neglected in spectral modelling.
Although nebular continuum and line emission has been adopted in several evolutionary synthesis models (e.g. Leitherer et al. 1999;Schaerer 2002;Mollá, García-Vargas & Bressan 2009;Martín-Manjón et al. 2010), no consistent nebular prescription has been implemented in inversion population synthesis codes until recently. Indeed,  presented in Fado the first population synthesis code to fit self-consistently both the stellar and nebular spectral components while adopting a genetic optimisation framework. Tests re-vealed that this code estimates the main stellar population properties of SF galaxies (e.g. mass, mean age, and mean metallicity) within an accuracy of ∼0.2 dex Cardoso, Gomes & Papaderos 2019;Pappalardo et al. 2021). For comparison, the typical level of deviations between input and inferred stellar properties when applying purely stellar population synthesis codes to evolved stellar populations with faint or absent nebular emission is of ∼0.15-0.2 dex (e.g. Cid Fernandes et al. 2005Fernandes et al. , 2014Ocvirk et al. 2006a,b;Tojeiro et al. 2007Tojeiro et al. , 2009; Koleva et al. 2009).
Moreover, Cardoso, Gomes & Papaderos (2019) (hereafter CGP19) compared Fado with a purely stellar population synthesis code using synthetic galaxy models for different star formation histories (SFHs; e.g. instantaneous burst, continuous, and exponentially declining) and different fitting configurations (e.g. with or without emission lines and including or excluding the Balmer and Paschen continuum discontinuities at 3646 and 8207 Å, respectively). This work showed that applying the public version of the purely stellar population synthesis code Starlight 1 (Cid Fernandes et al. 2005) to spectra with a relatively high nebular continuum can lead to the overestimation of the total stellar mass by as much as ∼2 dex and the mass-weighted mean stellar age up to ∼4 dex, whereas the mean metallicity and lightweighted mean stellar age can both be underestimated by up to ∼0.6 dex. Moreover, it was found that these stellar properties can still be recovered with Fado within ∼0.2 dex in evolutionary stages with severe nebular contamination. Pappalardo et al. (2021) (hereafter P21) continued this line of inquiry by adding the non-parametric purely stellar Steckmap (Ocvirk et al. 2006a,b) to the code comparison while also exploring the impact of varying spectral quality on the derived physical properties, finding similar results and trends. It is important to note that (a) these tests were carried out using the same evolutionary models (Bruzual & Charlot 2003) both in the creation of the synthetic spectra and in the spectral modelling and (b) the input synthetic composite stellar populations were built having a constant solar metallicity (Z ⊙ = 0.02). Thus, these uncertainties should be viewed as upper limits and might not necessarily directly translate into biases affecting observations. Also recently, Gunawardhana et al. (2020) studied the stellar and nebular characteristics of massive stellar populations in the Antennae galaxy using an updated version of Platefit (Tremonti et al. 2004;Brinchmann et al. 2004), capable of self-consistent modelling of the stellar and nebular continua. Spectral fitting provides estimates of the stellar and gas metallicities, stellar ages and electron temperature T e and density n e by taking as reference model libraries of HII regions built using the evolutionary synthesis code Starburst99 (Leitherer et al. 1999) and the photoionisation code Cloudy (Ferland et al. 1998(Ferland et al. , 2013. This work found the stellar and gas metallicity of the starbursts to be near solar and that the metallicity of the star-forming gas in the loop of NGC 4038 appears to be slightly richer than the rest of the galaxy. Using a different approach, López Fernández et al. (2016) presented a new version of the population synthesis code Starlight (Cid Fernandes et al. 2005) that combines optical spectroscopy with UV photometry. This work used a mixture of simulated and real CALIFA data (Sánchez et al. 2012) and found that the additional UV constraints have a low impact on the inferred stellar mass and dust optical depth. Although the mean age and metallicity of most galaxies remains unaffected by the additional UV spectral information, this work also showed that stellar populations of low-mass late-type galaxies are older and less chemically enriched than in purely-optical modelling. Werle et al. (2019) pursued further this line of inquiry by combining GALEX photometry with SDSS spectroscopy and found that the UV constraints lead to the increase of simple stellar populations (SSPs) with ages between ∼10 7 and 10 8 years in detriment to the relative contribution of younger and older populations, leading to slightly older mean stellar ages when weighted my mass. This redistribution of the SFH is particularly noticeable in galaxies in the low-mass end of the blue cloud. Later, Werle et al. (2020) adopted a similar approach to study earlytype galaxies in the same sample and found that the UV constraints broadens the attenuation, mean stellar age and metallicity distributions. Moreover, galaxies with young stellar populations have larger Hα equivalent widths (EWs) and larger attenuations, with the metallicity of these populations being increasingly lower for larger stellar masses. Although these three works successfully use UV spectral information to constrain the contribution of young stellar populations, one wonders if the nebular continuum in galaxies with relatively high specific SFRs still needs to be accounted for as intermediate-to-old stellar populations in purely stellar spectral synthesis, in which case the lack of nebular continuum treatment would still affect the inferred stellar properties.
Taking all these factors into consideration, it is possible that the current understanding of galaxy evolution, specifically that of SF galaxies in the local universe based on large-scale survey analysis, has been affected by a lack of an adequate nebular modelling prescription in previous spectral synthesis codes. To address this subject, this work aims to revisit the relation between key stellar properties (e.g. mass, mean age, and mean metallicity) of SF galaxies by comparing the results obtained with Fado and a representative of purely stellar population synthesis codes (i.e. Starlight) when applied to a well-studied large-scale survey such as SDSS (York et al. 2000).
This paper is organised as follows. Section 2 details the methodology adopted to extract and analyse the SDSS DR7 spectra using spectral synthesis. The main results of this work are presented in Section 3, which is particularly focussed on the relation between the main stellar properties of SF galaxies when adopting two different spectral modelling approaches. Finally, in Sections 4 and 5 the findings of this work are discussed and summarised, respectively. Unless stated otherwise, this work assumes H 0 = 70 km s −1 Mpc −1 , Ω M = 0.3 and Ω Λ = 0.7.

Methodology
The population synthesis codes Fado and Starlight were applied to the galaxy sample from SDSS Data Release 7 (Abazajian et al. 2009). This data set contains multi-band spectrophotometric data for 926246 objects, each with a single-fibre integrated spectrum with the wavelength range of 3800-9200 Å at a resolution of R ∼ 1800-2200, with the fibre covering ∼5.5 kpc of the central region of an object at z ∼ 0.1. There are several reasons to utilise this dataset: (a) photometric completeness, (b) relatively wide redshift coverage (0.02 ≲ z ≲ 0.6), (c) uniform spectral calibration, and (d) wide-range of topics related to galaxy formation and evolution that have been tackled with this data (e.g. massmetallicity relation, Tremonti et al. 2004; stellar mass assembly, Kauffmann et al. 2003a;galaxy bimodality, Strateva et al. 2001;Blanton et al. 2003;Kauffmann et al. 2003b;Baldry et al. 2004;host-AGN relation, Kauffmann et al. 2003c;Hao et al. 2005; environmental effects on the colour-magnitude relation Hogg et al. 2004;Cooper et al. 2010).
The raw spectra were corrected for extinction using the colour excess correction factor E(B − V) computed based on the Schlegel, Finkbeiner & Davis (1998) dust maps, considering the objects sky position using the right ascension and declination provided by the survey, combined with the correcting factor to E(B − V) of Schlafly & Finkbeiner (2011). The full correction equation reads, where F corrected λ and F raw λ are the corrected and raw fluxes as a function of wavelength λ, respectively, and A V is the extinction in the V-band. The adopted extinction curve (also known as the reddening law) was Cardelli, Clayton & Mathis (1989), with R V = 3.1. The spectra were then converted to the observer restframe and rebinned so that ∆λ = 1 Å.
As detailed in , the main task of the population synthesis code Fado is to reproduce the observed spectral energy distribution (SED) through a linear combination of spectral components (e.g. individual stellar spectra or SSPs) as expressed by: where F λ is the flux of the observed spectrum, N ⋆ is the number of unique spectral components in the adopted base library, M i,λ 0 is the stellar mass of the i th spectral component at the normalisation wavelength λ 0 , L j,λ is the luminosity contribution of the i th spectral component, A V is the V-band extinction, q λ is the ratio of A λ over A V , S (v ⋆ , σ ⋆ ) denotes a Guassian kernel simulating the effect of stellar kinematics on the spectrum, with v ⋆ and σ ⋆ representing the stellar shift and dispersion velocities, respectively, Γ λ (n e , T e ) is the nebular continuum computed assuming that all stellar photons with λ ≤ 911.76 Å are absorbed and reprocessed into nebular emission, under the supposition that case B recombination applies, A neb V is the nebular V-band extinction, and N(v η , σ η ) denotes the nebular kinematics kernel, with v η and σ η representing the nebular shift and dispersion velocities, respectively.
It is important to note that all publicly available population synthesis codes before Fado aimed to reconstruct the observed continuum using only the first term on the right-hand side of Equation 2, meaning, using a purely stellar scheme to reconstruct the overall observed SED (more details in . Indeed, one of the innovative features of Fado is the inclusion of the second term, which represents the relative spectral contribution of the nebular continuum inferred from the Lyman continuum production rate of the young stellar component and scaling with the intensity of prominent Balmer lines. This makes Fado a tool specially designed to self-consistently model the stellar and nebular continuum components of SF and starburst galaxies.
Subsequently spectral fitting was carried out with both Fado 2 and Starlight 3 between 3400 and 8900 Å using a base library with 150 SSPs from Bruzual & Charlot (2003) with 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). This stellar library contains 25 ages (between 1 Myr and 15 Gyr) and six metallicities (Z = 0.0001; 0.0004; 0.004; 0.008; 0.02; 0.05 = Z ⊙ × {1/200; 1/50; 1/5; 2/5; 1; 2.5}) and corresponds to an expanded version of the base adopted in CGP19, with the addition of Z ⊙ /200 and Z ⊙ /50 at the lower end of the metallicity range. The spectra were modelled assuming a Calzetti et al. (2000) extinction law, which was originally constructed taking into consideration integrated observations of SF galaxies. Moreover, the V-band extinction, stellar velocity shift and dispersion free parameters where allowed to vary in both codes within the ranges of A V = 0-4 mag, v ⋆ = -500-500 km/s and σ ⋆ = 0-500 km/s, respectively. Identical parameter ranges for the velocity shift and dispersion were adopted for the nebular component in the Fado runs. Moreover, two complete runs of the whole sample were performed with Starlight while changing the initial guesses in the parameter space in order to evaluate the formal errors (e.g. Ribeiro et al. 2016;Cardoso, Gomes & Papaderos 2016. Moreover, a pre-fitting analysis was carried out with Starlight using a smaller base library with 45 SSPs for 15 ages (between 1 Myr and 13 Gyr) and three metallicities (Z = {0.004; 0.02; 0.05} = Z ⊙ × {1/5; 1; 2.5}) that comes with the distribution package of Starlight (c.f. Cid Fernandes et al. 2005). The results of this preliminary step were used to create individual spectral masks in order to exclude the most prominent emission lines in the main fitting. This was achieved by subtracting the continuum of the best-fit model from the observation and, from the resulting residual spectrum, by fitting Gaussians to the lines. The masking of these spectral regions ensures that the Starlight results are more robust when it comes to SF galaxies, a procedure similar to that adopted by Asari et al. (2007) and Ribeiro et al. (2016). At the same time, Fado performs the masking of the most prominent emission lines as part of one of its built-in prefitting routines, as detailed in . Three additional spectral regions are masked in both codes which are associated with small bugs in the adopted evolutionary models: 6845-6945, 7165-7210 and 7550-7725 Å (cf. Bruzual & Charlot 2003 for more details).
In addition to masking, spectral modelling was carried out while taking into consideration several spectral flags provided by the SDSS survey that are included in the mask array of each 'fit' file. These flags mark spectral regions with a wide range of potential problems (e.g. no observation, poor calibration, bad pixel, or sky lines) and, therefore, must be excluded from spectral synthesis. The following individual flags provided by the survey were adopted in this work: 'Bad pixel within 3 pixels of trace', 'Pixel fully rejected in extraction', 'Sky level > flux + 10*(flux error)', and 'Emission line detected here' (with the later being adopted only for Starlight, for reasons previously mentioned). Choosing which individual flags to include when every pixel is important in a pixel-by-pixel code becomes an exercise in compromise. On the one hand, masking spectral regions affected by bad pixels or sky noise leads to more accurate estimates of the physical properties inferred through spectral synthesis. On the other hand, using all spectral flags computed by an automatic pipeline in any large-scale survey can lead to severe over-flagging and, thus, the removal of important spectral features. The chosen flags reflect this balancing act. Figure 1 shows an example of spectral fits from both codes for the object '52912-1424-151', a particularly noisy spectrum with S/N(λ 0 ) = 4.4. Black, red, and blue lines on the main panel represent the input, Starlight and Fado best-fit spectra, respectively. The right-hand side panels illustrate the best-fit SFHs based on the luminosity contribution L 4020 of the selected SSPs at the normalisation wavelength λ 0 = 4020 Å (top panels) and on their corresponding mass contributions (bottom panels). Although the mass distributions obtained with both codes differ very little, with most mass coming from SSPs older than 1 Gyr, Starlight best-fit solution includes more SSPs younger than 10 Myr compared to Fado when it comes to the light distribution. Given the relatively noisy nature of the observation, this type of variations are not necessarily indicative of a strong divergence between the codes. More reliable trends are found when grouping galaxies in populations with similar spectral characteristics, as explored in Sections 3 and 4.

Results
The following analysis is a comparison between the results from Fado and Starlight, particularly regarding the stellar properties of SF galaxies (e.g. mass, mean age, and mean metallicity). Although these results could be compared to the vast number of previous works that analysed the SDSS (e.g. Brinchmann et al. 2004;Cid Fernandes et al. 2005;Tojeiro et al. 2009), this endeavour is outside of the scope of this work for several reasons.
For one, this work is a continuation of CGP19 and P21, in that it follows similar methodologies specifically interested in looking into how (not) modelling the nebular continuum in SF galaxies can impact their inferred physical properties. Moreover, it is notably difficult to ascertain the specific details regarding how previous works dealt with spectral extraction and modelling (e.g. dust treatment, spectral masking, and evolutionary ingredients). This obstacle makes direct comparisons difficult and potentially misleading, since no clear route is available to quantify the assumptions behind the different methodologies.
Last but not least, there is also the issue of aperture effects and corrections. Several works have issued cautionary remarks when attempting to apply photometric-based aperture corrections to physical properties inferred through spectral synthesis (e.g. Richards et al. 2016;Gomes et al. 2016a,b;Green et al. 2017), whilst others claim that aperture-free properties such as SFR can still be inferred (e.g. Duarte Puertas et al. 2017). Since the main objective is to compare results between two different modelling approaches, no aperture corrections are necessary. Therefore, the following discussion focus only on the galaxy surface area covered by the spectroscopic fibre.

Sample selection
The analysis of the results is based on four main samples of galaxies. The Main Sample (MS) is defined by a cut in apparent magnitude in the r−band of 14.5 ≤ m r ≤ 17.77 and in redshift of 0.04 ≤ z ≤ 0.4. The first criterion takes into consideration photometric completeness of the survey, whereas the second aims to remove low-redshift interlopers while also assuring that the most prominent optical emission lines (e.g. Hβλ4861; [OIII]λ5007; [OI]λ6300; Hαλ6563; [NII]λ6583; [SII]λλ6717,6731) fall within the observed wavelength range. Duplicates where removed by selecting objects classified as galaxies ('specClass = 2') from the 'SpecObj' table created by the survey, which lists fibres categorised by 'SciencePrimary' (i.e. primary observation of the object). These criteria lead to the selection of 613592 objects (∼66% of the DR7 galaxy sample).
The Star-Forming Sample (SFS) is defined by the criteria of MS in combination with: EW(Hα) ≥ 3 Å, signal-to-noise at the normalisation wavelength of S/N(λ 0 ) ≥ 3, and the Kauffmann et al. (2003c) demarcation line. Similar to the criteria adopted by Asari et al. (2007), these conditions select for emission-line galaxies with relatively good spectral quality located in the 'SF' locus in the [NII]λ6583/Hα and [OIII]λ5007/Hβ emission-line diagnostic diagram (Baldwin, Phillips & Terlevich 1981). This selects 195479 objects (∼31.9% of MS). Moreover, the Intensively Star-Forming Sample (ISFS) is defined by the criteria of SFS with a cut of EW(Hα) ≥ 75 Å. This isolates galaxies in A&A proofs: manuscript no. FADO_SDSS_v11_-_final_version_-_publisher which the nebular continuum contribution is particularly high. This selects 8051 objects (∼1.3% of MS).
Finally, the Passive Sample (PS) is defined by the criteria of MS in combination with: EW(Hα) ≤ 0.5 Å and S/N(λ 0 ) ≥ 3. Similarly, Mateus et al. (2006) defined 'passive' or 'lineless' galaxies based on EW(Hα) and EW(Hβ) ≤ 1 Å, whereas Cid Fernandes et al. (2011) defined 'passive' as galaxies with EW(Hα) and EW(NII) ≤ 0.5. This selects 103510 objects (∼16.9% of MS). Figure 2 shows the distributions for the different samples of the signal-to-noise at the normalisation wavelength S/N(4020 Å), the flux of the Hα emission line F(Hα) and the Hα equivalent width EW(Hα). Moreover, Figure 3 displays the location of the four samples in the Baldwin, Phillips & Terlevich (1981) diagnostic diagram, whereas Figures A.1-A.8 show Fado and Starlight fit results for two randomly selected galaxies from each sample in the same format adopted in Figure 1. Figure 4 shows on the left-hand side panels the distributions of the total stellar mass ever-formed M ever ⋆ (top) and presently available M present ⋆ (bottom) (i.e. after correcting for mass loss through regular stellar evolution) for the MS and SFS galaxy populations. Results show that both codes yield very similar mass distributions. Indeed, the difference between the average total stellar mass from Starlight and Fado is ∼0.02 and 0.01 dex for both masses when it comes to the MS and SFS, respectively. Figure 5 compares the ISFS and PS distributions and shows similar trends, with the difference between the average total stellar mass from Starlight and Fado being ∼0.04 and 0.03 dex, respectively.

Total stellar mass, mean age and mean metallicity distributions
The age and metallicity distributions of the stellar populations contributing to the best-fit solution can be summarised in terms of their mean values, which can be interpreted as a first or-der moment of the overall stellar age and metallicity of a given galaxy. Following Cid Fernandes et al. (2005), the mean logarithmic stellar age weighted by light or mass can be defined as, respectively, where t i is the age of the i th stellar element (i.e. SSP, in this case) in the base library, and γ i and µ i are its light and mass relative contributions, respectively. The term µ i refers to the mass fraction that takes into consideration the amount of stellar matter returned to the interstellar medium through stellar evolution, thus being related to M present ⋆ . Congruously, assuming that Z i represents the relative metallicity contribution of the i th SSP to the best-fit solution, then the mean logarithmic stellar metallicity weighted by light and mass can be defined as, respectively, Applying these definitions, Figure 4 shows the distribution of the mean stellar age ⟨log t ⋆ ⟩ (centred panels) and mean stellar metallicity log⟨Z ⋆ ⟩ (right-hand side panels) estimated with Fado (blue) and Starlight (red). The MS distributions show that, overall, the Fado and Starlight results are once more rather similar, with a few noticeable differences: (a) the metallicity distributions of Fado are broader than that of Starlight, (b) the ⟨log t ⋆ ⟩ M absolute maxima differ by ∼0.1 dex between Fado and Starlight, and (c) the log⟨Z ⋆ ⟩ L and log⟨Z ⋆ ⟩ M absolute maxima differ by ∼0.1 dex between codes. Results for the SFS are relatively clearer: (a) the absolute maxima in the ⟨log t ⋆ ⟩ L differ by ∼0.1 dex between Fado and Starlight, with the former estimating slightly higher mean ages, (b) Fado displays a broader ⟨log t ⋆ ⟩ M distribution than Starlight, and (c) Fado also displays broader metallicity distributions than Starlight, with Starlight preferring values closer to solar metallicity.
It is worth noting that the overall MS mass and age distributions from Starlight illustrated in Figure 4 are qualitatively in agreement with Mateus et al. (2006), which analysed SDSS DR 7 with Starlight. Moreover, the broadening of the age and metallicity distributions going from Starlight to Fado in MS is also qualitatively similar to that reported by Werle et al. (2020) when going from purely-optical spectroscopy to UV photometry plus optical spectroscopy using Starlight, even though that work is focussed on early-type galaxies.
Comparing results from ISFS and PS in Figure 5 shows that: (a) Fado and Starlight mean age and metallicity distributions for the PS are very similar, with the main difference being the mass-weighted age distribution of Fado being narrower and with higher maxima in relation to Starlight, (b) Starlight lightweighted age (metallicity) distribution peaks at a lower (higher) left-hand side panels), mean stellar age ⟨log t ⋆ ⟩ in years (centred panels) and mean stellar metallicity log⟨Z ⋆ ⟩ in solar units (right-hand side panels). Blue and red colours represent results from Fado and Starlight, respectively, for the Main Sample (darker colours) and Star-Forming Sample (lighter colours). Top and bottom mass panels display, respectively, the total mass ever-formed throughout the life of the galaxy M ever ⋆ and the presently available total stellar mass M present ⋆ . Moreover, top and bottom age and metallicity panels represent their light-and mass-weighted versions, respectively. value than Fado when it comes the ISFS, and (c) mass-weighted age and metallicity ISFS distributions in both codes are much wider than their light-weighted counterparts.
Interestingly, the light-weighted age distributions for MS in Figure 4 show two peaks around ∼8.7 and 9.7 dex for both codes, which seem to visually match their SFS and PS distributions, respectively. This could be interpreted at first glance as evidence for the bimodal distribution of galaxies (e.g. Gladders & Yee 2000;Strateva et al. 2001;Blanton et al. 2003;Kauffmann et al. 2003b;Baldry et al. 2004;Mateus et al. 2006). Furthermore, the MS mass-weighted age distribution from Starlight displays two peaks around ∼9.8 and 10.1 dex which, once again, match its SFS and PS distributions. However, the MS mass-weighted age distribution from Fado shows only one prominent peak around ∼10.1 dex, which is congruent with its PS distribution. In addition, the SFS mass-weighted age distribution from Fado (a) is much broader than that of Starlight and (b) displays a conspicuous relative maximum at ∼9.1 dex of unknown origin, also present in the MS distribution.
At least two factors could be behind the differences between the codes regarding these particular features: (a) the different mathematical (e.g. noise treatment and minimisation procedure) and physical approaches of the two codes (e.g. nebular continuum treatment) and (b) the intrinsic degeneracies between the SSPs in the adopted library (e.g . Faber 1972;Worthey 1994;Cid Fernandes et al. 2005). In fact, it is worth considering that CGP19 found that nebular contamination in synthetic SF galaxies leads Starlight to favour best-fit solutions where most of the light comes from a combination of very young and very old SSPs (i.e. bimodal best-fit solutions). A dedicated study exploring, for instance, how variations of the SSPs in the stellar library (e.g. Leitherer et al. 1999;Bruzual & Charlot 2003;Sánchez-Blázquez et al. 2006;Mollá, García-Vargas & Bressan 2009) impact the distributions of these stellar properties would help to clarify the sources behind the MS distribution features observed in each code.
To put the results of Figures 4 and 5 into context, it is also worth noting that Cid Fernandes et al. (2005) showed that the original version of Starlight can recover the stellar mass, mean age and mean metallicity within ∼0.1 and ∼0.2 dex when modelling purely stellar synthetic spectra with S/N∼10. Later, Cid Fernandes et al. (2014) used simulations of CALIFA data (Sánchez et al. 2012) and found that the same estimated stellar properties can be affected by uncertainties between ∼0.1 and 0.15 dex related to noise alone and ∼0.15 and 0.25 dex when it comes to changes in the base models. Similar accuracy values were found in CGP19 regarding both Fado and Starlight. Among other factors, these biases can be attributed The average MS, SFS, ISFS and PS population values for mass, age and metallicity are gathered in Table 3.2 for each code. The differences between codes for total mass, mean age and mean metallicity can go up to ∼0.06 for MS and ∼0.12 dex for SFS, respectively. The most significant difference in MS relates to the average light-weighted age, whereas for SFS the most prominent discrepancies between codes rest on the lightweighted age and metallicity and amount to ∼0.12 and 0.08 dex, respectively. Moreover, the ISFS and PS results in this table show stellar properties differences between codes in the ranges of ∼0.02-0.17 dex for ISFS and ∼0-0.05 dex for PS. Although the main discrepancy in ISFS falls again on the average lightweighted age (∼0.17 dex), the main difference between codes for PS has to do with the average mass-weighted age (∼0.05 dex). The possible origins for these differences are explored in Subsection 3.3. Finally, the working definition for the mass for the remainder of this work is that of M present ⋆ .

Light and mass distributions
As the previous results have shown, the mean age ⟨log t ⋆ ⟩ and mean metallicity log⟨Z ⋆ ⟩ parameters are powerful tools to evaluate the stellar properties of galaxies, particularly when grouped in populations of the same type. However, these parameters also lack the detailed information required to understand some of the differences between Fado and Starlight reported in the previous subsection. The same would be true for any ad-hoc binning schema that can be adopted to downscale the best-fit population vector (PV) into more easily manageable parameters. Indeed, a better strategy to study the code differences for each galaxy population is to look into the full information encoded in the best-fit PVs regarding the light and mass contributions of each stellar element.
With this in mind, Figure 6 shows the light (top panels) and mass relative contributions (bottom panels) of each stellar element (i.e. SSP) in the adopted base library as a function age. The relative contributions of the six metallicities in the base are summed up along each age step (represented by the black squares). These results can be viewed as the light and mass SFHs of each galaxy population, with Figure A.9 showing their cumulative versions.
Several results are worth discussing. Firstly, the relative light contributions of young-to-intermediate SSPs with t < 10 8 yr increases with increasing EW(Hα) in both codes, following the Notes. Average values µ and corresponding standard deviations σ of the total stellar mass (in solar masses M ⊙ ) ever-formed M ever ⋆ and presently available M present ⋆ , mean stellar age (in years) weighted by light ⟨log t ⋆ ⟩ L and mass ⟨log t ⋆ ⟩ M , and mean stellar metallicity (in solar units Z ⊙ = 0.02) weighted by light log⟨Z ⋆ ⟩ L and mass log⟨Z ⋆ ⟩ M for the Main, Star-Forming, Intensively Star-Forming and Passive Samples, as defined in Subsection 3.1. PS→MS→SFS→ISFS sequence. The opposite trend is seen in old SSPs with t > 10 9 yr. This is to be expected, since young bright stellar populations dominate their older counterparts in terms of light output in ISFS and SFS galaxies, whereas the bulk of stellar content in PS is rather old and considerably less luminous. Secondly, light distribution results show a peak around t ∼ 10 9 yr in both codes, ranging between ∼5 (ISFS) and 25% (PS). Asari et al. (2007) reported a similar hump around 1 Gyr in the SFR when displayed as a function of stellar age. The authors carried tests with Starlight and found that hump disappears after changing the stellar library of SSPs from STELIB (Bruzual & Charlot 2003), as in this work, to MILES (Sánchez-Blázquez et al. 2006). The comparison between Fado and Starlight is particularly revealing when it comes to the relative contribution of young and intermediate aged SSPs. Results show that SSPs with ages t < 10 8 yr contribute with more light in Starlight than in Fado, which is increasingly noticeable with increasing EW(Hα) (PS→MS→SFS→ISFS). Indeed, Figure A.9 shows that the light contribution in Starlight is greater by 5.41% than in Fado at t = 10 7 yr for the ISFS, reaching 9.11% around t = 10 9 yr. This means that Starlight overestimates the relative light contribution of young stellar populations in relation to Fado (or vice versa), thus accounting for the ∼0.12 and 0.17 dex difference in ⟨log t ⋆ ⟩ L between the codes for the SFS and ISFS populations, respectively (Table 3.2).
When it comes to the mass distributions, the differences between the codes are considerably less pronounced. The main reason for this lies in the fact that, in the local universe, most stellar mass is locked into older stellar populations. This is best exemplified by the increasing importance of t > 10 10 yr SSPs with decreasing EW(Hα), following the ISFS→SFS→MS→PS sequence. Apart from the peak around t ∼ 10 9 yr already noted, it is worth observing that Starlight results in the ISFS show higher mass contributions from SSPs at t ∼ 3 · 10 8 yr than Fado. This seems to be offset by an excess of ∼5% of SSPs with t > 10 10 yr, which could explain why the average mean stellar age of Starlight for the SFS (∼9.7 dex) and ISFS (∼9.47 dex) documented in Table 3.2 is serendipitously identical to that of Fado. In general, these differences between Fado and Starlight when it comes to both the light and mass distributions are somewhat similar to those reported by Werle et al. (2019) when adopting a new version of Starlight that incorporates UV photometry.
Exploring further the information encoded in the PVs, Figure 7 shows the light (top panels) and mass relative contributions (bottom panels) of each stellar element in base as a function metallicity. In contrast to Figure 6, the relative contributions of the 25 age steps in the base are now summed along the six metallicities in the base library (represented by the black squares). Other remaining plot details are similar to those of Figure 6.
The light distributions displayed on the top panels show for both codes that the relative contributions of the lowest metallicities increase with increasing EW(Hα) following PS→MS→SFS→ISFS, with the opposite trend at the highest metallicities. The inversion point occurs somewhere between 2Z ⊙ /5 and Z ⊙ (i.e. log(Z ⋆ /Z ⊙ ) = −0.398 and 0, respectively). Moreover, the best-fit solutions from Fado have higher light contributions from the two lowest metallicities than Starlight, with SSPs with Z ⊙ /200 and Z ⊙ /50 (i.e. log(Z ⋆ /Z ⊙ ) ≃ −2.301 and -1.699, respectively) contributing ∼20-25% in Fado and ∼15-20% in Starlight. This is best displayed in Figure A.10, which shows that the difference in the relative contribution between Fado and Starlight is already of 8.05% at Z ⊙ /200 and reaches 13.51% by Z ⊙ /50 for ISFS. This indicates that Starlight overestimates the metallicity of SF galaxies in comparison to Fado (or vice versa). However, this interpretation is strongly tempered by the particularly large standard deviations observed for SFS and ISFS populations.
As a side note, CGP19 found using synthetic galaxy spectra that Starlight systematically underestimates the light-weighted mean meallicity by up to ∼0.6 dex with decreasing age of the galaxy (i.e. increasing EW(Hα)) for t < 10 9 yr, whereas the mass-weighted can be underestimated for 10 7 < t < 10 9 yr by up to ∼0.6 dex or overestimated by up to ∼0.4 dex for t < 10 7 yr. Similar results where found in P21. Although the details of these trends depend on the SFH of the models (e.g. instantaneous or continuous), adopted fitting methodology (e.g. including or excluding the emission lines and the Balmer and Paschen discontinuities) and S/N, the important fact to bear in mind is that these tests were carried out assuming a constant solar stellar metallicity of Z ⊙ = 0.02. Therefore, metallicity results detailed in Figures 4 and 7 cannot be easily compared to the results of CGP19 or P21. In contrast, the overall mean age trends documented in this work are compatible to those found in CGP19 or P21.
Finally, the bottom panels of Figure 7 shows that mass distributions differ very little between codes, in a similar vein to those in Figure 6. However, it is interesting to note the negligible contribution of two lowest metallicities to the mass when it comes to the PS galaxies, in contrast to both SFS and ISFS. Coupled with the trends observed in the mass distributions as a function age, these results point again to the idea of a bimodal distribution of galaxies (e.g. Gladders & Yee 2000;Strateva et al. 2001;Blanton et al. 2003;Kauffmann et al. 2003b;Baldry et al. 2004).

Relations between mass, mean age and mean
metallicity in star-forming galaxies The mass-metallicity relation displayed in Figure 8 shows an increasing divergence between the codes in light-weighted metallicity as EW(Hα) increases and metallicity decreases. This is illustrated by the somewhat steeper linear regression of Fado in comparison with Starlight and the difference between the average values for ISFS, SFS and PS populations. In contrast, both codes show similar results when it comes to the mass-weighted metallicity. Both trends are in keeping with the discussion of Figure 7.
The mass-age relation in Figure 9 also shows interesting results. For instance, Starlight shows a systematic light-weighted age underestimation in comparison with Fado with increasing EW(Hα), as shown by the increasing divergence between codes following the sequence PS→SFS→ISFS. This is particularly interesting since the age distributions from both codes are rather similar for SFS. At the same time, Fado shows a mass-weighted mean age distribution that is broader and with a lower absolute maximum when compared to Starlight, even though the average mean stellar ages of the ISFS and SFS samples for each code are very similar, as noted in Table 3.2 and in Subsection 3.2.
Finally, the age-metallicity relation in Figure 10 gives another perspective to the age and metallicity differences between codes seen in Figures 8 and 9. The contour lines alone clearly illustrate the wider age and metallicity distributions from Fado in comparison with Starlight. The mass-weighted linear regression lines indicate the overall SFS distributions are intrinsically different between codes, even though the average values for each population are rather similar. The Fado results in Figure  2 10 are particularly interesting in this regard, since light-weighted metallicity increases with age but its mass-weighted counterpart decreases with its corresponding age. Indeed, the distributions of ⟨log t ⋆ ⟩ M and log⟨Z ⋆ ⟩ M from Fado are significantly broader that those from Starlight (as illustrated in Figure 4), with relative maxima at ⟨log t ⋆ ⟩ M ∼ 9.1 yrs and log⟨Z ⋆ ⟩ M ∼ −2.1, weighting towards an anti-correlation between age and metallicity. This anticorrelation can be attributed, at least partly, to the well-known age-metallicity degeneracy, with young metal-rich stellar populations being indistinguishable of old metal-poor populations from the point of view of spectral synthesis (e.g. Faber 1972;O'Connell 1980;Bressan, Chiosi & Tantalo 1996;Pelat , 1998Cid Fernandes et al. 2005). Although this is especially noticeable when it comes to the SFS, the mass-weighted agemetallicity relation presented in Figure A.13 suggests that this trend is also present in the MS.

Impact of the nebular continuum
An important factor to consider in the interpretation of the results presented in Section 3 is the potential impact of the nebular continuum modelling approach in each code (or lack thereof) has on the estimated stellar properties of SF galaxies. Moreover, one wonders if this effect can be distinguish from (a) the intrinsic uncertainties (e.g. adequacy of age and metallicity coverages) and degeneracies associated with the adopted physical ingredients (e.g. SSPs, extinction, and kinematics) and (b) the different mathematical methods adopted in each code (e.g. Metropolis algorithm coupled with simulated annealing in Starlight and genetic differential evolution optimisation in Fado). With this in mind, one of the main objectives of the methodology presented in Section 2 was to focus on the impact of nebular continuum modelling in SF galaxies by reducing the number of variables in the code comparison, following a fitting strat- . Light (left-hand side panels) and mass-weighted (right-hand side panels) mean stellar age ⟨log t ⋆ ⟩ as a function of total stellar mass M ⋆ for the Star-forming Sample. Other legend details are identical to those in Fig. 8. egy similar to that of CGP19+P21. However, there are important methodological differences between this study and CGP19+P21 that prevent a straightforward comparison between works, such as: (a) the synthetic galaxies analysed in CGP19+P21 have constant solar stellar metallicity, (b) the current work includes two extra sub-solar metallicities in the adopted stellar library, and (c) the most prominent emission lines were masked in Starlight using individual spectral masks in this work, whereas CGP19+P21 adopted a general mask built from Starlight tests using SDSS observations (Cid Fernandes et al. 2005).
Notwithstanding, the question regarding the impact of the nebular continuum modelling approach on the inferred stellar properties remains. In order to address this issue from a different angle, Fado was applied to the SFS and ISFS galaxies in a purely stellar mode (c.f. Gomes & Papaderos 2017) using the same input spectra and spectral fitting setup as in Starlight. The objective is to model the spectra using only a combination of stellar components, similarly to Starlight, and compare it with the previous Fado results in which the stellar and nebular spectral continua were fitted self-consistently. . Light (left-hand side panels) and mass-weighted (right-hand side panels) mean stellar metallicity log⟨Z ⋆ ⟩ as a function of mean stellar age ⟨log t ⋆ ⟩ for the Star-Forming Sample. Other legend details are identical to those in Fig. 8. Figure A.14 compares the distributions of the stellar properties using Fado in 'purely stellar mode' (STmode) with the results presented in Section 3 for Fado in 'full-consistency mode' (FCmode) and Starlight. Several interesting results are worth noting: (a) Fado-STmode mass distributions are very similar to those of Fado-FCmode, (b) Fado-STmode light-weighted mean age distribution is close to that of Starlight for SFS and falls between Starlight and Fado-FCmode for ISFS, and (c) Fado-STmode mass-weighted metallicity distributions are slightly more skewed to higher values than Fado-FCmode for both SFS and ISFS, a trend similar to that observed in the Starlight results.
On the one hand, the fact that the Fado-STmode results are in general closer to Fado-FCmode than to Starlight suggests that the code differences documented in Section 3 are dominated by the fundamental mathematical and statistical differences between codes (e.g. different minimisation procedures), with different physical ingredients and methods (e.g. different nebular continuum modelling approaches and emission-line masking) playing a secondary role. However, it is reasonable to think that these two factors are likely interconnected, if not mutually dependent.
On the other hand, the Fado-STmode light-weighted age distributions for SFS and ISFS seem to skew from the Fado-FCmode distributions towards those of the Starlight. This is more clearly illustrated in Figure A.15, which compares the light and mass distributions of the PVs of Fado-STmode with those of Starlight and Fado-FCmode. This figure shows that Fado-STmode overestimates for ISFS the contribution of SSPs younger than ≤ 10 7 yr (10 9 ) by ∼5.74% (0.88%) in relation to Fado-FCmode, while also overestimating the contribution of SSPs with metallicities ≤ Z ⊙ /200 (Z ⊙ /50) by ∼9.68% (8.7%). These results show that the nebular continuum modelling approach impacts the Fado results similarly to the light-weighted age trends presented in Subsection 3.3 and, therefore, indirectly suggest that the SFS Starlight results are likely affected by the lack of a nebular continuum modelling recipe (even if its impact is mixed with other uncertainties inherent to population synthesis). One way to rigor-ously quantify this impact would be to apply Starlight in 'selfconsistent mode' to these samples of SF galaxies, which obviously is not an option, even considering the code version first presented in López Fernández et al. (2016).
These results highlight an obvious yet elusive idea worth considering during the development of the next generation of spectral synthesis codes. The comparison of new codes, with evermore relevant physical ingredients (e.g. self-consistent dust or AGN treatment), with their older counterparts will necessarily be increasingly complex due to the introduction of more statistical uncertainties and degeneracies between the physical ingredients. Parallel tests with corresponding increasingly physical complexity using synthetic spectral data (e.g. CGP19; P21) will continue to be essential in such an endeavour.

Potential implications
Assuming that purely stellar modelling in fact leads to an overestimation of mean stellar metallicity with decreasing mass or increasing EW(Hα), this means that the mass-metallicity relation presented in Figure 8 could become steeper with increasing redshifts as gas reservoirs are both increasingly more abundant and less chemically enriched. This raises the further doubt of how well the synthesis models can mimic the stellar content at large distances. Most commonly adopted evolutionary models are based on stellar libraries of stars in the solar vicinity (e.g. Leitherer et al. 1999;Bruzual & Charlot 2003;Vázquez & Leitherer 2005;Sánchez-Blázquez et al. 2006;Mollá, García-Vargas & Bressan 2009;Röck et al. 2016), which might not be representative of the stellar populations at high-redshifts, especially when it comes to extremely low metallicities and Population III stars (e.g. Schaerer 2002). The fact that AGN and nebular continua dilute absorption features (e.g. Koski 1978;Cid Fernandes, Storchi-Bergmann & Schmitt 1998;Moultaka & Pelat 2000;Kauffmann et al. 2003c;Vega et al. 2009;Cardoso, Gomes & Papaderos 2016 further exacerbates this problem from the point of view of spectral synthesis. However, this concern is counterbalanced with the increasingly sophisticated theoretical synthesis models (e.g. Coelho et al. 2007;Mollá, García-Vargas & Bressan 2009;Leitherer et al. 2010;Coelho 2014;Stanway & Eldridge 2019;Coelho, Bruzual & Charlot 2020) which can help bridge the gap towards better hybrid evolutionary models.
Moreover, a systematic mean stellar age underestimation when adopting a purely stellar modelling approach has repercussions for the current interpretation of the physical properties of young stellar populations. As noted by Reines et al. (2010), estimating the physical properties (e.g. age and metallicity) of star clusters through population synthesis can only be reliably accomplished when modelling both the nebular continuum and emission-lines. The same is true for galaxies with relatively high specific SFRs in the local universe, such as dwarfs (e.g. Papaderos et al. 1998;Papaderos & Östlin 2012) and 'green peas' (e.g. Izotov, Guseva & Thuan et al. 2011;Amorín et al 2012), or at higher redshifts (e.g. Zackrisson, Bergvall & Leitet 2008;Schaerer & de Barros 2009).
From a different perspective, variations in the emission-line measurements could impact SFR estimations and, thus, indirectly affect estimations of the cosmic SF history ). The fact that Fado measures emission-lines after subtracting the nebular continuum means that the EWs of the Balmer lines will always be greater than those based on a purely stellar modelling approach. At the same time, SFR estimators based on emission-line fluxes are not expected to change since the modelling (or not) of the nebular continuum does not impact the measured fluxes.

Summary and conclusions
The population synthesis code Fado  was applied to the main sample of galaxies from the SDSS (York et al. 2000) DR 7 (Abazajian et al. 2009) with the aim of re-evaluating the relations between the main stellar properties of galaxies (e.g. mass, mean age, and mean metallicity), particularly those of SF galaxies. The main reason for such re-analysis is the fact that Fado is the first publicly available population synthesis tool to self-consistently model both the stellar and nebular continua. In fact, previous studies adopted purely stellar population synthesis codes (e.g. Starlight, Cid Fernandes et al. 2005) to infer the physical properties of galaxies, regardless of their spectral and morphological type (e.g. Kauffmann et al. 2003a;Panter, Heavens & Jimenez 2003;Cid Fernandes et al. 2005;Asari et al. 2007;Tojeiro et al. 2009).
Comparing the physical properties inferred by the population synthesis codes Fado and Starlight for four distinct galaxy samples: Main Sample (i.e. general population of galaxies), Star-Forming (EW(Hα)>3), Intensively Star-Forming (EW(Hα)>75), and Passive (EW(Hα)<0.5), shows that: • Mass distributions for the different galaxy samples are similar between Fado and Starlight. • Mean age and mean metallicity distributions in SFS from Fado are broader than those of Starlight, especially when weighted by mass. Moreover, the average light-weighted age of Starlight is lower by ∼0.17 dex than Fado for ISFS galaxies, whereas the light-weighted metallicity of Starlight is ∼0.13 dex higher than Fado. • Even though both codes show very similar mass, age and metallicity distributions for the PS population of galaxies, Fado displays a narrower and higher mass-weighted age peak than Starlight.
• Cumulative light distributions of the best-fit PVs as a function of age show for ISFS that the contribution of SSPs with t < 10 7 yr is greater by 5.41% in Starlight than in Fado, reaching 9.11% for SSPs younger than t < 10 9 yr, which means that Starlight overestimates the relative light contribution of young stellar populations in comparison with Fado (or vice versa). Moreover, cumulative light distributions as a function of metallicity indicate that the light difference between Fado and Starlight for ISFS is ∼8.05% and 13.51% for Z ⋆ ≤ Z ⊙ /200 and Z ⊙ /50 (i.e. the lowest metallicities in the adopted stellar library), respectively. This means that Fado underestimates the metallicity when compared to Starlight (or vice versa). Meanwhile, mass distributions as a function of age and metallicity from both codes are very similar for all samples, except that Starlight once again fits more SSPs with ages t < 10 9 yr than Fado for both SFS and ISFS. • Comparing the different relation permutations between total mass, mean age and mean metallicity shows all previous results in a more general way: (a) light-weighted age and metallicity differences between codes increases with increasing EW(Hα) and decreasing total mass, (b) mass-weighted age and metallicity differences between codes are minimal, except for the age underestimation of Starlight in relation to Fado when it comes specifically to PS galaxies, and (c) SFS age and metallicity distributions of Fado are broader than those of Starlight.
These results indicate that the nebular continuum modelling approach significantly impacts the inferred stellar properties of SF galaxies, even if the negative effects of a purely stellar modelling approach are mixed with other uncertainties and degeneracies associated with population synthesis synthesis. For instance, this work found that the modelling of the nebular continuum with Fado yields steeper light-weighted mass-metallicity correlation and a flatter light-weighted mass-age correlation when compared to Starlight, a purely stellar population synthesis code. Among other potential implications, this means that the stellar populations of low-mass galaxies in the local universe with relatively high specific SFRs are both more metal poor and older than previously thought. These results are particularly relevant in light of future high-resolution spectroscopic surveys at higher redshifts, such as the 4MOST-4HS (de Jonh et al. 2019) and MOONS (Cirasuolo et al. 2011Maiolino et al. 2020), for which the fraction of intensively SF galaxies is expected to be higher.  (Figures A.5 & A.6) show in particular the differences between Fado and Starlight, with the relative contribution of SSPs younger than 100 Myr in Starlight being much more prominent than those in the Fado results. Figures A.9 and A.10 show the cumulative versions of Figures 6 and 7 with the light and mass contributions being summed with increasing age and metallicity, respectively. The fact that the assembly histories for the PS are practical identical between codes is an important sign that both codes estimate similar SFHs for galaxies with negligible to none nebular contribution to the continuum (EW(Hα)<0.5 Å). Figures A.11, A.12 and A.13 shows the relations between the total stellar mass, mean stellar age and mean stellar metallicity for the MS of galaxies for both codes. The contour and histograms show that the distributions for the three stellar properties in both codes are similar, with the only interesting difference being that Fado displays a slightly larger light-weighted mean metallicity distribution than Starlight.  Cumulative light (top panels) and mass (bottom panels) relative contributions of the stellar library elements as a function of their age (left-hand side panels) and metallicity (right-hand side panels). Full and dashed lines represent Fado results in 'full-consistency mode' (FCmode) and 'purely stellar mode' (STmode), respectively. Black squares on the left and right-hand side panels represent the age and metallicity coverages of the adopted stellar library, respectively. Other legend details are identical to those in Fig. 6.