Impact of an AGN featureless continuum on estimation of stellar population properties

TL;DR

This study investigates how the featureless AGN continuum affects the accuracy of stellar population property estimates in galaxy spectra, highlighting the need for proper modeling of AGN contributions.

Contribution

It demonstrates the significant biases introduced by neglecting AGN power-law components in spectral synthesis, emphasizing the importance of including AGN features for accurate galaxy property recovery.

Findings

Neglecting AGN continuum causes overestimation of stellar mass by ~2 dex.

Ignoring AGN leads to overestimated stellar ages and underestimated metallicities.

Proper modeling of AGN contribution reduces uncertainties in derived galaxy properties.

Abstract

The effect of the featureless power-law (PL) continuum of an active galactic nucleus (AGN) on the estimation of physical properties of galaxies with optical population spectral synthesis (PSS) remains largely unknown. With this in mind, we fit synthetic galaxy spectra representing a wide range of galaxy star formation histories (SFHs) and including distinct PL contributions of the form $F_{\nu} \propto \nu^{-\alpha}$ with the PSS code STARLIGHT to study to which extent various inferred quantities (e.g. stellar mass, mean age, and mean metallicity) match the input. The synthetic spectral energy distributions (SEDs) computed with our evolutionary spectral synthesis code include an AGN PL component with $0.5 \leq \alpha \leq 2$ and a fractional contribution $0.2 \leq x_{\mathrm{AGN}} \leq 0.8$ to the monochromatic flux at 4020 \AA. At the empirical AGN detection threshold $x_{\mathrm{AGN}}\simeq 0.26$ that we previously inferred in a pilot study on this subject, our results show that the neglect of a PL component in spectral fitting can lead to an overestimation by $\sim$2 dex in stellar mass and by up to $\sim$1 and $\sim$4 dex in the light- and mass-weighted mean stellar age, respectively, whereas the light- and mass-weighted mean stellar metallicity are underestimated by up to $\sim$0.3 and $\sim$0.6 dex, respectively. Other fitting set-ups including either a single PL or multiple PLs in the base reveal, on average, much lower unsystematic uncertainties of the order of those typically found when fitting purely stellar SEDs with stellar templates, however, reaching locally up to $\sim$1, 3 and 0.4 dex in mass, age and metallicity, respectively. Our results underscore the importance of an accurate modelling of the AGN spectral contribution in PSS fits as a minimum requirement for the recovery of the physical and evolutionary properties of stellar populations in active galaxies.

Figures

Equations

$$F^{\text{CSP}}(\lambda ,t,Z)={\int }_0^t\Psi (t-t^{\prime })F^{\text{SSP}}(\lambda ,t^{\prime },Z)\text{ d}{t}^{\prime },$$

$$F_{\lambda }=x_{\star }\frac{F_{\lambda }^{\star }}{F_{{\lambda }_0}^{\star }}+x_{\mathrm{AGN}}\frac{F_{\lambda }^{\mathrm{AGN}}}{F_{{\lambda }_0}^{\mathrm{AGN}}},$$

$$F_{\lambda }=F_{\lambda }^{\star }+\frac{x_{\mathrm{AGN}}}{x_{\star }}{F}_{{\lambda }_0}^{\star }{\left(\frac{\lambda }{{\lambda }_0}\right)}^{\alpha -2},$$

Full text

1 Introduction

The mass assembly, star formation and chemical evolution history of galaxies are fundamental to the understanding of the physical and spectrophotometric properties of galaxies in the local universe. In recent decades, numerous studies have employed spectral synthesis with the aim of addressing these issues. This technique saw its initial realisation with the works by Whipple (1935) and Baade (1944a, b). They have attempted to categorise groups of stars within local galaxies according to their spectral types and found that, predominantly, O and B stars were in the disk and G, K, and M stars were in the bulge. This was later recognised to be a signature of the different stellar populations (population I, II, and III stars) composing distinct structural components in galaxies. These works planted the seed of the modern concept for modelling the spectral continuum of galaxies as due to the linear superposition of individual star spectra or star clusters in galaxies. Historically, spectral synthesis has followed two main approaches: population spectral synthesis (PSS) and evolutionary spectral synthesis (ESS).

On the one hand, also known as population synthesis or inversion technique, PSS aims at decomposing the observed spectrum of a galaxy into its main elementary building blocks, such as individual stars and/or groups of stars of a given age, metallicity, and stellar initial mass function (IMF). Per definition, this technique yields a discretised approximation to the star formation and chemical evolution history (SFH and CEH, respectively) of a galaxy. The PSS approach saw its inception with the works by Morgan (1956), Wood (1966), and Faber (1972) and has over the years been subject to significant development (e.g. Bica 1988; Pelat 1997, 1998; Moultaka et al. 2004; Heavens, Jimenez & Lahav 2000; Heavens et al. 2004; Cid Fernandes et al. 2005; Ocvirk et al. 2006a, b; Tojeiro et al. 2007; MacArthur et al. 2009; Koleva et al. 2009). In the last two decades, spectral synthesis has transitioned from the modelling of a few observables, such as colours, fluxes, and/or equivalent widths (EWs) of absorption lines (e.g. Wood 1966; Faber 1972; Worthey 1994; Kauffmann et al. 2003) to the more powerful pixel-by-pixel ($\lambda$-by-$\lambda$) fitting technique that exploits the full information in current medium- to high-resolution spectroscopic data (e.g. Cid Fernandes et al. 2005; Koleva et al. 2009; Tojeiro et al. 2007).

On the other hand, also known as the direct approach, ESS allows for the estimation of physical properties of galaxies, such as the stellar mass, age, and metallicity, on the basis of empirically founded assumptions on the IMF, chemical evolution, and star formation rate (SFR). This approach was pioneered by Tinsley (1968) and Spinrad & Taylor (1972) and since then has been a thriving topic of research with the creation of increasingly sophisticated models (e.g. Krüger, Fritze-v. Alvensleben & Loose 1995; Fioc & Rocca-Volmerange 1997; Leitherer et al. 1999; Zackrisson et al. 2001; Bruzual & Charlot 2003; Anders et al. 2004; Le Borgne et al. 2004; Maraston 2005; Coelho et al. 2007; Mollá et al. 2009; Kotulla et al. 2009; Vazdekis et al. 2010).

Whereas there is broad consensus that these techniques have led in the past decade to a great leap forward in our understanding of the assembly history of normal galaxies, their application to active galactic nucleus (AGN) remains controversial and problematic. As a matter of fact, the chief goal of spectral modelling of AGN spectra in the 70s and 80s has been the recovery of the AGN emission after fitting and removal of the underlying stellar spectral energy distribution (SED). The main assumption behind this approach was that the underlying host galaxy was composed solely of old and metal-rich stellar populations that are in general well represented by an elliptical galaxy spectral template. For instance, the original work by Koski (1978) showed, using photometric scans of a sample of Seyferts (named after Carl Keenan Seyfert; Seyfert 1943) and narrow-line radio galaxies, that a non-stellar component could be approximated by a featureless continuum (FC) with a power-law (PL) function $F_{\nu}\propto\nu^{-\alpha}$ in the UV-optical range with $\alpha$ between 0.2 and 2.8 and an average value of $\langle\alpha\rangle=1.0\pm 0.5$. In this study, the AGN was found to provide on average $\sim$30% of the continuum at the H$\beta$ wavelength.

The accuracy of these methods for removing the underlying stellar SED has generally been limited by various sources of uncertainty (as extensively discussed in Ho, Filippenko & Sargent 1997, and more recently in Ho 2008) and, from the spectral synthesis point of view, the techniques employed were not meant to investigate, or were not capable of exploring, the SFH and CEH of galaxies. As a consequence, it might be questioned whether these studies could accurately retrieve the pure SED of the AGN, in particular for galaxies where stellar emission dominates. Even though more elaborate methods were developed for the starlight removal, such as off-nuclear spectrum subtraction within the same galaxy (Storchi-Bergmann, Baldwin & Wilson, 1993) or a linear combination of the spectra of different galaxies (Ho, Filippenko & Sargent, 1997, 2003; Storchi-Bergmann, Cid Fernandes & Schmitt, 1998), they generally did not have the reconstruction of the stellar mass assembly history of active galaxies as prime objective.

Concerning the PL contribution, it is well-established in the framework of the AGN unified model (Antonucci 1993; Netzer 2015) that it constitutes the dominant component in the SED optical continuum of type 1 AGNs, such as quasars and Seyfert 1 galaxies. However, its relative contribution to the UV-optical in type 2 AGN, such as Seyfert 2 and low-ionisation nuclear emission-line region (LINER) galaxies (Heckman 1980), remained controversial, largely owing to the limited capability of current spectral fitting codes.

Indeed, even though a substantial body of work has been devoted to this subject, a review of the literature reveals a wide set of discordant conclusions, in some cases drawn using similar methods. Part of the confusion might be attributed to poorly described or justified methodological assumptions and the assessment of the propagation of uncertainties in the physical properties of both the estimated stellar and AGN contribution, besides sample selection effects.

For example, observations of the optical stellar features, such as the Mgb band ($\sim$5200Å), have shown in some cases a certain degree of dilution in these regions by a PL component (e.g. Koski 1978; Dressler 1984; Nelson & Whittle 1995; Serote Roos et al. 1998; Cid Fernandes, Storchi-Bergmann & Schmitt 1998; Boisson et al. 2000; Moultaka & Pelat 2000; Kauffmann et al. 2003; Garcia-Rissmann et al. 2005; Vega et al. 2009). Conversely, it was pointed out that the EWs of the infrared Ca II triplet (hereafter CaT at 8498, 8542, and 8662 Å) were similar or even larger than in normal galaxies (e.g. Terlevich, Díaz & Terlevich 1990; Nelson & Whittle 1995), as they are more suitable for a stellar population and kinematical analysis. Vega et al. (2009) applied the PSS code Starlight to active galaxies and found that Seyfert 1 galaxies have EW(CaT) between $\sim$1.5 and $\sim$7.5 Å with a median value of 4.8 Å and that Seyfert 2 lie between $\sim$1.5 and $\sim$10.5 Å with a median value of 6.5 Å. Their results suggest that Seyfert 2 show no sign of dilution of EW(CaT) as compared to normal galaxies, even though optical absorption lines such as the CaII K (CaK) band at 3933 Å are much weaker than in old, bulge-like stellar populations. Evolutionary models in an EW(CaT)–EW(CaK) diagram suggest that young stellar populations are responsible for the dilution of optical lines in active galaxies. However, the authors concluded that non-stellar contributions can reach $\sim$85% in Seyfert 1, $\sim$50 % in Seyfert 2 and starburst galaxies, and $\sim$32% even for normal galaxies by applying Starlight with multiple PLs with $\alpha=1$–2.

The controversy on the PL contribution to type 2 objects remains until today. For example, in studies by Ho, Filippenko & Sargent (1995, 1997, 2003) and Ho (2008) it was argued that samples analysed by some other groups were, in general, not large enough and unbiased, preventing an accurate separation of nuclear from global properties. Even though other surveys, such as SDSS (York et al. 2000), have surpassed statistically the number of galaxies of smaller surveys (e.g. the Palomar sample; Ho, Filippenko & Sargent 1995), these larger surveys have the disadvantage of large projected aperture sizes, where the continuum emission in type 2 AGN could be significantly diluted by star formation (Kauffmann et al. 2003). This might explain why some works have found almost no PL contribution in type 2 AGN (Cid Fernandes, Storchi-Bergmann & Schmitt 1998; Schmitt, Storchi-Bergmann & Cid Fernandes 1999, while others found that a PL could account for up $\sim$90% of the optical flux (e.g. Koski 1978; Vega et al. 2009).

For instance, Cid Fernandes, Storchi-Bergmann & Schmitt (1998) have studied the stellar content of active galaxies with long-slit spectroscopy of 38 active and 4 normal galaxies with the goal of detecting the fractional contribution of a FC relative to the stellar emission. This study detected dilution of stellar absorption lines by a PL in most of the galaxies with broad-line emission, whereas almost no PL contribution was found for most of the type 2 Seyferts in their sample. Likewise, a follow-up study in Schmitt, Storchi-Bergmann & Cid Fernandes (1999) investigated the contribution by an AGN PL to the optical using flux continuum ratios and absorption-line EWs on spectra of nuclear regions of 20 local Seyfert 2 galaxies. The conclusions drawn was that the contribution from stars younger than 10 Myr and a PL rarely exceeds 5% and that the high FC contribution (up to $\sim$70%) found by Koski (1978) was likely overestimated because of the use of an elliptical template spectrum for the nuclear stellar population.

More recently, Cid Fernandes et al. (2004) applied a $\lambda$-by-$\lambda$ spectral synthesis approach with an early version of the PSS code Starlight (Cid Fernandes et al. 2005) to constrain the SFH of 79 Seyfert 2 galaxies in the 3500–5200 Å range while simultaneously fitting an additional PL continuum with varying spectral index $\alpha$. The authors found stellar populations of all ages and a PL  mean contribution of $\sim$20–30% and maximum contribution of $\sim$63% to the monochromatic flux at 4020 Å. Although the authors warned that the PL component might be due to scattered light from a hidden AGN or a young and dusty starburst, its importance was clearly documented in this PSS modelling study. Likewise, Benítez et al. (2013) found that a PL component provides up to $\sim$30% to the monochromatic luminosity at 5100 Å by applying Starlight on 10 nearby intermediate-type AGN from SDSS.

Some contradictory results were also found by Eracleous & Halpern (2001), who characterised NGC 3065 as a LINER with broad Balmer emission lines coming from an accretion disk. The authors modelled the continuum of that galaxy with a linear combination of a PL non-stellar component with starlight described by different spectra of elliptical S0 galaxies. This work found that the SED continuum of NGC 3065 can be modelled without the need for a PL if the spectrum of NGC 4339 is taken as a template for subtracting the stellar contribution, whereas other templates would imply a non-stellar component contributing up to $\sim$10% of the SED continuum at H$\alpha$ and H$\beta$ and $\la$15% at [Oii]$\lambda$3727.

As many of the results obtained by fitting a mixture of stellar and PL templates to observed spectra remain divergent and the estimation of stellar population properties with this technique is still largely uncharted territory111Whether a second and dominant blue featureless continuum (FC2) is present in type 2 AGNs due to compact starbursts (e.g. Cid Fernandes & Terlevich 1995; Heckman et al. 1995; Tran 1995), instead of coming only from the AGN synchrotron emission, is irrelevant to our work, which focusses on the recovery of the stellar population properties (e.g. SFHs and CEHs) of the host galaxy in case a PL (or AGN-like) contribution is present., it appears to be worthwhile to supplement existing work with modelling of synthetic composite spectra with stellar and AGN components of known constitution for the sake of evaluating the capability of SED fitting codes to retrieve the input physical and evolutionary characteristics of a galaxy (e.g. SFH, CEH, and AGN contribution). The very few works existing in this regard (e.g. Bon, Popović & Bon 2014; Hayward & Smith 2015; Ciesla et al. 2015; Cardoso, Gomes & Papaderos 2016) yield a variety of conclusions.

Bon, Popović & Bon (2014) simulated a mock galaxy sample of 7000 integrated spectra of Seyfert 2 galaxies and used the PSS code ULySS (Koleva et al. 2009) to recover the underlying stellar contribution and PL component. Their main results show that the stellar populations characteristics can be retrieved whenever the signal-to-noise ratio is higher than 20 and if the stellar contribution is more than 10% of the total flux.

A different methodology was adopted by Hayward & Smith (2015), where hydrodynamical simulations with full radiative transfer modelling were used to create mock SEDs with known observational and physical characteristics (e.g. photometric bands, V-band extinction, stellar mass, dust luminosity and mass, and SFR). The authors in turn applied the magphys code (da Cunha, Charlot & Elbaz 2008) with ultraviolet to millimetre photometry in an attempt to recover the stellar population properties. Overall, results showed that most physical parameters are well estimated, although it was found that in certain cases (e.g. major mergers) AGN contamination leads to a stellar mass overestimation of $\sim$0.03 dex at AGN activity peaks.

At variance to the previous studies, a strong dependence on the estimated stellar properties on the specifics of SED fitting was reported by Ciesla et al. (2015). These authors fitted broadband photometric SEDs with cigale (Noll et al. 2009) assuming exponentially decreasing and delayed SFHs, solar-metallicity Maraston (2005) SSPs, Calzetti et al. (2000) dust extinction law, and dust remission. Their results showed that exclusion of AGN optical templates from the fit can lead to the overestimation of the stellar mass by up to $\sim$150% and of the SFR by up to $\sim$300%, depending on AGN spectral type and SFH. Both overestimations increase with the AGN overall contribution to the total infrared luminosity $x^{\mathrm{IR}}_{\mathrm{AGN}}$. Inclusion of AGN templates in the fit can lead to the overestimation of the AGN light fraction up to $\sim$100%, depending on the AGN spectral shape and $x^{\mathrm{IR}}_{\mathrm{AGN}}$. However, the estimated AGN light fraction tends roughly to its true value with increasing $x^{\mathrm{IR}}_{\mathrm{AGN}}$. Quite importantly, the authors found uncertainties in mass up to $\sim$40% ($\sim$0.15 dex) and $\sim$40–50% on the SFR, depending on the AGN spectral shape and $x^{\mathrm{IR}}_{\mathrm{AGN}}$.

Our work builds upon and extends a pilot investigation of biases in optical PSS modelling of early-type active galaxies showing a high Lyman continuum (LyC) photon escape fraction presented in Cardoso, Gomes & Papaderos (2016) (hereafter CGP16). The observational motivation behind this study comes from the estimation of a high ($>$90%) LyC photon escape fraction in early-type active galaxies by Papaderos et al. (2013) and Gomes et al. (2016), which, according to our current knowledge, host a super-massive black hole in their nuclei with expected accretion-powered activity. More generally, it is conceivable that the model framework adopted in this study is applicable to any other galaxy spheroid or spheroidal component (e.g. classical bulges), where, in the absence of absorbing gas of sufficient density, the bulk of the LyC radiation from the AGN accretion disk may escape without producing in situ optical nebular emission. Motivated by such considerations, the main objective of this study is to shed further light into the assembly history of active galaxies by quantifying potential uncertainties in the recovery of their SFHs and CEHs and characteristic stellar population properties with state-of-the-art PSS.

This paper is organised as follows. Section 2 presents our library of synthetic spectra with stellar and AGN components (Subsection 2.1) and the adopted PSS modelling approach (Subsection 2.2). Sections 3 and 4 provide an overview and discussion of our results. Finally, Section 5 summarises the main findings and conclusions of this work.

2 Methodology

2.1 Synthetic spectra library

A library of synthetic SEDs was computed with the ESS code Rebetiko (Papaderos & Gomes, in prep.; hereafter PG). These SEDs are a time- and metallicity-dependent linear combination of coeval and chemically homogeneously groups of stars known as simple stellar populations (SSPs), which is a term introduced by Renzini (1981). The flux of a composite stellar population (CSP) $F^{\mathrm{CSP}}_{\lambda}(t)$ representing the overall stellar content of a galaxy of age $t$ can be computed by linearly combining SSPs of different ages and metallicities according to an assumed SFR $\Psi(t)=\mathrm{d}M_{\star}/\mathrm{d}t$,

[TABLE]

where $F^{\textrm{SSP}}(\lambda,t^{\prime},Z)$ represents the SSP spectral library as a function of the wavelength $\lambda$, age $t^{\prime}$, and metallicity $Z$ (see e.g. Tinsley 1980; Walcher et al. 2011; Cerviño 2013; Conroy 2013 for reviews on spectral synthesis).

This work adopts Bruzual & Charlot (2003) solar-metallicity ($Z_{\odot}=0.02$) SSPs with a Chabrier (2003) IMF to create synthetic CSP spectra for a fixed metallicity and a wide range of SFH parameterisations that are generally assumed to approximate those of the different Hubble types (e.g. Rocca-Volmerange & Guiderdoni 1988; Gavazzi et al. 2002; Bruzual & Charlot 2003). Figure 1 shows the adopted SFR functions: exponentially declining (red lines) with e-folding timescales of 0.001, 1, and 3 Gyr from left to right, respectively; continuous (green line), and delayed (blue lines) with different star formation peaks as a function of the look-back time (e.g. Gavazzi et al. 2002; PG).

The library comprises 716 spectra with ages between 1 Myr and 15 Gyr for each SFH complemented with a wide range of physical properties, such as SFHs weighted by light and mass and CEHs with constant solar metallicity. These purely stellar models are used in this work to provide baseline results to which active galaxy models are compared to when applying PSS.

Another library of more complex models was created by adding simple AGN continua models to the CSPs. The optical AGN FC is assumed to be well represented by a PL defined as $F_{\nu}\propto\nu^{-\alpha}$ (e.g. Oke, Neugebauer & Becklin 1970; O’Connell 1976; Koski 1978). This approximation has been widely adopted in spectral synthesis (e.g. Goerdt & Kollatschny 1998; Schmitt, Storchi-Bergmann & Cid Fernandes 1999; Kauffmann et al. 2003; Cid Fernandes et al. 2004; Moultaka 2005) and photoionisation studies (e.g. Ferland & Netzer 1983; Stasińska 1984a, b; Mathews & Ferland 1987; Veilleux & Osterbrock 1987), in which the PL index is commonly within $\alpha=0.5$–2.

The flux $F_{\lambda}$ normalised at wavelength $\lambda_{0}$ is the linear combination of purely stellar model $F^{\star}_{\lambda}$ and AGN $F^{\mathrm{AGN}}_{\lambda}$ continua,

[TABLE]

where $x_{\star}$ and $x_{\mathrm{AGN}}$ are the fractional contributions of stellar and AGN monochromatic fluxes at $\lambda_{0}=4020$ Å, respectively. Equation 2 is constrained by the normalisation condition $x_{\star}+x_{\mathrm{AGN}}=1$ at $\lambda_{0}$.

It is useful to rewrite Equation 2 as

[TABLE]

where the total stellar flux $F^{\star}_{\lambda}$ is unchanged and it is assumed that $F^{\mathrm{AGN}}_{\lambda}\propto\lambda^{\alpha-2}$, following $F^{\mathrm{AGN}}_{\nu}\propto\nu^{-\alpha}$. We adopted AGN PL continua with $\alpha=0.5$, 1, 1.5 and 2 and $x_{\mathrm{AGN}}=0.2$, 0.4, 0.6 and 0.8 to study the effects of AGN PL index and fractional contribution variations in the estimation of stellar population properties with PSS. These wide ranges assure us that no strong prior assumptions are made concerning the AGN optical shape and relative contribution. Figure 2 shows the adopted AGN PL continua with $\alpha=0.5$, 1, 1.5, and 2 corresponding to blue, green, yellow, and red lines, respectively. This figure shows that the AGN PL becomes bluer with decreasing $\alpha$.

As an example, Fig. 3 shows SEDs resulting from combining an instantaneously formed stellar population with 100 Myr and solar metallicity with AGN PLs with $\alpha$ and $x_{\mathrm{AGN}}$ variations. Purely stellar spectra are represented by black lines and models with AGN PLs with $\alpha=0.5$, 1, 1.5, and 2 for $x_{\mathrm{AGN}}=0.4$ (top panel) and $x_{\mathrm{AGN}}=0.2$, 0.4, 0.6, and 0.8 for $\alpha=1.5$ (bottom panel) are represented by blue, green, yellow, and red lines, respectively.

This figure illustrates two facts of special relevance to this study. Firstly, the normalisation wavelength $\lambda_{0}$ defines how the AGN continuum shape affects the underlying stellar continuum since Eq. 3 implies a flux increase when $\alpha<2$ larger than $F_{\star,\lambda_{0}}$ for $\lambda<\lambda_{0}$ and smaller than $F_{\star,\lambda_{0}}$ for $\lambda>\lambda_{0}$. Thus, adopting typical values of $\lambda_{0}$ such as 4020 Å (e.g. Cid Fernandes et al. 2004) or 5050 Å (e.g. Goerdt & Kollatschny 1998) means that the AGN continuum approximated as a PL for $\alpha<2$ has a slope similar to that of young SSPs. Secondly, the dilution effect of absorption-line features due to the addition of a FC (e.g. Koski 1978; Serote Roos et al. 1998; Moultaka & Pelat 2000; Kauffmann et al. 2003) depends on adopted $\lambda_{0}$, $\alpha,$ and $x_{\mathrm{AGN}}$.

Figure 4 illustrates this last point by presenting SEDs normalised at $\lambda_{0}$ for an instantaneous burst SFH with 100 Myr and solar metallicity (black line) combined with AGN continua for fixed $\alpha=1.5$ and $x_{\mathrm{AGN}}=0.2$, 0.4, 0.6, and 0.8 (blue, green, yellow, and red lines, respectively). This figure shows that $\lambda_{0}$ defines a demarcation point in the spectra that separates regions suffering from blueing or reddening. Moreover, the dilution of the underlying stellar continuum induced by the AGN PL increases with increasing $x_{\mathrm{AGN}}$.

2.2 Application of population synthesis

Both purely stellar and active galaxy SEDs were modelled with the latest public distribution222Starlightv04: http://astro.ufsc.br/starlight/ of Starlight (Cid Fernandes et al. 2005) to investigate the ability of PSS to infer the SFHs and CEHs of mock galaxies. The code Starlight is widely used (e.g. Mateus et al. 2006; Stasińska et al. 2008; Cid Fernandes et al. 2010, 2011; Kehrig et al. 2012; Benítez et al. 2013; Papaderos et al. 2013; Cid Fernandes et al. 2014; Stasińska et al. 2015; Gomes et al. 2016) and can be regarded as a good representative of the state of the art. It is important to bear in mind the main objective of this work is to quantify potential biases on fundamental physical properties of active galaxies estimated with purely-stellar PSS that might rise owing to the introduction of a simple AGN model, as in Ciesla et al. (2015).

The adopted base library comprises Bruzual & Charlot (2003) SSPs with 25 ages (illustrated as the vertical black lines on the top of the right-hand side panels of Fig. 5) and 4 metallicities ($Z=0.2$, 0.4, 1 and 2.5 $Z_{\odot}$). This library is similar to that adopted by Asari et al. (2007) with the following notable differences: first, better coverage of the younger evolutionary stages and, second, without the two lowest metallicities that are not expected to contribute significantly to the intrinsic degeneracies. A base library with a finer age and metallicity coverage would in principle help reduce uncertainties associated with poorly sampled evolutionary stages, although at a high computational cost. For instance, base libraries with 100 and 300 SSPs (i.e. the maximum number of base elements that Starlight can handle) take $\sim$110 and 1000 seconds to perform a single fit to one of these synthetic spectra in a Intel® Core™  i7 CPU 870 @ 2.93GHz workstation desktop. Thus, the adopted base library should be regarded as a compromise between computational time and a fine age and metallicity coverage. The spectral fitting is made between 3400 and 8900 Å, following common practice in local galaxy studies (e.g. Asari et al. 2007; Ribeiro et al. 2016) and to take advantage of the wavelength coverage of the STELIB stellar library (Le Borgne et al. 2003) adopted in the Bruzual & Charlot (2003) evolutionary models.

Moreover, the stellar kinematics and V-band extinction $A_{V}$ were fixed to zero and no pixel clipping method was adopted. The reason for keeping stellar kinematics and extinction fixed to input values in the Starlight fits is to facilitate an easier comparison between fits to purely stellar and active galaxy SEDs and to better isolate the potential PSS biases owing to the addition of the AGN PL continuum.

The AGN models were processed by Starlight for three fitting set-ups:

1. (i)

without any PL contribution in the base 2. (ii)

including in the base a single PL with same $\alpha$ as that embedded in the input spectra 3. (iii)

including in the base PLs with $\alpha=0.5$, 1, 1.5 and 2, leaving to Starlight the choice to pick up the PL which best matches the data.

Set-ups (i) and (ii) may be regarded as special cases in PSS modelling. The former corresponds to the worst-case scenario where fitting of an active galaxy SED is attempted with purely stellar templates, whereas the latter represents the ideal situation of the base library containing a single PL that is identical to that integrated in the input spectrum. Conversely, set-up (iii) can be viewed as a conservative fitting approach in the situation of not having prior knowledge of the AGN optical continuum shape (e.g. Ho & Kim 2009; Roche et al. 2016). Hereafter, these approaches are referred to as no power law (npl), single power law (spl), and multiple power laws (mpl). In the following, focus is given to results derived for an instantaneous burst SFH. However, similar trends are seen for the other SFHs, results for which are discussed when appropriate.

As an example, Figure 5 shows the average results for 10 Starlight fits for npl with the best fit (blue line) to the input synthetic active galaxy spectrum (red line) with an stellar population instantaneously formed with 100 Myr and solar metallicity (green line) and an AGN PL with $\alpha=1.5$ and $x_{\mathrm{AGN}}=0.2$ (yellow line) in the main panel. The bottom panel presents the residuals after subtracting the fit from the input spectrum and the right-hand side panels show the estimated SFHs in terms of the light and mass fraction (top and bottom, respectively).

3 Results

Figure 6 shows the ratio in the currently available total stellar mass $M_{\star}$ between the PSS (output) and ESS (input) values as a function of age $t$, with increasing $x_{\mathrm{AGN}}$ from left- to right-hand side panels, respectively. The results from npl, spl, and mpl are presented on the top, middle, and bottom rows, respectively (subsequent figures follow a similar panel configuration, except when explicitly stated otherwise). Moreover, PSS results for purely stellar spectra results are represented by black lines, while active galaxy spectra with $\alpha=0.5$, 1, 1.5, and 2 are represented by blue, green, yellow, and red lines, respectively. The results for purely stellar models show a maximum mass overestimation up to $\sim$0.5 dex, which are likely linked to a poor age coverage of the SSPs in the adopted base library.

The results in Fig. 6 for npl show a systematic mass overestimation up to $\sim$3.5 dex correlated with increasing $\alpha$ and $x_{\mathrm{AGN}}$ and with decreasing age $t$, with a maximum at $t\sim 10$ Myr. Figure 12 in Appendix A shows similar results for other SFHs. Results for spl and mpl reveal local mass overestimations within $\sim$1 dex. This overestimation increases with increasing $x_{\mathrm{AGN}}$, is somewhat independent from $\alpha$ and becomes more severe in models with ages where spectral synthesis of purely stellar spectra also show considerable uncertainties.

One question of considerable interest concerns the origin of this severe stellar mass overestimation. Results show that Starlight has to compensate for the lack of the PL component in the base in a purely stellar modelling configuration. This leads to a nearly bimodal mixture of young and old stellar populations, as illustrated in the SFH panels of Fig. 5. It is to be expected that in such a modelling approach the blue AGN continuum is accounted for by a non-realistic mixture of stellar populations, where very young SSPs would always dominate the light, whereas the contribution of old SSPs is almost negligible in terms of light. Notwithstanding this fact, these old SSPs provide the bulk of stellar mass, thus leading to the observed mass overestimation.

Figure 7 shows the difference in the light-weighted mean stellar age $\langle\log t_{\star}\rangle_{L}$ between PSS and ESS values as a function of age $t$. The results of purely stellar models show mean stellar age uncertainties within $\sim$0.1 dex. Cid Fernandes et al. (2005) showed that both light-weighted mean stellar age and metallicity could be recovered with Starlight within a typical uncertainty of $\sim$0.2 dex using mock data with S/N=10. This uncertainty is considered in this work to be a fiducial value.

The results in Fig. 7 for npl show a systematic age over- and underestimation of up to $\sim$1.5 and $\sim$2 dex for young ($t\lesssim 100$ Myr) and old ($t\gtrsim 100$ Myr) evolutionary models, respectively. This trend is accentuated with increasing $x_{\mathrm{AGN}}$ and originates from an increasing dilution of stellar absorption features, as seen in Fig. 4. Figure 13 shows similar results for other SFHs. Moreover, the results for spl show an age over- and underestimation of up to $\sim$0.3 dex, which appear to be correlated with ages where purely stellar fits also show considerable uncertainties. In addition, there is a clear age underestimation at $\sim$20 Myr that increases with increasing $\alpha$ and $x_{\mathrm{AGN}}$. The results for mpl show uncertainties and trends similar to those found for spl. However, the age of old evolutionary models ($t>1$ Gyr) can be underestimated by as much as $\sim$1.5 dex. This feature becomes more prominent with increasing $x_{\mathrm{AGN}}$ for $\alpha=0.5$. The results in Fig 11 suggest that this trend is due to a poor estimation of $x_{\mathrm{AGN}}$ for this fitting set-up, which also explains the slight mass underestimation for this age and AGN contribution.

Figure 8 illutrating results for mass-weighted mean stellar age $\langle\log t_{\star}\rangle_{M}$ show overall uncertainties that are larger than their light-weighted counterparts. Synthesis results of purely stellar models show a maximum age overestimation up to $\sim$2 dex. Moreover, the results for npl show a systematic age overestimation up to $\sim$4 dex with $\alpha$ and decreasing $t$. Figure 14 shows similar results for other SFHs. This trend is similar to that found in the top row of Fig. 6. Results for spl and mpl show age overestimations up to $\sim$3 and 2 dex, respectively, in ages where purely stellar models display considerable uncertainties, similarly to total stellar mass results in Fig. 6.

Figure 9 shows the difference in the light-weighted mean stellar metallicity $\log\langle Z_{\star}\rangle_{L}$ between PSS and ESS values as a function of age $t$. The results for the purely stellar models show metallicity over- and underestimations within $\sim$0.2 dex correlated with evolutionary ages with light-weighted age under- and overestimations, respectively. This is likely due to age-metallicity degeneracies that prevent a clear distinction between young metal-rich and old metal-poor stellar populations or vice versa (Faber 1972; O’Connell 1980; Renzini & Buzzoni 1986; Bressan, Chiosi & Tantalo 1996; Pelat 1997, 1998; Cid Fernandes et al. 2005). This problem arises after changes of absorption line features induced by variations in both age and metallicity (e.g. Worthey 1994), which is aggravated in this case by the dilution of the absorption features by the PL.

The results in Fig. 9 for npl show a systematic underestimation of metallicity by up to $\sim$0.7 dex with decreasing $\alpha$ and age $t$ and increasing $x_{\mathrm{AGN}}$. Figure 15 shows similar results for other SFHs. The metallicity underestimation plateau at -0.7 dex seen for $\alpha=2$ and $x_{\mathrm{AGN}}=0.8$ when $t<1$ Gyr corresponds to a metallicity of $0.2\,Z_{\odot}$, which is the lowest metallicity found in the base. Thus, this result is due to the construction of the base library and has no physical meaning. In addition, this result is due to an increasing dilution of the absorption features with increasing $x_{\mathrm{AGN}}$, which artificially makes the stellar continuum increasingly metal poorer. The results for spl display random metallicity uncertainties up to $\sim$0.15 dex predominately located around ages where purely stellar synthesis results also show considerable uncertainties. These roughly increase with increasing $\alpha$ and $x_{\mathrm{AGN}}$. Moreover, the results for mpl show an increasing metallicity overestimation up to $\sim$0.2 dex with increasing $\alpha$ and $x_{\mathrm{AGN}}$ for $t\gtrsim 1$ Gyr, mirroring the age underestimation at the same evolutionary stages shown in Fig. 7.

Figure 10 illustrating results for the mass-weighted mean stellar metallicity $\log\langle Z_{\star}\rangle_{M}$ show that light- and mass-weighted metallicities display similar trends with $t$, $\alpha$ and $x_{\mathrm{AGN}}$ for all set-ups. Indeed, the results for npl show again that the metallicity can be underestimated by up to $\sim$0.7 dex. Figure 15 shows similar results for other SFHs. Results also show that the mass-weighted properties tend to have larger uncertainties than their light-weighted counterparts. The reason for this is that small light variations translate into large mass variations when applying the mass-to-light ratio to the light fractions vector of the SSPs to determine mass related properties.

Figures 17–21 in Appendix A show plots of spectral synthesis results for a continuous SFH analogous to those presented in Figures 6–10. These results show that stellar uncertainties induced by the added AGN PL component are relatively lower for a continuous SFH than those estimated for an instantaneous burst SFH. This happens because short phases of significant spectrophotometric evolution in the first few hundred Myr are largely washed out in the case of continuous star formation when looking at the time evolution of light- and mass-weighted stellar properties. With respect to an instantaneous burst of a given age and metallicity, the SSP that best approximates this evolutionary stage is not necessarily included in the base library. In this case, and assuming the best case scenario in which there are almost no degeneracies among the stellar populations, adjacent SSPs in age and/or metallicity would have to be used, which unavoidably induces strong local deviations between the true and estimated evolutionary quantities. Thus, results obtained for the instantaneous burst SFH should be viewed as upper limits to the uncertainties on the analysed stellar properties.

4 Discussion

A substantial body of work devoted to the exploration of the assembly history of galaxies relies on automated PSS modelling of large extragalactic data sets (e.g. Kauffmann et al. 2003; Asari et al. 2007; Cid Fernandes et al. 2010). A common practice in these studies is the prior exclusion from the PSS modelling of galaxies classified as Seyfert on the basis of BPT emission-line diagnostics, since both the spectroscopic characteristics of these systems (e.g. broad emission lines and dilution of stellar features by the featureless AGN PL) and code-specific limitations generally hinder a reliable separation of the non-thermal and stellar SEDs and the extraction from the latter of key physical and evolutionary properties (e.g. stellar mass, mean stellar age and metallicity, and SFH). However, as pointed out in Papaderos et al. (2013) and Gomes et al. (2016) (see also CGP16), standard BPT classification diagnostics become inapplicable in the case of virtually gas-evacuated galaxies where the bulk of the Lyman continuum radiation from an AGN eventually escapes without locally producing detectable optical line emission. This, together with dilution of nuclear emission-line EWs by the stellar component along the line of sight, may readily prevent detection of accretion-powered nuclear activity in an early-type galaxy (Papaderos et al. 2013; Gomes et al. 2016), which is then classified as retired/passive and included in automated PSS studies. The same may obviously apply to old classical bulges, many of which show faint nebular emission despite hosting a super-massive black hole whose accretion-powered energy release may result in significant contamination of the optical spectrum by an AGN PL component. Our pilot study in CGP16 has first addressed the question of the detectability of an AGN PL embedded within an old instantaneously formed Lyman-photon leaking galaxy and explored the effect that this PL component may have on optical PSS modelling studies. This analysis has revealed that a PL contributing up to $\sim$26% ($\equiv x_{\rm AGN}^{\rm threshold}$) of the monochromatic luminosity at 4020 Å generally evades detection both from visual inspection and PSS modelling of a galaxy spectrum and introduces a substantial bias in the physical properties of the stellar component obtained through spectral synthesis. Here, we extend this study by modelling with the PSS code Starlight an extensive grid of synthetic stellar SEDs that trace the spectral evolution of galaxies for a wide range of SFHs over an age span between 1 Myr and 15 Gyr. A comparison of the physical properties obtained for the stellar component from PSS modelling with the input values constrained by the synthetic SEDs has allowed to examine in detail potential biases in PSS modelling for different fitting set-ups (Sect. 3).

Of special importance is the fact that, regardless of the spectral index $\alpha$ of the PL component integrated in the synthetic input SEDs, PSS fitting with purely stellar templates yields at $x_{\rm AGN}^{\rm threshold}$ an overestimation of stellar mass $M_{\star}$ by up to two orders of magnitude for instantaneously formed stellar populations of age $\la 10^{8}$ yr. Even though this bias is becoming gradually smaller with increasing age of the input SEDs, it still impacts $M_{\star}$ estimates by a factor of $\sim$2 over several Gyr of galactic evolution. For the same set of synthetic SEDs, PSS fits overestimate (underestimate) the light-weighted stellar age by up to $\sim$1 dex of stellar populations younger (older) than 10 Myr, whereas the mass-weighted stellar age is overestimated by up to $\sim$3 dex for young ages and a factor $\geq$2 for old ($\sim$5 Gyr) stellar populations. As for the mass-weighted stellar metallicity, our study indicates a systematic underestimation by a factor of $\sim$3 for young ages with this bias decreasing yet still present beyond an age of $\sim$1 Gyr. Such biases are also well documented for SFHs involving a prolonged star formation process (e.g. continuous or exponentially decreasing), suggesting a principal trend for purely stellar PSS fits to compensate for the presence of a PL spectral component through a severely overestimated contribution from old, low-metallicity stars.

Regarding the impact of an AGN PL on the stellar continuum, stellar population biases correlate on first order with $x_{\mathrm{AGN}}$ and on second with $\alpha$. These are mainly due to the combined effect of the dilution of stellar absorption features by the AGN continuum and a degeneracy between the AGN PL and young SSPs. The addition of a PL leads to a stellar age underestimation given the decrease of the Dn4000 Å break strength and simultaneously to a stellar metallicity overestimation due to the dilution of Balmer absorption features. Moreover, Balmer absorption features become weaker with age and their EWs are inversely correlated with Dn4000 for continuous SFHs (Kauffmann et al. 2003).

A question that naturally arises in view of these modelling biases pertains to the interpretation of the nature of evolved Lyman-photon-leaking galaxies hosting an AGN. Taken at face value, PSS fits might prompt the conclusion that these systems are extraordinarily massive and have assembled quasi-monolithically early on, which would be consistent with the much lower level of chemical enrichment as compared to local galaxies of equal $M_{\star}$. A possibly significant population of such old, ultra-massive galaxies erroneously included in galaxy statistics could then impact determinations of the cosmic SFH.

5 Conclusions

Synthetic spectra of active galaxies for solar metallicity and multiple SFHs were created with the ESS code Rebetiko and used to investigate how a simple AGN FC model impacts the estimation of stellar population properties with the state-of-the-art PSS code Starlight. The AGN continuum model is defined by its PL index $\alpha$ and AGN fractional contribution $x_{\mathrm{AGN}}$ to the $\lambda_{0}=4020$ Å monochromatic flux. The stellar spectral contribution is computed assuming solar metallicity SSPs from Bruzual & Charlot (2003) in both the creation of the synthetic spectra with ESS and the application of PSS. The accuracy of the estimated physical properties depends on the model age, $\alpha$, $x_{\mathrm{AGN}}$ and, more importantly, the approach adopted for dealing with the AGN continuum:

Excluding any PL component in the fit can lead to uncertainties on stellar mass by up to $\sim$3.5 dex, light- and mass-weighted mean stellar age up to $\sim$2 and $\sim$4 dex, respectively, and on both light- and mass-weighted mean stellar metallicity up to $\sim$0.7 dex. 2. 2.

Including a single or multiple PL components in the fit leads to uncertainties on stellar mass overestimation up to $\sim$1 dex, light- and -mass-weighted mean stellar age up to $\sim$1.5 and $\sim$3 dex, respectively, and light- and mass-weighted mean stellar metallicity up to $\sim$0.2 and $\sim$0.4 dex, respectively.

The uncertainties on the stellar population properties estimated with PSS are weakly dependent on the SFH, in particular when no AGN model is included in the base of elements. These results might lead to the misinterpretation that evolved Lyman-photon-leaking galaxies hosting an AGN are particularly massive, metal poor, and formed monolithically when the universe was still very young. Hence, these results show the importance of accounting for AGN spectral contributions when applying state-of-the-art PSS to active galaxies for a viable estimation of their physical properties, such as star formation and chemical evolution histories.

Acknowledgements.

This work was supported by Fundação para a Ciência e a Tecnologia (FCT) through national funds (UID/FIS/04434/2013) and by FEDER through COMPETE2020 (POCI-01-0145-FEDER-007672). We acknowledge support from the FCT through national funds (PTDC/FIS-AST/3214/2012) and by FEDER through COMPETE (FCOMP-01-0124-FEDER-029170). We thank the anonymous referee for valuable comments and suggestions. LSMC acknowledges funding by FCT through the grant CIAAUP-01/2016-BI in the context of the FCT project UID/FIS/04434/2013 & POCI-01-0145-FEDER-007672. JMG is supported by the fellowship SFRH/BPD/66958/2009 funded by FCT (Portugal) and POPH/FSE (EC) and by the fellowship CIAAUP-04/2016-BPD in the context of the FCT project UID/FIS/04434/2013 & POCI-01-0145-FEDER-007672. PP is supported by FCT through Investigador FCT contract IF/01220/2013/CP1191/CT0002 and by a Ciência 2008 contract, funded by FCT/MCTES (Portugal) and POPH/FSE (EC). We also acknowledge the exchange programme “Study of Emission-Line Galaxies with Integral-Field Spectroscopy” (SELGIFS, FP7-PEOPLE-2013-IRSES-612701), funded by the EU through the IRSES scheme.

Appendix A Additional resources

Figures 11 and 22 show the AGN fractional contribution $x_{\mathrm{AGN}}$ difference between PSS and ESS as a function of model age $t$ for models with instantaneous burst and continuous SFHs, respectively, with increasing $x_{\mathrm{AGN}}$ from left- to -right-hand side panels. Top and bottom row panels show results for the spl and mpl fitting set-ups (Section 2.2). On the one hand, the results for spl show that the accuracy of the estimated $x_{\mathrm{AGN}}$ increases with $\alpha$ and $x_{\mathrm{AGN}}$. Moreover, there is a underestimation bump between $10^{8}$ and $10^{9}$ yr that is attenuated with increasing $x_{\mathrm{AGN}}$. On the other hand, the results for mpl show a systematic $x_{\mathrm{AGN}}$ under- and overestimation for $t\lesssim 10^{8}$ and $t\gtrsim 10^{8}$ yr with increasing $x_{\mathrm{AGN}}$, respectevely. This underestimation increases with decreasing $\alpha$. Results for both spl and mpl show that the estimation of $x_{\mathrm{AGN}}$ is sensitive to the underlying stellar continuum shape.