© The Authors 2025

1. Introduction

Studying the importance of modelling the nebular emission from galaxies has been an active topic for a number of years (e.g. Krueger et al. 1995; Izotov et al. 1997; Papaderos et al. 1998, 2002; Anders & Fritze-v. Alvensleben 2003; Guseva et al. 2007; Schaerer & de Barros 2009; Reines et al. 2010; Atek et al. 2011; Papaderos & Östlin 2012; Cardoso et al. 2019, 2022; Gunawardhana et al. 2020; Pappalardo et al. 2021; Miranda et al. 2023). In particular, its aim has primarily been to analyse galaxies hosting very strong star-forming (SF) activity, the so-called extreme emission-line galaxies (EELGs). Among these systems, it has been shown that for blue compact dwarfs (BCDs) and HII galaxies, it is fundamental to model the stellar and nebular emission to retrieve accurate physical properties (Krueger et al. 1995; Izotov et al. 1997; Papaderos et al. 1998; Guseva et al. 2001, 2007). In addition to BCDs, this analysis has also been applied to other SF luminous compact objects, such as green peas (Izotov et al. 2011). On the modelling side, Anders & Fritze-v. Alvensleben (2003) updated previously existent evolutionary synthesis models (from Schulz et al. 2002) by including both continuum and line emission from the nebular gas. They found that for the case of low metallicity and young ages, the nebular emission can contribute as much as 60% to the observed broad-band fluxes.

The nebular contribution is highly sensitive to the age of the stellar population, showing a rapid decline as the population ages. For instance, several works have demonstrated that the equivalent width (EW) of Hα, a suggested tracer of the nebular contribution, decreases to relatively low levels around 10 Myr after a star formation burst (e.g. Leitherer et al. 1999; Cid Fernandes et al. 2011; Gomes & Papaderos 2017; Pappalardo et al. 2021). Consequently, stellar populations older than 10 Myr are expected to exhibit a relatively low nebular contribution. Beyond age, metallicity also plays a crucial role, being inversely related with the production of ionizing photons per unit mass (e.g. Weilbacher & Fritze-v. Alvensleben 2001; Cid Fernandes et al. 2011).

In recent decades, EELGs have been studied extensively from low-redshifts up to the epoch of reionisation (e.g. Salzer et al. 1989; Terlevich et al. 1991; Izotov et al. 1997; Cardamone et al. 2009; Amorín et al. 2010; Atek et al. 2011; Lagos et al. 2011; van der Wel et al. 2011; Amorín et al. 2015; Smit et al. 2015; Brunker et al. 2020; Endsley et al. 2021; Breda et al. 2022; Boyett et al. 2022; Simmonds et al. 2023; Llerena et al. 2024). These systems are generally characterised by high specific star formation rates, low stellar masses, sub-solar metallicities, and small amounts of dust (noting that these properties are similar to those of local HII galaxies). One consequence of these properties are the high EWs of emission lines. In fact, Hα, Hβ and [OIII]λ5007 EWs are commonly used to identify these systems and although there is no clear definition of an EW limit to identify EELGs, the considered limit varies in the range between hundreds and thousands of angstroms (e.g. Cardamone et al. 2009; van der Wel et al. 2011; Amorín et al. 2015; Endsley et al. 2021; Boyett et al. 2022, 2024). It should be noted that although in most cases the ionising radiation leading to the high EWs of EELGs has its origin on star formation, it can also originate in active galactic nuclei (AGNs, e.g. Cardamone et al. 2009; Davis et al. 2024; Llerena et al. 2024; Boyett et al. 2024). In the case of SF galaxies, the high EWs of EELGs point towards very prominent SF activity, which also implies a considerable nebular emission. In fact, for an EW(Hα) = 1000 Å, the nebular line emission leads to an enhancement of the r band magnitude of ∼0.7 (Papaderos et al. 2022), thus highlighting the relevance of the nebular emission in the case of galaxies with such high EW values.

Galaxies hosting such a strong SF activity must have a significant nebular emission and carefully accounting for this emission arises as a fundamental aspect of their analysis. This was one of the motivations for the development of the spectral fitting code known as Fitting Analysis using Differential evolution Optimization (FADO; Gomes & Papaderos 2017). This is a population spectral synthesis code that self-consistently fits the optical spectrum of a galaxy, considering both the stellar and nebular continuum. This tool was used to investigate both stellar and nebular emission when deriving the physical properties of galaxies, using synthetic spectra (Cardoso et al. 2019; Pappalardo et al. 2021) and also real data (Breda et al. 2022; Cardoso et al. 2022; Lagos et al. 2022; Miranda et al. 2023; Guo et al. 2024).

Breda et al. (2022) analysed a sample of ∼400 EELGs from SDSS-DR7 (York et al. 2000; Abazajian et al. 2009) with the objective of studying their physical properties and the fundamental correlations between them. In addition to these two goals, this study also clearly highlighted the importance of incorporating nebular emission when modelling these highly SF objects.

Miranda et al. (2023) examined how the additional consideration of the nebular emission affects the star-forming main sequence (SFMS). These authors used a sample of ∼180 000 SF galaxies from SDSS-DR7 with a mean redshift of 0.07. The results revealed no significant impact in the determination of the SFMS due to the extra consideration of the nebular emission. This result was attributed to the fact that, at low redshifts, normal SF galaxies (galaxies in the SFMS) do not experience a sufficiently high SF activity so that the nebular emission contributes significantly to the overall emission.

It thus seems that for galaxies hosting intense SF activity, accounting for the nebular emission is fundamental to accurately retrieve their physical properties. Conversely, for normal low-redshift SF galaxies, the nebular emission is not as relevant. Thus, it would be important to define a criterion to identify galaxies for which considering the nebular emission is indispensable. This is even of greater importance at increasingly higher redshifts, as a larger fraction of the galaxy population is expected to be experiencing a significant SF activity, when compared to galaxies in the local Universe.

The star formation history of the Universe shows that there was an epoch of rising SF activity that reached its maximum around z ∼ 1–3 and it has been steadily decreasing towards z = 0 (Madau & Dickinson 2014). Additionally, several works have studied the SFMS in a wide range of redshifts using different methods applied to data collected at various wavelengths (see Speagle et al. 2014 and Popesso et al. 2023 for a compilation of results). The results show that the slope of the SFMS does not evolve significantly with redshift. In contrast, the normalisation increases with redshift. This means that as the redshift increases, the general level of SF activity occurring in galaxies also increases.

Finally, the rest-frame EW of emission lines, which for SF galaxies is associated with the relative level of SF activity occurring within a galaxy, has also been shown to increase with redshift (e.g. Fumagalli et al. 2012; Faisst et al. 2016; Mármol-Queraltó et al. 2016; Reddy et al. 2018; Sun et al. 2023). The increase in EW with redshift can become quite significant, reaching more than one order of magnitude at z ∼ 4. It is worth noting that a trend has also been observed between EW and stellar mass, with lower stellar mass galaxies having higher EWs (Salzer et al. 1989; Krueger et al. 1995; Fumagalli et al. 2012). This is a consequence of the fact that the ratio between the present and the average past star-formation rate (SFR) increases with decreasing stellar mass (Brinchmann et al. 2004). Hence, the EW value that is reached considering the redshift evolution depends on the stellar mass of the galaxy.

All these points suggest that at higher redshifts, the importance of modelling the nebular emission should be even greater, due to the higher level of SF activity occurring within galaxies. In fact, this warning has been highlighted in some works (e.g. Schaerer & de Barros 2009, 2010; Atek et al. 2011). Furthermore, this can also be inferred from the evolution with redshift of the fraction of EELGs (e.g. Boyett et al. 2024). Hence, it is decisive to understand the role of nebular emission especially at higher redshifts, in order to identify for which galaxies modelling the nebular emission is fundamental and how the fraction of these objects evolves with redshift.

The importance of taking into consideration the nebular emission has also been the subject of interest due to the revolutionary observations of high-redshift galaxies by the James Webb Space Telescope (JWST). Namely, studies have highlighted the substantial importance of incorporating nebular emission in templates used to select and characterise the physical properties of high-redshift galaxies (e.g. Larson et al. 2023; McKinney et al. 2023; Wang et al. 2024). Also, JWST has enabled the census and characterisation of the physical properties of EELGs at z > 3 (Davis et al. 2024; Boyett et al. 2024; Llerena et al. 2024). The results have shown highly SF systems, which may have had a crucial role in the reionisation of the Universe, with measured EWs beyond several hundred angstroms. Galaxies with such high EW measurements have also been observed in other studies that did not focus on EELGs (Caputi et al. 2024; Roberts-Borsani et al. 2024).

JWST observations also revealed the existence of some elusive objects whose properties have been difficult to explain in light of our current knowledge. There have been several works addressing these objects and, in some cases, models with significant nebular emission have been used to explain their properties. For example, the recent discovery of the Lyα-emitting galaxy JADES-GS+53.12175−27.79763 at z = 5.943 (Saxena et al. 2024), has sparked debate regarding the origin of its emission. Cameron et al. (2024) proposes that the nebular lines and UV continuum emission of this galaxy can be explained by an initial mass function (IMF) that is between 10 to 30 times more top heavy than what is usually assumed. On the other hand, Li et al. (2024) proposes that a combination of a young and metal poor stellar population, plus a low-luminosity active galactic nucleus can also reproduce the observed emission. The different explanations lead to significantly different contributions of the nebular emission to the observed UV flux, with more than 80% for the first case and 24% for the latter case. JWST has also revealed a population of compact and extremely red objects, named ‘little red dots’ (Labbé et al. 2023; Matthee et al. 2024), whose nature is still unclear. However, models with significant nebular continuum emission have been proposed to explain the peculiar features of these objects (e.g. Iani et al. 2024; Pérez-González et al. 2024). All in all, the groundbreaking high-quality data from JWST is defying our understanding of galaxy and stellar population models and the role of the nebular emission is pivotal to explore these matters.

The significance of our study is manifold. Primarily, it will enable the identification of the population of galaxies for which the nebular emission plays a significant role at different redshifts. Moreover, some low-redshift EELGS, specifically extremely metal-poor BCDs, are local analogues of high-redshift galaxies. Therefore, accurately characterising their properties is fundamental to further comprehend the high-redshift Universe. Also, it is essential to understand how the fraction of galaxies for which the nebular emission is significant evolves with redshfit. Due to recent observation of numerous high-redshift galaxies by JWST, it has been pointed out that it is necessary to understand what are the best templates and procedures to model these objects. This work could also provide valuable insights to this question, specifically to the matter of the importance of taking into consideration the nebular emission when carrying out spectral fitting. These questions will become even more fundamental to address in light of future high-redshift surveys, such as MOONS (Multi-Object Optical and Near-infrared Spectrograph; Cirasuolo & MOONS Consortium 2016; Cirasuolo et al. 2020; Maiolino et al. 2020).

In this work, we relate the contribution of the nebular to the total optical emission (hereafter referred to simply as the nebular contribution) with some commonly used tracers, such as the EW of Balmer lines, sum of emission line fluxes, SFR, and specific SFR (sSFR = SFR/M_*). Then, we study how the nebular contribution impacts the determination of physical properties of galaxies, specifically, the mass, age and stellar metallicity. Finally, we estimate the evolution with redshift of the fraction of galaxies for which we ought to take into account the nebular emission when carrying out spectral modelling.

The paper is organised as follows. Sect. 2 explains the sample selection. In Sect. 3, we describe the spectral fitting tools used to obtain the physical properties and nebular contribution estimates. Section 4 presents the estimation of the nebular contribution and its relationship with possible tracers and some physical properties. In Sect. 5, we study the impact of considering both stellar and nebular continuum on the physical properties of galaxies. In Sect. 6, we discuss the implications of these results. Lastly, in Sect. 7 the main results and implications of this work are summarised. We use a cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_M = 0.3, and Ω_Λ = 0.7, and we consider a Chabrier (2003) IMF to estimate the star formation rate.

2. Sample

To select our sample, we considered the results from the application of FADO to the data in SDSS-DR7 from Cardoso et al. (2022). This work provides measurements of emission line properties, plus estimates of physical and evolutionary properties for 926 246 galaxies from SDSS-DR7. The spectral fitting procedure used to obtain this data set is fully described in Cardoso et al. (2022).

To define our sample, we started by selecting galaxies that were properly fit with FADO by selecting the galaxies within the 0 < χ² < 1.5 interval. Then, to ensure that the ionising radiation is the product of star formation (and not of an AGN), we selected galaxies classified as SF following the BPT diagram (Baldwin et al. 1981). We also ensured that, for the selected galaxies, FADO detected the following emission lines: [OII] (λ = 3727, 3729 Å), [OIII] (λ = 4363, 4959, 5007 Å), Hβ, HeI (λ = 5876 Å), [NII] (λ = 6548, 6584 Å), Hα, and [SII] (λ = 6717, 6731 Å). In the end, we obtained a sample of 173 457 galaxies, 19% of the original sample.

Our aim was to obtain a sample that spans a wide range of values of both SF activity and nebular contribution. The EWs of Balmer lines are used in the literature as a tracer of both these quantities (e.g. Papaderos et al. 1998; Guseva et al. 2007; Atek et al. 2011). Hence, we considered the EW(Hα) as a tracer of the properties mentioned above to select our sample.

We started by dividing our sample into bins of EW(Hα), as shown in the left panel of Fig. 1 with the blue bars representing the distribution of our selected SF sample. The sample is highly skewed, with ≈95% of the sample in the EW(Hα) < 100 Å interval. However, for our study, we wanted a wide coverage of EW(Hα) values, from the lowest to the most extreme. Since for the highest EW(Hα) values there are not many galaxies in each bin, we included all of them in our final sample. On the other hand, for lower EW(Hα) bins, the number of galaxies is large; thus we randomly selected some of these objects in order to end up with a well-rounded sample of 500 galaxies evenly distributed in bins of EW(Hα). Figure 1 shows, for our final sample, the EW(Hα) distribution (orange bars in the left panel) and the redshift distribution (right panel).

Fig. 1. Characterisation of the samples considered in this work. Left panel: EW(Hα) distribution for the selected SF galaxies, in blue, and for the final sample, in orange. Right panel: redshift distribution of the final sample.

There are two important points regarding our final sample that are worth discussing: 1) it does not follow the overall distribution of the selected SF galaxies and 2) it has a small number of galaxies compared to the selected SF galaxies. Regarding the first point, this is consequence of the objective of this work. Given that we want a sample of galaxies that cover the widest possible range of SF activity and nebular contribution, there is no need to follow the distribution of the selected SF galaxies. In fact, if we had followed the overall distribution, we would have a final sample with a significant discrepancy between the number of galaxies with low and high EW(Hα) without changing the range covered by the sample. In relation to the second point, the size of our sample is limited by the number of galaxies available for the highest EW(Hα) bins. Additionally, since we are already selecting all of the available galaxies for the highest EW(Hα) bins, increasing the size of the sample would mean increasing the number of galaxies with low EW(Hα) values. As for the previous point, this would not increase the range covered by the sample.

3. Method

In general, the study of the importance of modelling the nebular emission when estimating the physical properties of galaxies has been put in terms of the level of SF activity of the galaxies. However, an absolute criterion required to characterise a galaxy as hosting an intense or low SF activity could be ambiguous. Firstly, because this is usually a relative characterisation which compares the properties of a galaxy relative to a reference value. For example, this characterisation can be done by comparing the position of a galaxy in the star formation rate-stellar mass diagram with the reference SFMS at that redshift (e.g. Popesso et al. 2023). Secondly, the reference point used to compare the properties of a galaxy may change with redshift. In the case of the SFMS, it is well established the increase in the normalisation with increasing redshift. In this way, defining the necessity of modelling the nebular emission in terms of the nebular contribution could be a better approach.

We used the spectral fitting tool FADO to fit the spectra of the galaxies in our sample. FADO, with its self-consistent approach to the fitting of the optical spectrum of galaxies, provides a way to obtain a stellar and nebular model that are consistent with each other. FADO is capable of ensuring this self-consistent approach to spectral fitting through the combination of its coherent modelling of the most important nebular features, such as the Hα and Hβ luminosities and EWs and the continuum around the Balmer and Paschen jumps. The nebular spectrum itself is computed using standard photoionisation prescriptions (e.g. Krueger et al. 1995) considering the contribution of two-photon, free-free and free-bound emission (see Gomes & Papaderos 2017, for further details). FADO also offers the possibility of being used in pure stellar mode, meaning it can be invoked such that the fit is done assuming only a stellar continuum. The mathematical concept of FADO, in particular the search for the best-fitting solution using differential evolution optimisation (Storn & Price 1996), is a further asset. Thus, FADO is an appropriate tool to obtain a model of the stellar and nebular continuum of the galaxies in our sample in order to estimate the fraction of the nebular emission.

With this work, we aim to study how considering the nebular emission from galaxies impacts the estimation of their physical properties. To do so, we also need a spectral fitting code that lacks the self-consistency concept of FADO (i.e. a code driven solely by χ² minimization between model and observed continuum) and invoked such as to only take stellar emission into account. We used STARLIGHT¹ (Cid Fernandes et al. 2005), as representative of this type of codes. In our application, STARLIGHT only considers a linear combination of simple stellar populations (SSPs) to create a model spectrum which is fitted to the observed spectrum using a minimisation process carried out through the combination of simulated annealing with the Metropolis and Markov chain Monte Carlo (MCMC) techniques. These two codes have been compared in the past with the aim of understanding the role of nebular emission in spectral fitting (Gomes & Papaderos 2017; Cardoso et al. 2019, 2022; Breda et al. 2022), thereby supporting our methodology.

Before the spectral fitting, we corrected the spectra for Galactic extinction using the dust maps from Schlegel et al. (1998), plus the correction factor from Schlafly & Finkbeiner (2011), and we considered the Cardelli et al. (1989) extinction curve. Then, we converted the spectra to the rest-frame and rebinned it such that Δλ = 1 Å. This rebinning was performed since it is recommended for STARLIGHT. Although not required for FADO, the spectra needed to be resampled so that both codes could be applied consistently.

For the application of FADO and STARLIGHT, we considered the λ = 3000–9000 Å spectral range with a spectral basis composed of 171 SSPs, derived from Bruzual & Charlot (2003) considering 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), and the Calzetti et al. (2000) extinction law for dust attenuation corrections. The spectral basis contains 57 ages (t = [0.5 Myr, 13 Gyr]) and 3 metallicities (Z = 0.2, 0.4, 1 Z_⊙). Additionally, we allowed in the fitting the V-band extinction to vary between 0 and 4 mag. Both the spectral basis and the V-band extinction interval are different relatively to the spectral fitting parameters used in Cardoso et al. (2022). Regarding the masking of the most prominent emission lines, FADO performs this process automatically as part of its internal procedure; whereas for STARLIGHT, these spectral regions were manually identified and provided as input for the fitting. The masked regions are consistent between both codes.

Regarding the selected spectral basis, we chose a set of reasonable SSPs for the modelling of galaxies with both low and significant nebular emission. Moreover, since young SSPs contribute more significantly to the overall nebular emission, we considered all the SSPs available in Bruzual & Charlot (2003) in the critical period between 0.5–15 Myr. To ensure that our conclusions do not depend on the chosen spectral basis, we repeated the analysis using a spectral basis with 150 SSPs containing 25 ages (t = [1 Myr, 15 Gyr]) and 6 metallicities (Z = 0.005, 0.02, 0.2, 0.4, 1, 2.5 Z_⊙). The obtained results are comparable using both spectral bases and we can confirm that the results derived in this work are independent of the chosen spectral basis.

With the aim of obtaining an adequate error estimation for the parameters considered in this work, we perturbed the observed spectrum and repeated the FADO run 50 times for each galaxy. To perturb the spectrum, we considered that the error associated with each flux measurement followed a Gaussian distribution centred around zero with standard deviation equal to the error in the flux measurement. After perturbing the spectra and repeating the FADO run, we consider the error as the standard deviation of the values obtained for the runs performed for each galaxy.

As products of the application of FADO, we will use the continuum models for the total, stellar and nebular components. We will also use the flux and EW measurements for several emission lines and the estimated physical properties of galaxies, specifically the currently available and total ever formed stellar mass (M_curr and M_ever, respectively), light- and mass-weighted stellar ages (t_L and t_M, respectively) and light- and mass-weighted stellar metallicities (Z_L and Z_M, respectively). From the application of STARLIGHT, we will use the aforementioned physical properties of galaxies.

One important aspect to take into consideration is the fact that SDSS is a fibre-based survey. This means that in some cases the fibre aperture (3 arcsec) might not be sufficient to completely cover the galaxy. Since the aim of this work is to understand the impact of considering the nebular emission in spectral modelling, and not to accurately estimate the physical properties of galaxies as a whole, we do not correct for aperture effects. Additionally, even if we corrected the spectra for aperture effects, in the end we would be applying FADO and STARLIGHT to the same spectrum and comparing the different estimates obtained by the two codes. Thus, the results presented in this work do not depend on aperture effects. Henceforth, when we mention galaxy properties we are referring to the properties of the galaxy region covered by the fibre.

Another aspect that is worth highlighting is that the EW measurements of FADO are obtained following a different continuum definition in relation to what is commonly used in the literature. FADO measures the EW considering that the continuum is the monochromatic continuum level at the central wavelength of the line, thus accounting for possible absorption features and stellar velocity dispersion. However, many works in the literature measure the EW based on a continuum that is estimated considering the average level of the continuum adjacent to the considered spectral line, which is called the pseudo-continuum. In Miranda et al. (2023), these two approaches to estimate the EW were compared using FADO and data from the MPA-JHU catalogue. It was found that the approach followed in FADO to measure EWs leads to estimates that are on average 0.1 dex (∼25%) higher relative to the other approach. Hence, EW measurements obtained based on the pseudo-continuum level need to be scaled by this factor to be comparable with the relations and results derived in this work, since we use the EW estimates from FADO.

4. Nebular contribution, its tracers, and the physical properties of galaxies

The upper panel of Fig. 2 shows the observed spectrum (black line), and the total, stellar and nebular continuum (blue, green and red lines, respectively) fitted by FADO, for one of the galaxies in our sample. In the bottom panel, the nebular contribution as a function of wavelength is plotted for that same galaxy.

Fig. 2. Model fitted by FADO to a galaxy in our sample with Xneb = 10%. Upper panel: observed spectrum (black line) and total, stellar and nebular continuum fitted by FADO (blue, green and red lines, respectively). The shaded regions represent the uncertainty in the estimated models. Bottom panel: ratio between the nebular and total continuum fitted by FADO as a function of wavelength.

The definition of median nebular contribution (X_neb) used throughout this work is obtained by computing the ratio between the nebular and total continuum at each wavelength, and then calculating the median value of this ratio in the 3000–9000 Å interval. It is worth considering that the contribution of the nebular to the total continuum changes significantly with wavelength (see the bottom panel of Fig. 2 and also Izotov et al. 2011). This is particularly noticeable in the region before and after the Balmer and Paschen discontinuities, with a sudden and conspicuous drop in the contribution. Nonetheless, considering a median level of nebular contribution provides a way to have an estimate of the importance of the nebular emission to the total observed emission of a galaxy in the optical regime.

In Fig. 3, we show the X_neb distribution for the galaxies in our sample. We also show the nebular contribution at 3600 Å, blueward of the Balmer discontinuity, and at 8300 Å, after the Paschen discontinuity. These two wavelengths sample the regions where the nebular contribution in the optical spectral range is close to its maximum and minimum, respectively. As expected, the distribution of the nebular contribution at 3600 Å extends towards higher values relative to the median X_neb, whereas the nebular contribution at 8300 Å has more occurrences around lower values.

Fig. 3. Distribution of the estimated median nebular contribution (blue) and at 3600 Å and 8300 Å (orange and green, respectively) for the galaxies in our sample.

The galaxies in our sample have an average nebular contribution of X_neb = 8%. The distribution shows that we have an adequate coverage of the parameter space, even if many galaxies (38%) have X_neb < 5% and the number of galaxies drops dramatically as X_neb increases, with only 14% having X_neb > 15%.

4.1. Tracers of the nebular contribution

The analysis carried out in this work is highly demanding in terms of quality of the spectrum and continuum detection. This is not always possible, particularly at higher redshifts. Hence, it is important to understand which nebular contribution tracers, among the most commonly used in the literature, are more adequate, and to quantify this relation. This will enable the possibility of having an estimate of the optical nebular contribution simply by using a tracer.

We considered nebular contribution tracers related with emission line properties, namely the EW of Hα and Hβ, considering for each case only the galaxies where these lines had S/N > 3 (for Hα all galaxies meet the criteria whereas for Hβ, 480 galaxies are selected). We also considered the sum of the flux of the most prominent emission lines in the optical (Breda et al. 2022), specifically: [OII] (λ = 3727, 3729 Å), [OIII] (λ = 4363, 4959, 5007 Å), Hβ, HeI (λ = 5876 Å), [NII] (λ = 6548, 6584 Å), Hα, and [SII] (λ = 6717, 6731 Å). In this case, we did not apply the S/N > 3 restriction since some of these lines are weak and the sample would be significantly reduced. Finally, we also considered as tracers the current SFR and sSFR.

We started by correcting the fluxes of all emission lines for the intrinsic extinction using the Balmer decrement (Hα/Hβ) and following the procedure described in Miranda et al. (2023). To estimate the current SFR, we used the Hα flux as a tracer. First, we calculated the Hα luminosity and then, to obtain the SFR, divided it by the Kennicutt (1998) conversion factor, calibrated for a Chabrier (2003) IMF: η(Hα) = 10^41.31 erg s⁻¹ M${}_{\odot }^{-1}$ yr (Miranda et al. 2023). There are other conversion factors that could have been used (e.g. Wilkins et al. 2019). However, using a different calibration would change the SFR estimate for all galaxies by the same amount and either to lower or higher values. Thus, this would only impact the relations derived in this work in term of normalisation, but not in terms of the observed scatter. The sSFR is obtained by computing the ratio between the estimated SFR and the currently available stellar mass estimated by FADO.

It is noteworthy that the Kennicutt (1998) conversion factor between SFR and Hα luminosity assumes a continuous SFR over more than 80 Myr and that the metallicity of the galaxy is similar to the solar metallicity. For galaxies with sub-solar metallicities which deviate substantially from these assumptions, such as EELGs, this leads to an overestimation of the SFR (Weilbacher & Fritze-v. Alvensleben 2001; Papaderos et al. 2023). Although many of the galaxies in our sample have sub-solar metallicities, they are not extremely metal-poor. Thus, even though our SFR estimates are likely affected by errors, these should not be substantial enough to significantly affect the relationships established in this study.

In Fig. 4, X_neb is plotted as a function of different potential nebular contribution tracers. For all tracers studied, there is a positive correlation with X_neb even if for the sum of fluxes and SFR, the scatter is significant. Indeed, the Pearson correlation coefficient for these two tracers further points to a weak correlation with X_neb.

Fig. 4. Relation between Xneb and the different tracers considered in this work: EW of Hα (upper left), EW of Hβ (upper middle), sSFR (upper right), sum of the flux of the most prominent emission lines (lower left) and SFR (lower right). In the upper left of each panel, we present the Pearson correlation coefficient assuming a linear fit to the data. For the selected best tracers of Xneb, the red line is the derived best linear fit and the resulting equation is presented in the legend. The median error of the data is presented in the bottom right corner. The axes are in logarithmic scale.

Among the selected tracers of the nebular contribution, the EW of Hα and Hβ clearly stand out from the others, having a strong linear correlation and minor scatter. The EW is the ratio between the line flux and the continuum level. For SF galaxies, this means that it is also a measure of the relation between the current SFR, more closely represented by the line flux, and the average SFR throughout the lifetime of the galaxy, approximately represented by the continuum level. Therefore, this result is understandable and in line with the literature (Papaderos et al. 1998; Guseva et al. 2007; Atek et al. 2011).

Finally, the sSFR also shows a clear correlation with the nebular contribution, even though the scatter is larger than the one for the EW of Hα and Hβ. Since there is a significant scatter for the SFR, but not so large for the sSFR, this means that the stellar mass contributes significantly to the observed scatter. Galaxies of significantly different stellar masses can have the same SFR, hence the nebular emission associated with the SF activity will also be the same. However, the contribution of the nebular emission to the total emission will depend on the stellar mass, since the stellar mass is connected to the level of the stellar emission. By combining the SFR with the stellar mass, this factor is removed and consequently the scatter becomes lower.

The clear correlation between the sSFR and X_neb suggests that the former can be used as a tracer of the latter. However, there are some disadvantages associated. First, our sSFR estimates were obtained from the Hα luminosity, and might not be directly comparable with sSFR estimates obtained from other methods. Also, whereas the EW of Hα and Hβ can be directly measured, the sSFR depends on the method used to estimate both SFR and stellar mass, thus increasing the uncertainties associated with its estimate. Consequently, the sSFR is not the ideal tracer of the nebular contribution. Nevertheless, it can still be useful to provide an indication of the approximate level of the nebular contribution. For the specific case of galaxies that only have photometric data, and thus might not be possible to estimate either the Hα or Hβ EW, their sSFR can be estimated and used to have an approximate estimate of the level of nebular contribution.

In summary, from the analysed nebular contribution tracers, our results show that the EW of Hα and Hβ are the most appropriate, whereas the sSFR emerges as a possible alternative, with the previously mentioned caveats. Given the observed linear correlation between X_neb and these tracers, we applied a linear fit to the data and the results are presented in the upper panels of Fig. 4. The best-fitting equations obtained are the following:

$$\begin{array}{cc}\hfill log({\mathrm{X}}_{\mathit{neb}}[\%])& =0.98\pm 0.05\times log({\mathrm{EW}}_{\mathrm{H}\alpha }[\mathrm{\AA }])-1.77\pm 0.13,\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill log({\mathrm{X}}_{\mathit{neb}}[\%])& =1.03\pm 0.03\times log({\mathrm{EW}}_{\mathrm{H}\beta }[\mathrm{\AA }])-1.19\pm 0.05,\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill log({\mathrm{X}}_{\mathit{neb}}[\%])& =0.50\pm 0.03\times log(\mathrm{sSFR}[{\mathrm{yr}}^{-1}])+4.80\pm 0.24.\hfill \end{array}$$(3)

These equations allow for the estimation of X_neb without the need for carrying out spectral fitting. This is particularly useful for high-redshift galaxies for which (in most cases) the continuum has not been detected and the nebular emission might play a significant role.

In addition to the EW of Hα and Hβ, we also plotted X_neb as a function of the EW of other optical emission lines. Among the considered ones, [OIII]λ5007 and HeIλ5876 stand out as the emission lines with EWs that are best correlated with X_neb, probably due to being related with young massive stars ionising the gas. The [OIII]λ5007 emission line is commonly used to select EELGs and the HeIλ5876 emission line is relatively strong in low-metallicity SF galaxies, thus it is interesting that these lines show a strong relation with X_neb. Before fitting the relation, we applied the same S/N> 3 criteria as before, selecting 466 and 397 galaxies for [OIII]λ5007 and HeIλ5876, respectively. In Fig. 5, we show these correlations and the resulting fits.

Fig. 5. Relation between Xneb and EW of [OIII]λ5007 (left panel) and HeIλ5876 (right panel) in log–log scale. The red line is the derived best fit to the data and the resulting equation is presented in the legend. The median error of the data is presented in the bottom right corner.

The relation between X_neb and the EW([OIII]λ5007) is not linear, being better described by a second-order polynomial. This suggests that the connection between the EW([OIII]λ5007) and the EW of Hα and Hβ also follows a second-order polynomial. Possibly, this is related with the dependence on metallicity of the [OIII]λ5007 emission line. However, this analysis is beyond the scope of the present paper and should be continued in a future work. In the case of the EW(HeIλ5876), there is a linear relation with X_neb. The best-fitting equations obtained are as follows:

$$\begin{array}{cc}\hfill log(& {\mathrm{X}}_{\mathit{neb}}[\%])=0.08\pm 0.01\times log{({\mathrm{EW}}_{[\mathrm{OIII}]\lambda 5007}[\mathrm{\AA }])}^2\hfill \\ \hfill & +0.23\pm 0.04\times log({\mathrm{EW}}_{[\mathrm{OIII}]\lambda 5007}[\mathrm{\AA }])-0.24\pm 0.04,\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill log({\mathrm{X}}_{\mathit{neb}}[\%])& =0.96\pm 0.02\times log({\mathrm{EW}}_{\mathrm{He}\lambda 5876}[\mathrm{\AA }])-0.20\pm 0.03.\hfill \end{array}$$(5)

4.2. Physical properties of galaxies and the nebular contribution

It is also interesting and relevant to understand how the nebular contribution relates to different physical properties of galaxies, namely the stellar mass, age, and metallicity, along with the gaseous metallicity. We used the currently available stellar mass and light-weighted age and metallicity estimates from FADO. We estimate the gaseous metallicity, Z_gas, by calculating the tracer R₂₃ = ([OII]λ3727+[OIII]λ4959, 5007)/Hβ, using the flux measurements from FADO, and following the relation defined in Nakajima et al. (2022).

In Fig. 6, the relation between X_neb and the currently available stellar mass, light-weighted age and metallicity and gaseous metallicity is plotted. The stellar mass and age and gaseous metallicity show a negative correlation with X_neb. This shows that, as expected, high nebular contributions are particularly associated with young, low-mass and metal-poor galaxies. Moreover, considering the obtained relations between the nebular contribution and sSFR, stellar mass and gaseous metallicity we observe the same “downsizing” evidence reported in Mannucci et al. (2010), where high stellar mass galaxies have lower sSFR and higher gaseous metallicities.

Fig. 6. Relation between Xneb and currently available stellar mass (upper left), light-weighted stellar age (upper right), light-weighted stellar metallicity (bottom left) and gaseous metallicity (bottom right). In the upper right of each panel, we present the Pearson correlation coefficient assuming a linear fit to the data. The median error of the data is presented in the bottom left corner.

Our analysis does not reveal a clear correlation between the stellar metallicity and X_neb, which is also indicated by the low value of the Pearson correlation coefficient. For the relation between the gaseous metallicity and X_neb, the Pearson correlation coefficient points towards a rather moderate to strong correlation. However, the scatter is still large enough so that it is not an ideal tracer of X_neb. In the end, we can hypothesise that the gaseous metallicity could be used to obtain a broad indication of the nebular contribution level.

5. Nebular component impact threshold

One of the aims of this work is to derive a threshold above which the contribution of the nebular to the total emission becomes significant, and considering it leads to non-negligible differences in the derived physical properties of galaxies. For this analysis, we calculated the difference between the estimates of FADO and STARLIGHT for the different physical properties. Then, we divided our sample according to the EW(Hα) estimated by FADO, which in Sect. 4 we have shown to be a strong tracer of the nebular contribution, in the following four bins: EW(Hα) < 100 Å, 100 ≤ EW(Hα) < 500 Å, 500 ≤ EW(Hα) < 1000 Å and EW(Hα) ≥ 1000 Å.

Figures 7, 8, and 9 show the distribution of the logarithm of the ratio between FADO (FD) and STARLIGHT (ST) estimates of different physical properties of galaxies: stellar mass, age and metallicity, respectively. The distributions are normalised such that the area under the curve is equal to one, thus allowing to directly compare the distributions of the different EW(Hα) bins. The distributions will help identify when FADO obtains larger estimates relative to STARLIGHT (positive side of the plots) and the opposite case (negative side of the plots). Additionally, we show as a vertical dashed line the median value of each distribution.

Fig. 7. Distribution of the logarithmic difference between the currently available stellar mass (upper panel) and total ever formed stellar mass (bottom panel) estimated by FADO (FD) and STARLIGHT (ST), for four EW(Hα) bins: EW(Hα) < 100 Å (blue), 100 ≤ EW(Hα) < 500 Å (orange), 500 ≤ EW(Hα) < 1000 Å (green) and EW(Hα) ≥ 1000 Å (red). The vertical dashed lines represent the median value of each distribution.

Fig. 8. Same as Fig. 7, but for the light- and mass-weighted stellar age (upper and lower panels, respectively).

Table 1 shows the percentage of objects for which either FADO or STARLIGHT obtains an estimate that is significantly different relative to the estimate of the other code. This will complement the information given by the previously mentioned distributions and give further insights into how the distributions change with varying level of nebular contribution. FADO and STARLIGHT can retrieve the studied physical properties within a typical accuracy of 0.2 dex. This value is an upper limit of the tested accuracy of both codes when retrieving the analysed physical properties of galaxies (Gomes & Papaderos 2017; Cardoso et al. 2019; Pappalardo et al. 2021; Cid Fernandes et al. 2005, 2014). Hence, we consider that FADO and STARLIGHT obtain significantly different estimates when they are larger than 0.2 dex.

5.1. Stellar mass

Figure 7 shows the evolution with the level of nebular contribution of the distribution of the difference between FADO and STARLIGHT for the currently available and total ever formed stellar mass. The distributions are similar for both quantities and evolve in a comparable manner and can be analysed together.

For the first interval, EW(Hα) < 100 Å, the distributions are centred around zero, with less than 15% of galaxies being significantly different between the two codes. As we move to the higher EW(Hα) intervals, the distribution becomes broader, with an increasing number of galaxies being significantly different between the two codes, and the peak around zero becomes less prominent. However, the distribution does not broaden uniformly, it moves towards the negative side of the plot, which represents STARLIGHT obtaining larger estimates than FADO. Although for the 100 ≤ EW(Hα) < 500 Å interval there is an increase in the percentage of galaxies significantly different between the two codes in both directions, even if with a slight tendency for STARLIGHT obtaining higher estimates than FADO, in the two intervals of higher EW(Hα) a trend is clearly established. In the 500 ≤ EW(Hα) < 1000 Å and EW(Hα) ≥ 1000 Å intervals, STARLIGHT estimates a higher stellar mass for the vast majority of galaxies (more than 70%), while for only a few galaxies does it estimate lower stellar masses (less than 7%). This is evident in the distributions, particularly with the appearance of a peak in the region of STARLIGHT estimating higher values, around log(M^FD/M^ST) = −1.5 dex.

The results show that for EW(Hα) ≥ 500 Å, there is a clear trend for STARLIGHT to obtain higher stellar mass estimates relative to FADO, with a mean difference of ∼0.6 dex. Considering only the most extreme cases, EW(Hα)≥1000 Å, then the mean difference is ∼1 dex and can reach values up to 2 dex. Thus, we can identify the value EW(Hα) = 500 Å as the threshold for which a statistically significant difference between the two codes is identifiable.

5.2. Stellar age

Figure 8 shows the evolution with the level of nebular contribution of the distribution of the difference between FADO and STARLIGHT for the light- and mass-weighted stellar age. In this case, the distribution of these quantities evolves differently, therefore we will analyse them separately.

Considering the light-weighted stellar age, in the first three EW(Hα) bins, FADO estimates higher values relative to STARLIGHT for most galaxies. On the other hand, for the EW(Hα) ≥ 1000 Å interval, the percentage of objects where one of the codes obtains significantly higher estimates relative to the other is comparable. Looking at the distributions, it seems that except for the EW(Hα) < 100 Å interval, they become increasingly bimodal with increasing EW(Hα), with one of the peaks moving towards higher values and the other towards lower values. For the highest EW(Hα) bin, one of the peaks is around zero while the other is at $log(t_{\mathrm{L}}^{\mathrm{FD}}/t_{\mathrm{L}}^{\mathrm{ST}})\simeq 2$ dex. This suggests that as the EW(Hα) increases, there is a population of galaxies for which FADO and STARLIGHT estimates become compatible, while for another population FADO obtains increasingly larger estimates relative to STARLIGHT.

Considering the mass-weighted stellar age, for the first two EW(Hα) bins there is a similar percentage of objects which either one of the codes obtains significantly larger estimates relative to the other. For the next two EW(Hα) bins, STARLIGHT obtains larger estimates than FADO for at least 54% of the objects. Analysing the distributions for further insights, it is seen that for all EW(Hα) bins there is a peak around zero, showing that for a subsample of galaxies both codes obtain compatible estimates. However, as the EW(Hα) increases, a peak arises and becomes increasingly more pronounced in the region of the distribution corresponding to FADO obtaining lower estimates than STARLIGHT, with differences larger than 2 dex. This feature is also observed in Fig. 7 for the EW(Hα) ≥ 1000 Å interval.

In summary, regarding the light-weighted stellar age, for galaxies with EW(Hα) < 1000 Å there is a trend for FADO to obtain larger estimates compared to STARLIGHT, with a mean difference of ∼0.7 dex. Contrarily, for mass-weighted stellar ages, the trend is for STARLIGHT to estimate larger values when EW(Hα) ≥500 Å, obtaining a mean difference of ∼0.9 dex. However, the distribution becomes bimodal, with FADO estimating significantly higher values for a subset of galaxies. Therefore, these limits represent thresholds for a significant difference between the estimates of the two codes.

5.3. Stellar metallicity

Figure 9 shows the evolution with the level of nebular contribution of the distribution of the difference between FADO and STARLIGHT for the light- and mass-weighted stellar metallicity. Both quantities have close distributions and have an analogous evolution, so we will analyse them together.

Fig. 9. Same as Fig. 7, but for the light- and mass-weighted stellar metallicity (upper and lower panels, respectively).

For all defined EW(Hα) intervals, the percentage of galaxies for which one of the codes estimates significantly different metallicities relative to the other only reaches 50%. Nonetheless, as the EW(Hα) increases, there is a slight tendency for FADO to obtain higher metallicity estimates. In fact, for EW(Hα) ≥ 500 Å FADO obtains larger estimates for around 30% of the sample. This fact can also be seen in the distributions, which tend to shift to the part of the plot corresponding to FADO estimating higher values than STARLIGHT, and in the mean value evolution. However, a peak around zero remains in all distributions, highlighting that FADO and STARLIGHT measurements are compatible for a significant proportion of the galaxies. To further emphasise this point, we find that the mean differences between the two codes are 0.03 and 0.003 dex for the light- and mass-weighted metallicity, respectively. Even for the EW(Hα) ≥ 1000 Å interval, the mean differences only reaches 0.09 dex, well within the defined uncertainty level.

Summarising, there is a trend for FADO obtaining higher metallicities relative to STARLIGHT. However, it does not reach the point where FADO estimates higher values for most of the galaxies or that the mean differences are larger than the uncertainty level, even when considering the EW(Hα) ≥ 1000 Å interval. Hence, there is no clear threshold above which we observe a significant difference between the estimates of the two codes.

6. Discussion

With the determination of the best tracers of the nebular contribution and the study of the impact of considering or not the nebular emission for different physical properties of galaxies, we can now evaluate the importance of taking into account the nebular emission depending on its contribution to the total emission. Based on these results, we discuss the implications for future analysis of both low- and high-redshift galaxies. In Fig. 10, the logarithm of the median differences between FADO and STARLIGHT estimates of the stellar mass, age and metallicity for the considered four EW(Hα) bins are shown as function of EW(Hα) and X_neb, thus summarising the results presented in Figs. 7, 8, and 9.

Fig. 10. Logarithm of the median differences between FADO (FD) and STARLIGHT (ST) estimates of the stellar mass (diamonds), age (squares) and metallicity (triangles) for the previously considered EW(Hα) bins. Points in each bin slightly shifted for clearer view.

6.1. The importance of including nebular emission in spectral modelling

From our comparison between FADO and STARLIGHT, we observe that for both currently available and total ever formed stellar mass, STARLIGHT tends to obtain larger estimates. This effect becomes increasingly more significant with the increase of the nebular contribution and is most noticeable when EW(Hα) ≥ 500 Å, where for more than 70% of our galaxies, considering the nebular emission leads to lower stellar mass estimates, with an average difference above 0.6 dex. As the EW(Hα) increases beyond this threshold, the differences can reach up to 2 dex. Thus, we can identify the value EW(Hα) = 500 Å as a clear threshold for which the neglect of nebular emission in spectral models leads to significant differences in the derived stellar masses.

For the light-weighted stellar age, we find that for EW(Hα) < 1000 Å, FADO tends to obtain larger estimates relative to STARLIGHT, while for EW(Hα) ≥ 1000 Å the estimates are compatible. Looking to the mass-weighted stellar age, the results show that the differences between the estimates of the two codes become significant for EW(Hα) ≥ 1000 Å, with STARLIGHT obtaining larger estimates. For lower values, there is also a slight tendency for STARLIGHT to obtain larger estimates.

In this way, we can establish EW(Hα) = 1000 Å as a threshold for clear differences in the estimation of the stellar age due to neglecting the nebular emission. However, since considering the nebular emission impacts stellar mass estimates for EW(Hα) ≥ 500 Å, then the derived stellar ages could also be impacted at this level, due to a change in the fitted stellar populations. We will return to this point and further discuss it.

Relative to the light- and mass-weighted metallicity, our analysis suggests that there is no clear difference between the estimates of the two codes. Even if FADO has a tendency to estimate higher values for EW(Hα) ≥ 500 Å, these differences never become so significant that a clear impact due to considering the nebular emission can be identified. Nonetheless, since for EW(Hα) ≥ 500 Å the stellar mass estimates are impacted as consequence of accounting for the nebular emission, which in turn should impact the stellar age estimates, due to a change in the fitted stellar populations, as previously mentioned, then this change in stellar populations should also impact the estimated stellar metallicities. Moreover, the well-known degeneracy between stellar age and metallicity further suggests that above this limit, the stellar metallicity estimates must be affected as result of taking into consideration the nebular emission (see, for example, Conroy 2013, for a discussion of the influence of the age–metallicity degeneracy in spectral fitting).

The evolution of the difference between FADO and STARLIGHT estimates with varying level of nebular contribution that we observe in this work is qualitatively similar to the results from Gomes & Papaderos (2017), where FADO and STARLIGHT were applied to synthetic spectra. Our results can be understood as a consequence of the additional consideration of the nebular emission. The impact of including the nebular emission when modelling the optical spectrum of a galaxy has two main origins: 1) the relative contribution of the nebular to the total continuum and 2) the flat spectral shape of the nebular continuum (in the wavelengths between the Balmer and Paschen discontinuities). Both these points can be seen in Fig. 2. Regarding the first point, by dividing the total continuum in the contribution from the nebular and stellar continuum, then the stellar continuum will be at a lower level relative to when the nebular continuum is neglected, since in the latter case the total continuum corresponds to the stellar continuum. In relation to the second point, since the nebular continuum is flat, a spectral fitting tool that does not take it into consideration will select an older stellar population in order to account for the flat continuum. Such nebular continuum modelling effects have been pointed out in several works (e.g. Izotov et al. 2011; Gomes & Papaderos 2017; Breda et al. 2022; Cardoso et al. 2022; Miranda et al. 2023).

Considering these effects, the inclusion of the nebular emission leads to a clear impact on the stellar mass estimates. As mentioned in the previous paragraph, considering both stellar and nebular emission leads to a lower stellar continuum and the selection of younger stellar populations, which have lower mass-to-light ratios. Thus, for galaxies where the nebular emission is significant, then the stellar mass estimates will be lower.

On the other hand, for the stellar age and metallicity the interpretation is more complex. The previously mentioned impacts on the stellar continuum also imply changes on the fitted stellar populations, which in turn will lead to different estimates of the stellar age and metallicity. As a matter of fact, this effect was studied by Cardoso et al. (2022) where it was shown that considering the nebular emission leads to the estimation of higher light-weighted stellar ages, but compatible mass-weighted stellar ages. This was attributed to differences in the light and mass contributions of stellar populations younger than 1 Gyr as consequence of taking into account the nebular emission. For the light- and mass-weighted stellar metallicities, analogous results were observed. These results highlight that the extra inclusion of the nebular emission leads to different stellar populations being employed to fit the observed spectrum of a galaxy. Hence, it can be inferred that changes in the fitted stellar populations are one of the factors driving the differences in the stellar age and metallicity estimates as consequence of also considering the nebular emission.

In contrast with the results from Cardoso et al. (2022), we observe that considering the nebular emission can also significantly impact the mass-weighted stellar age estimates. This difference is related to the fact that in our work we observe these deviations for galaxies with EW(Hα) ≥ 1000 Å, whereas in Cardoso et al. (2022) the highest EW(Hα) interval is EW(Hα) ≥ 75 Å, containing very few galaxies with EW(Hα) ≥ 150 Å. In fact, if we consider our EW(Hα) < 100 Å interval, we see that our results are comparable with the ones from Cardoso et al. (2022). This difference in behaviour that arises for the galaxies with the highest nebular contributions means that the nebular contribution can reach such a level that the impact on the selected stellar populations is dramatic, thus influencing both light- and mass-weighted physical properties.

The discussion we undertake in this section is affected by two important factors that are worth discussing further: (1) the age–metallicity degeneracy and (2) the intrinsic differences between the two considered codes, FADO and STARLIGHT. The age–metallicity degeneracy relates to the fact that both old and metal rich stellar populations impact the spectra in similar ways, namely they tend to redden the spectra. With respect to the second point, it has been highlighted that the inner workings of different spectral synthesis tools might have a non-negligible effect on the derived physical properties (Breda et al. 2022; Cardoso et al. 2022; Miranda et al. 2023). Specifically, the minimisation algorithm used to find the best fitting model could play a major role in differences between spectral fitting codes.

It is foreseeable that our results are impacted by an interplay of these two factors, specifically the different codes can reach different best-fitting models due to their different fitting algorithms and this effect can be exacerbated by the age–metallicity degeneracy. For example, one of the codes could select a set of older stellar populations, thus leading to an older mean stellar age, whereas the other code could select a set of stellar population with higher metallicity, leading to a higher mean stellar metallicity. In the end, both fits can be equally good in terms of their χ², but lead to different estimates of the physical properties. Despite these difficulties, by fixing the variables related to the spectral fitting and conducting a judicious analysis of the results, namely considering a conservative error interval for the physical parameters, it is still achievable to derive conclusions.

As a consistency check of our results, we made use of the capability of FADO to be used in pure stellar mode. By running FADO in pure-stellar mode (PS) and comparing the results to the ones previously obtained for the full-consistency mode (FC), we can remove from the analysis the uncertainties related to the intrinsic differences between FADO and STARLIGHT. Furthermore, with this comparison we can also establish if the inclusion of the nebular emission is in fact the main driver of the differences obtained in Sect. 5.

To do this test, we applied FADO to our sample with the same parameters as the ones described in Sect. 3, but this time in PS mode. The results are summarised in Fig. 11, where we compare the estimates of the considered physical properties of galaxies obtained by FADO in full-consistency mode (FD_FC) and in pure-stellar mode (FD_PS). The full distribution of the physical properties investigated in Fig. 11 are reported in the appendix (Figs. A.1, A.2, and A.3).

Fig. 11. Logarithm of the median differences between FADO in full-consistency mode (FDFC) and in pure-stellar mode (FDPS) estimates of the stellar mass (diamonds), age (squares) and metallicity (triangles) for the previously considered EW(Hα) bins. Points in each bin slightly shifted for clearer view.

There is a trend for the differences between the FC and PS modes to be mainly within the uncertainties when considering the lower EW(Hα) bins, but they increase and become significant when EW(Hα) ≥ 500 Å. The stellar mass and age tend to be overestimated due to not considering the nebular emission, reaching ∼0.75 dex differences for EW(Hα) ≥ 1000 Å. In contrast, the stellar metallicity tends to be underestimated when neglecting the nebular emission, although at a lower level when compared to the other physical properties.

Furthermore, by comparing Figs. 10 and 11 we can have hints regarding the impact of the intrinsic differences between the two codes. We see that there is relatively good agreement between the two figures for all the properties except for the light-weighted stellar age. In this case, for the comparison between FADO and STARLIGHT, FADO retrieves higher estimates relatively to STARLIGHT and the difference decreases with increasing EW(Hα) and are negligible for EW(Hα) ≥ 1000 Å. On the other hand, for the comparison between FADO in FC and PS mode, there is a negligible difference between the estimates for the lower EW(Hα) bins that becomes significant with increasing EW(Hα), towards FADO in PS mode obtaining higher estimates. This points to some internal difference between FADO and STARLIGHT that particularly affects the light-weighted stellar age.

This analysis shows that, even if intrinsic differences between the two codes play a role, particularly in the light-weighted stellar age estimates, the main conclusions about the relevance of considering the nebular emission when fitting the spectra and deriving galaxy properties remain valid. Moreover, the previously identified EW(Hα) ≥ 500 Å threshold also remains applicable.

Summarising, the most evident impact of considering the nebular emission is on the derived stellar masses, leading to lower estimates. Nevertheless, an impact on the derived stellar ages and metallicities can also be inferred, considering the implications on the fitted stellar populations of a stellar spectrum at different levels and with different shapes. Thus, we reach the conclusion that for galaxies with EW(Hα) ≥ 500 Å the impact of considering the nebular emission is clear and significant. Considering Eq. (1), this threshold corresponds to X_neb ∼ 8%, meaning that for galaxies with a median nebular contribution above this value it is fundamental to model both stellar and nebular emission. Using Eqs. (2) and (3), we obtain the threshold in terms of EW(Hβ) and sSFR: EW(Hβ) ≃ 110 Å and log(sSFR [yr⁻¹]) ≃ −7.8, respectively. However, we note that it is possible that galaxies with median nebular contributions lower than 8% are impacted by the inclusion of the nebular emission, but this is at a level which is more subtle and difficult to assess.

Considering the different method used by FADO to measure the EW (see Sect. 3 for details), then we need to scale our threshold so that it is comparable to works that calculate the EW based on a pseudo-continuum estimate. Thus, the EW(Hα) = 500 Å threshold defined in this work corresponds to an EW(Hα) ≃ 375 Å threshold for works that use this different EW definition.

6.2. Potential impact at low and high-redshifts

Considering X_neb = 8% as the threshold above which there is an impact on the derived physical properties of galaxies, it is interesting to discuss its relevance for optical spectroscopic studies. Given the aforementioned evolution of the SF activity throughout the history of the Universe, considering the nebular contribution will have a different impact depending on the considered epoch.

At low redshifts, it is not expected that the nebular emission will have an impact for most galaxies, as explained in Sect. 1. Revisiting the results from Miranda et al. (2023), they observed no significant impact in the SFMS due to the extra inclusion of the nebular emission and this result was attributed to the fact that most galaxies in the sample did not host sufficiently high SF activity in order to have a significant nebular contribution. Their sample consists of ∼180 000 SF galaxies from SDSS-DR7, with a median EW(Hα) = 34 Å and 6% of the sample has EW(Hα) > 100 Å, considerably below the EW(Hα) ≥ 500 Å threshold defined in this work. In this way, we can confirm their hypothesis that for galaxies belonging to the SFMS at low redshift, the nebular contribution is not high enough so that considering it impacts significantly the derived physical properties of galaxies.

On the other hand, Breda et al. (2022) observed significant differences in the estimated physical properties for a sample of EELGs due to the considering the nebular emission. The EW(Hα) distribution of their sample has an average value of ∼400 Å, with approximately half of the sample having EW(Hα) ≳ 500 Å. Therefore, at least 50% of the galaxies in their sample have an EW(Hα) above the threshold defined in this work and the nebular emission should not be neglected.

Thus, our results are aligned with previous studies of the subject and highlight that while for most SF galaxies at low redshift the nebular emission is negligible, there is a subset for which it is fundamental to model the nebular emission to obtain accurate estimates of the physical properties. Amongst this subset of galaxies, it is worth highlighting the extremely metal-poor EELGs, often regarded as local analogues of high-redshift galaxies. Although rare at low redshifts, recent works have proposed that they are much more numerous at higher redshifts and likely to be common during the Epoch of Reionisation (Endsley et al. 2023; Boyett et al. 2024). Considering the importance of these objects in understanding galaxy evolution in the young Universe and their role in the reionisation epoch, it is fundamental to accurately characterise their properties and to do that it is fundamental to carry out an adequate modelling of both stellar and nebular emission.

As discussed in Sect. 1, the general level of SF activity of the Universe was higher in the past. Considering specifically the evolution of EW(Hα) with redshift, we can infer how the nebular contribution evolves with redshift, since we related the two quantities in Sect 4. In Faisst et al. (2016), the evolution with redshift of the rest-frame EW(Hα) was studied for a sample with a median stellar mass of log(M_*/M_⊙) ∼ 9.8, obtaining the following empirical relation: EW(Hα)∝(1 + z)^1.8 at z < 2 and ∝(1 + z)^1.3 at 2 ≤ z < 6. Figure 12 shows the redshift evolution of the rest-frame EW(Hα) for different stellar mass intervals following the empirical relation derived in Faisst et al. (2016). We used the data from Cardoso et al. (2022) to normalise the EW(Hα) evolution with redshift. We started by selecting only SF galaxies and then dividing them in bins of stellar mass. After that, for each stellar mass bin, we calculated the average EW(Hα) and redshift. These values are used to parameterise the redshift evolution of EW(Hα) for each stellar mass bin.

Fig. 12. Evolution of the rest-frame EW(Hα) with redshift for different stellar mass intervals. We considered the EW(Hα) redshift evolution from Faisst et al. (2016). We also present, on the right-side axis, the median optical nebular contribution corresponding to the EW(Hα), calculated using the relation derived in this work. The black dotted line is the EW(Hα) = 500 Å threshold defined in this work.

From Fig. 12, we can estimate the average EW(Hα) for a galaxy with a specific stellar mass as a function of redshfit. Hence, we can also calculate at what redshift a galaxy with a specific stellar mass will have, on average, an EW(Hα) above the defined threshold of 500 Å. This will mean that the derived physical properties of that galaxy are significantly impacted due to considering both stellar and nebular emission.

Since we consider the same evolution of EW(Hα) with redshift regardless of the stellar mass, the redshift at which the threshold is crossed depends on the average EW(Hα) at z ∼ 0. Hence, the lower the stellar mass of the galaxy, the lower the redshift at which the threshold is crossed. Starting with the highest stellar mass interval, we observe that it does not reach the threshold in the considered redshift range. For the other stellar mass intervals, between log(M_* [M_⊙]) = 7 and 11 in steps of 1 dex, the threshold is reached between z ∼ 2 and 6. This means that galaxies in the stellar mass range between log(M_* [M_⊙]) = 7 and 11 will reach a point in the z ∼ 2–6 interval, where on average the nebular contribution in the optical is so significant that it cannot be neglected when carrying out spectral fitting. However, it is important to note that these results only hold on average and do not exclude the possibility that a galaxy of any mass experiencing an intense episode of SF activity surpasses the defined EW(Hα) threshold, thus leading to the necessity of accounting for both stellar and nebular emission.

Considering again the specific case of EELGs, works by Boyett et al. (2024) and Llerena et al. (2024) have conducted a census and characterisation of these objects at high-redshift (z ∼ 3–9) using data from JADES (Bunker et al. 2020; Rieke 2020; Eisenstein et al. 2023) and CEERS (Finkelstein et al. 2023, 2024), respectively. In these works, the EW(Hα) of the EELGs ranges mainly between 300 and 1500 Å, but goes up to ∼2000 Å. Consequently, for almost all these galaxies, the optical nebular emission is so significant that it must be properly accounted for in order to obtain well determined physical properties.

With the extraordinary capabilities of JWST, an increasing number of objects are being found at these redshifts and even higher. In the near future these numbers will keep increasing with the onset of new high resolution spectroscopy instruments, namely MOONS. This means that for a progressively higher number of galaxies, considering both stellar and nebular emission in a self-consistent way is fundamental to accurately recover their physical properties and thus further develop our understanding of the process through which galaxies evolve in the young Universe.

7. Summary and conclusions

The contribution of the nebular to the total optical emission has been shown to be highly relevant when modelling galaxies in periods of intense SF activity. At low redshifts, only a small fraction of galaxies host a sufficiently large SF activity to have a significant nebular contribution, the EELGs. In contrast, the increase in the general SF level of galaxies at higher redshifts hints that the nebular contribution should be relevant to take into consideration for a larger amount of galaxies. However, a threshold at which it becomes necessary to model the optical nebular contribution has not been established. Our goal is to define this threshold and discuss the implications for both low- and high-redshift galaxies.

We selected a sample of 500 SF galaxies from SDSS-DR7, spanning a varying level of SF activity based on the EW(Hα) (see Fig. 1). Next, we applied to our sample FADO, a code that self-consistently fits the stellar and nebular continuum of the optical spectrum of a galaxy, and STARLIGHT, another population synthesis code, which lacks self-consistency prescriptions and was run with a purely stellar base. With the stellar and nebular continuum fitted by FADO we estimate the median optical nebular contribution, X_neb, for our sample (see Fig. 3). With estimates of the physical properties of both codes, we study how they differ between each other depending on the median nebular contribution.

We started by evaluating different tracers of X_neb: the EW of Hα and Hβ, the sum of the fluxes of the most prominent optical emission lines, SFR and sSFR (see Fig. 4). We find that the EW of Hα and Hβ are good tracers of X_neb, having a tight positive linear correlation. The sSFR could also be used as an approximate tracer of X_neb as the relation is relatively tight; however, it is subject to caveats that need to be carefully taken into account. On the other hand, the sum of the fluxes of the most prominent optical emission lines and the SFR are not good tracers; despite showing a positive correlation with X_neb, the scatter of the relation is substantial. We also expanded this analysis to other emission lines, finding that the EW of [OIII]λ5007 and HeIλ5876 can also serve as possible tracers of X_neb (see Fig. 5).

Using the EW(Hα), we divided the sample in four intervals and studied how the differences between FADO and STARLIGHT (when estimating the physical properties of galaxies) are impacted with an increase in the value of EW(Hα) and, consequently, in the median nebular contribution (see Figs. 7, 8, and 9 and Table 1). Our results show that the self-consistent consideration of the nebular emission has a significant impact on the physical properties of galaxies when X_neb ≳ 8%, corresponding to EW(Hα) ≃ 500 Å, EW(Hβ) ≃ 110 Å and sSFR ≃ 10^−7.8 yr⁻¹. For the EW of [OIII]λ5007 and HeIλ5876, this corresponds to ∼390 Å and ∼14 Å, respectively. Given the different definition of FADO for the continuum, for works considering a pseudo-continuum, these EW values actually correspond to: EW(Hα)≃375 Å, EW(Hβ)≃80 Å, EW([OIII]λ5007)≃290 Å and EW(HeIλ5876)≃10 Å.

At the defined threshold, taking into account both stellar and nebular emission has a clear impact on the stellar mass, leading to lower estimates for more than 70% of the galaxies with an average 0.6 dex difference, rising up to 2 dex for the galaxies with the highest nebular contribution. For the stellar age and metallicity, the differences are not as straightforward as for the stellar mass. For the light-weighted stellar age, for X_neb ≲ 15% (EW(Hα) ≃ 1000 Å) there is a trend for obtaining older ages when considering the nebular emission, whereas for the mass-weighted stellar ages only for galaxies with X_neb ≳ 15% is there a significant difference when including the nebular emission, leading to younger ages. In the case of light- and mass-weighted stellar metallicities there is no clear difference on the estimates due to accounting for both stellar and nebular emission. However, the inclusion of the nebular emission and its subsequent impact on the stellar mass estimate, also affects the estimated stellar populations that constitute the galaxy. Thus, even though the consequence of taking into account the nebular emission on the stellar age and metallicity is not easily observable, it can be inferred from our results due to a change in the underlying fitted stellar populations of the galaxy.

With the aim of testing if the intrinsic difference between FADO and STARLIGHT has an impact in these results, we applied FADO in pure-stellar mode to our sample and compared to the estimates previously obtained (see Fig. 11). This analysis has demonstrated that the intrinsic differences between the two codes play a role, but the derived results regarding the importance of considering the nebular emission when fitting the spectra and deriving galaxy properties (including the thresholds mentioned above) remain valid.

Our results imply that at low redshifts, where the nebular contribution is low, most galaxies are not impacted by also considering the nebular emission. However, there is a subset of galaxies in a phase of very strong SF activity, the EELGs, for which the contribution of the nebular emission is non-negligible and impacts the inferred physical properties and stellar populations that compose the galaxy. At higher redshifts, there are several evidences for an increase of the general SF activity level of galaxies. Considering the EW(Hα) evolution with redshift, we find that galaxies with M_* = 10⁷ − 10¹¹ M_⊙ will, on average, reach and surpass the defined threshold at z ∼ 2–6 (see Fig. 12). This implies that considering the nebular emission is a vital aspect for most galaxies at these redshifts.

This work shows that accounting for both the optical stellar and nebular continuum is fundamental for EELGs at low redshifts and more generally for SF galaxies at high-redshifts, specifically at z > 2. In light of the recent JWST observations of high-redshift galaxies and future instruments aimed at exploring the young Universe (e.g. MOONS) this calls attention to a crucial aspect to take into consideration when analysing the physical and evolutionary properties of galaxies. In the future, this work should be expanded towards high-redshift galaxies to directly estimate their nebular contribution and its relevance in spectral fitting. Additionally, other commonly used spectral fitting tools that include nebular emission in their models, such as Prospector (Leja et al. 2017; Johnson et al. 2021) and Bagpipes (Carnall et al. 2018), should be included in upcoming studies of the nebular contribution of high-redshift galaxies.

Acknowledgments

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. H.M. acknowledges support from the Fundação para a Ciência e a Tecnologia (FCT) through the PhD Fellowship 2022.12891.BD. C. P. acknowledges support from DL 57/2016 (P2460) from the ‘Departamento de Física, Faculdade de Ciências da Universidade de Lisboa’. P.P. acknowledges support by FCT through Principal Investigator contract CIAAUP-092023-CTTI. R.C. acknowledges support from FCT through the Fellowship PD/BD/150455/2019 (PhD:SPACE Doctoral Network PD/00040/2012) and POCH/FSE (EC). P.L. (DOI 10.54499/dl57/2016/CP1364/CT0010) is supported by national funds through Fundação para a Ciência e a Tecnologia (FCT) and the Centro de Astrofísica da Universidade do Porto (CAUP). D.B. acknowledges support from FCT through the Fellowship UI/BD/152315/2021. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

References

Appendix A: Comparison of FADO estimates in pure-stellar and full-consistency mode

Figures A.1, A.2, and A.3 show the distribution of the logarithm of the ratio between the estimates obtained by FADO in full-consistency (FC) and pure-stellar (PS) mode of different physical properties of galaxies: stellar mass, age and metallicity, respectively. The distributions are normalised such that the area under the curve is equal to one. Additionally, in the plots, the positive values correspond to FC mode obtaining larger estimates relative to the PS mode and negative values correspond to the opposite case.

Fig. A.1. Distribution of the logarithmic difference between the currently available stellar mass (upper panel) and total ever formed stellar mass (bottom panel) estimated by FADO in full-consistency (FC) and pure-stellar (PS) mode, for four EW(Hα) bins: EW(Hα) < 100 Å (blue), 100≤EW(Hα) < 500 Å (orange), 500≤EW(Hα) < 1000 Å (green) and EW(Hα)≥1000 Å (red). The vertical dashed lines represent the median value of each distribution.

Fig. A.2. Same as Fig. A.1, but for the light- and mass-weighted stellar age (upper and lower panels, respectively).

Fig. A.3. Same as Fig. A.1, but for the light- and mass-weighted stellar metallicity (upper and lower panels, respectively).

All Tables

All Figures

	Fig. 1. Characterisation of the samples considered in this work. Left panel: EW(Hα) distribution for the selected SF galaxies, in blue, and for the final sample, in orange. Right panel: redshift distribution of the final sample.
In the text

	Fig. 2. Model fitted by FADO to a galaxy in our sample with Xneb = 10%. Upper panel: observed spectrum (black line) and total, stellar and nebular continuum fitted by FADO (blue, green and red lines, respectively). The shaded regions represent the uncertainty in the estimated models. Bottom panel: ratio between the nebular and total continuum fitted by FADO as a function of wavelength.
In the text

	Fig. 3. Distribution of the estimated median nebular contribution (blue) and at 3600 Å and 8300 Å (orange and green, respectively) for the galaxies in our sample.
In the text

Fig. 4. Relation between Xneb and the different tracers considered in this work: EW of Hα (upper left), EW of Hβ (upper middle), sSFR (upper right), sum of the flux of the most prominent emission lines (lower left) and SFR (lower right). In the upper left of each panel, we present the Pearson correlation coefficient assuming a linear fit to the data. For the selected best tracers of Xneb, the red line is the derived best linear fit and the resulting equation is presented in the legend. The median error of the data is presented in the bottom right corner. The axes are in logarithmic scale.

In the text

	Fig. 5. Relation between Xneb and EW of [OIII]λ5007 (left panel) and HeIλ5876 (right panel) in log–log scale. The red line is the derived best fit to the data and the resulting equation is presented in the legend. The median error of the data is presented in the bottom right corner.
In the text

	Fig. 6. Relation between Xneb and currently available stellar mass (upper left), light-weighted stellar age (upper right), light-weighted stellar metallicity (bottom left) and gaseous metallicity (bottom right). In the upper right of each panel, we present the Pearson correlation coefficient assuming a linear fit to the data. The median error of the data is presented in the bottom left corner.
In the text

Fig. 7. Distribution of the logarithmic difference between the currently available stellar mass (upper panel) and total ever formed stellar mass (bottom panel) estimated by FADO (FD) and STARLIGHT (ST), for four EW(Hα) bins: EW(Hα) < 100 Å (blue), 100 ≤ EW(Hα) < 500 Å (orange), 500 ≤ EW(Hα) < 1000 Å (green) and EW(Hα) ≥ 1000 Å (red). The vertical dashed lines represent the median value of each distribution.

In the text

	Fig. 8. Same as Fig. 7, but for the light- and mass-weighted stellar age (upper and lower panels, respectively).
In the text

	Fig. 9. Same as Fig. 7, but for the light- and mass-weighted stellar metallicity (upper and lower panels, respectively).
In the text

	Fig. 10. Logarithm of the median differences between FADO (FD) and STARLIGHT (ST) estimates of the stellar mass (diamonds), age (squares) and metallicity (triangles) for the previously considered EW(Hα) bins. Points in each bin slightly shifted for clearer view.
In the text

	Fig. 11. Logarithm of the median differences between FADO in full-consistency mode (FDFC) and in pure-stellar mode (FDPS) estimates of the stellar mass (diamonds), age (squares) and metallicity (triangles) for the previously considered EW(Hα) bins. Points in each bin slightly shifted for clearer view.
In the text

Fig. 12. Evolution of the rest-frame EW(Hα) with redshift for different stellar mass intervals. We considered the EW(Hα) redshift evolution from Faisst et al. (2016). We also present, on the right-side axis, the median optical nebular contribution corresponding to the EW(Hα), calculated using the relation derived in this work. The black dotted line is the EW(Hα) = 500 Å threshold defined in this work.

In the text

Fig. A.1. Distribution of the logarithmic difference between the currently available stellar mass (upper panel) and total ever formed stellar mass (bottom panel) estimated by FADO in full-consistency (FC) and pure-stellar (PS) mode, for four EW(Hα) bins: EW(Hα) < 100 Å (blue), 100≤EW(Hα) < 500 Å (orange), 500≤EW(Hα) < 1000 Å (green) and EW(Hα)≥1000 Å (red). The vertical dashed lines represent the median value of each distribution.

In the text

	Fig. A.2. Same as Fig. A.1, but for the light- and mass-weighted stellar age (upper and lower panels, respectively).
In the text

	Fig. A.3. Same as Fig. A.1, but for the light- and mass-weighted stellar metallicity (upper and lower panels, respectively).
In the text