SciELO - Scientific Electronic Library Online

 
vol.52 número1Halo occupation distributions of moderate X-ray AGNS formed through major and minor mergers in A Λ-CDM CosmologyUnveiling new associations in puppis índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados

Journal

Artigo

Indicadores

Links relacionados

  • Não possue artigos similaresSimilares em SciELO

Compartilhar


Revista mexicana de astronomía y astrofísica

versão impressa ISSN 0185-1101

Rev. mex. astron. astrofis vol.52 no.1 Ciudad de México Abr. 2016

 

Pipe3D, a pipeline to analyze integral field spectroscopy data: I. New fitting philosophy of FIT3D

S.F. Sánchez1 

E. Pérez2 

P. Sánchez-Blázquez3 

J.J. González1 

F.F. Rosález-Ortega4 

M. Cano-Díaz1 

C. López-Cobá1 

R. A. Marino5 

A. Gil de Paz5 

M. Mollá6 

A. R. López-Sánchez7 

Y. Ascasibar3 

J. Barrera-Ballesteros8 

1 Instituto de Astronomía, Universidad Nacional Autónoma de México, México. (sfsanchez@astro.unam.mx)

2Instituto de Astrofísica de Andalucía (CSIC), Granada, Spain.

3Departamento de Física Teórica, Universidad Autónoma de Madrid, Madrid, Spain.

4Instituto Nacional de Astrofísica, Óptica y Electrónica, Puebla, México.

5CEI Campus Moncloa, UCM-UPM, Departamento de Astrofísica y CC. de la Atmósfera, Facultad de CC. Físicas, Universidad Complutense de Madrid, Madrid, Spain.

6 Departamento de Investigación Básica, CIEMAT, Madrid, Spain.

7Australian Astronomical Observatory, Australia.

8Instituto de Astrofísica de Canárias (IAC), Tenerife, Spain.


Abstract:

We present an improved version of FIT3D, a fitting tool for the analysis of the spectroscopic properties of the stellar populations and the ionized gas derived from moderate resolution spectra of galaxies. This tool was developed to analyze integral field spectroscopy data and it is the basis of P 3D, a pipeline used in the analysis of CALIFA, MaNGA, and SAMI data. We describe the philosophy and each step of the fitting procedure. We present an extensive set of simulations in order to estimate the precision and accuracy of the derived parameters for the stellar populations and the ionized gas. We report on the results of those simulations. Finally, we compare the results of the analysis using FIT3D with those provided by other widely used packages, and we find that the parameters derived by FIT3D are fully compatible with those derived using these other tools.

Key Words: galaxies: ISM; techniques: spectroscopic

Resumen:

Presentamos una versión mejorada de FIT3D, una herramienta de ajuste para el análisis de las poblaciones estelares y el gas ionizado en espectros de galaxias de resolución intermedia. La misma se desarrolló para el análisis de datos de espectroscopía de campo integral y es la base de P 3D, un dataducto usado en el análisis de datos de los muestreos CALIFA, MaNGA y SAMI. Describimos la filosofía y los pasos seguidos en el ajuste, presentando un conjunto amplio de simulaciones con el fin de estimar la precisión de los parámetros derivados, mostrando el resultado de dichas simulaciones. Finalmente, comparamos el resultado del análisis con FIT3D y el obtenido mediante otros paquetes de uso frecuente, encontrando que los parámetros derivados son totalmente compatibles.

1. INTRODUCTION

The optical spectrum of a galaxy, or a part thereof, comprises information about the different components that emit or absorb light. Therefore, it can be considered as the net sum of different emitting sources, mostly stars and ionized gas, red- or blue-shifted and broadened by their particular kinematics, and attenuated by dust. All these components come together, co-added, to conform the observed spectrum, so they have to be decoupled in order to study their individual properties.

Fortunately, some of them present clear observational differences. While the stellar populations dominate the continuum emission, the ionized gas shines in a set of clearly defined emission lines at fixed rest-frame wavelengths defined by atomic physics. Several different tools have been developed to model the underlying stellar population, effectively decoupling it from the emission lines (e.g., Cappellari & Emsellem 2004a; Cid Fernandes et al. 2005; Ocvirk et al. 2006; Sarzi et al. 2006a; Sánchez et al. 2007a; Koleva et al. 2009; MacArthur et al. 2009; Walcher et al. 2011; Sánchez et al. 2011; Wilkinson et al. 2015). Most of these tools are based on the same principles. It is assumed that the stellar emission is the result of a single or a combination of different single stellar populations (SSP), or resulting from a particular star formation history (SFH). It is well known that the spectral energy distribution (SED) of simple stellar populations (chemically homogeneous and coeval stellar systems) depends on a set of first principles (e.g. initial mass function, star formation rate, stellar isochrones, metallicity, etc.), from which it is possible to generate the spectra of synthetic stellar populations. This technique, known as evolutionary synthesis modeling (e.g. Tinsley 1980), has been widely used to unveil the stellar population content of galaxies by reconciling the observed spectral energy distributions with those predicted by the theoretical framework.

Unfortunately, the variation of different physical quantities governing the evolution of stellar populations produces similar effects in the integrated light of those systems, leading to a situation in which the observational data are affected by hard to solve degeneracies, such as the ones involving age, metallicity, and extinction (e.g. Oconnell 1976; Aaronson et al. 1978; Worthey 1994; Gil de Paz & Madore 2002). Even so, the use of spectrophotometric calibrated spectra and the sampling of a wide wavelength range helps to break the degeneracy, and allows the derivation of reliable physical parameters by fitting the full spectral distribution with single stellar populations (Cardiel et al. 2003).

However, the simple assumption that a single stellar population describes well the SED of a galaxy is not valid even for early-type galaxies, less so for late-type ones. Most galaxies present complex star formation histories, with different episodes of activity, of variable intensity, time scales, and complex dust distributions. Therefore, a single stellar population does not reproduce well their stellar spectra. A different technique, known as full spectrum modeling, involving the linear combination of multiple stellar populations and the non-linear effects of dust attenuation, has been developed to reconstruct their stellar populations (e.g. Panter et al. 2003; Cid Fernandes et al. 2005; Tojeiro et al. 2007; Ocvirk et al. 2006; Sarzi et al. 2006b; Koleva et al. 2009; MacArthur et al. 2009).

In general, these reconstructions require a wide wavelength range to probe simultaneously the hot, young stars and the cool, old stars. They also require the best spectrophotometric calibration to disentangle the effects of age, metallicity, and dust attenuation. Although the different implementations of this technique have some differences, they are very similar in their basis. The information extracted from the multi stellar population modeling differs among implementations. In some cases, the luminosity (or mass) weighted ages and metallicities are derived, based on the linear combination of different models (e.g. Sarzi et al. 2006a). In other implementations the fraction of light (or mass) of different stellar populations is derived (e.g. Stoklasová et al. 2009; MacArthur et al. 2009). Finally, in other tools the information in the shape of the stellar continuum is considered unreliable and is removed prior to any further analysis (e.g. Ocvirk et al. 2006).

In addition to the actual composition of the stellar population, kinematic effects need to be taken into account. The stellar continuum should be redshifted by a certain velocity, broadened and smoothed to account for velocity dispersion, and attenuated due to a certain dust content. Some of the algorithms mentioned above perform this kinematic analysis prior to the decomposition of the underlying stellar population (e.g. Cid Fernandes et al. 2005), while in other algorithms it is a fundamental part of the computation (e.g. Sarzi et al. 2006 a).

Once the modelled stellar spectrum is determined, it is subtracted from the observed spectrum to provide a pure emission line spectrum that includes the information from the ionized gas. In some algorithms the stellar and ionized gas spectra are analyzed together, although in many cases a two-step procedure is used in the analysis (e.g. Sánchez-Blázquez et al. 2014a). Then, the main properties of the emission lines, including the intensity and kinematics, need to be derived. In general, given their larger signal to noise ratio, it is assumed that the information derived from emission lines is more accurate and stable than that derived for stellar populations (see Walcher et al. 2011, for a review of the state of the art). However, a proper subtraction of the stellar contribution is required prior to the measure of emission line ratios, such as the [ OIII ] /Hβ ratio. This subtraction allows to interpret the data in terms of physical processes.

In Sánchez et al. (2006b) we presented a new tool to perform all these analyses, mostly focused on the spectra obtained using Integral Field Spectroscopy (IFS), in the optical spectral range. This tool, named FIT3D, was first developed with the aim of analyzing the properties of the ionized gas through its emission lines, and has been used in several science publications and PhD theses. In the last few years FIT3D has increased its capabilities to derive a more reliable characterization of the underlying stellar population, and a good estimation of the errors of the derived parameters.

This article is the first of a series focused on the description of P 3D, a spectroscopic analysis pipeline developed to characterize the properties of the stellar populations and ionized gas in the spatially resolved data from optical IFU surveys, in particular CALIFA (Sánchez et al. 2012b), MaNGA (Bundy et al. 2015), and SAMI (Croom et al. 2012). P 3D uses FIT3D as its basic fitting package. In this article we present the new algorithms included in the distributed version of FIT3D. In § 2 we describe the new fitting philosophy and the details of the algorithms implemented. In § 3 we demonstrate, based on extensive simulations, the accuracy and precision of the parameters recovered for the stellar populations (§ 3.2), using a simulated single-burst stellar population (§ 3.2.1), or a fully simulated star formation and chemical enrichment history (§ 3.2.2). The reliability of the errors estimated for the derived parameters is described in § 3.2.3. Tests on the accuracy and precision of the parameters recovered for the emission lines and of the estimated errors are discussed in § 3.3 and 3.3.2, respectively. In § 3.4 we explore the precision in the parameters derived for both the stellar populations and the ionized gas, with more realistic simulations based on actual IFU data extracted from the CALIFA survey. A comparison with more broadly used fitting algorithms is presented in § 4. Finally, § 5 summarizes the main conclusions of this study.

2. NEW FITTING PHILOSOPHY OF FIT3D

The main goal of any analysis pipeline dealing with IFS galaxy data is to disentangle the main components comprising the spatially resolved spectra. As in the case of long-slit or single-fiber spectroscopic data, the two main components are: (i) the light coming from the stellar population, that dominates the continuum and absorption spectrum, and (ii) the emission of the ionized gas, that is observable as a set of emission lines at certain characteristic wavelengths. The two components are affected by dust attenuation (Av), produced by dust grains distributed across the galaxy (e.g. Tuffs et al. 2004). The disentanglement of these two main components is mandatory to properly study the properties of local and distant galaxies.

The rest-frame spectra of the two components are observed redshifted due the Doppler effect induced by the expansion of the Universe and the residual velocity of the galaxy, and further blue- or red-shifted by the particular kinematics of the target cell and component within the galaxy. In addition, the spectra are broadened by a characteristic velocity dispersion, accounting for the unresolved kinematics that comprises the perturbations with respect to the average velocity pattern, and the non-homogeneous motions (such as non-circular motions). In rotationally supported systems, like disks of spiral galaxies, the velocity dispersion is much lower than the asymptotic rotation velocity. On the other hand, in pressure supported systems, like bulge dominated early-type galaxies, the broadening from random orbits is much larger than the pure rotation velocity. Obviously, this picture is an oversimplification, since many other components may affect the kinematic pattern, like bars, outflows, interactions, and so on (e.g. Barrera-Ballesteros et al. 2014). Irrespectively of the particular kinematic pattern, its effect has to be taken into account by any analysis pipeline.

Former versions of FIT3D address this problem by (i) modeling the continuum with a linear combination of synthetic stellar populations (SSP), and (ii) modeling the emission lines with a set of single Gaussian functions. Details of the procedures adopted to fit both components are given in the sections below. In the case of the stellar populations, their weights or eigenvalues were derived without boundary constraints (Sánchez et al. 2006 b, 2011). However, if any of the weights derived was found to be negative, the corresponding SSP was removed from the template list, and the fitting procedure was repeated until only positive weights were derived. Under this procedure, only a small number of SSPs remain in the final iteration, limiting the extracted information to the average properties of the underlying stellar population, like the luminosity or mass weighted ages and metallicities. This procedure was also very sensitive to the effects of noise, and did not allow to estimate the errors of any of the derived quantities. Despite these limitations, this approach provides a reliable subtraction of the underlying stellar population.

In the current version of FIT3D9 we adopt a Monte-Carlo approach (MC), in which the previous procedure is iterated for a randomized version of the input spectrum taking into account the error vector. The individual weight vector is stored in each iteration. The final model of the stellar component is built using the mean weights derived along the MC sequence. In addition, this method provides realistic uncertainties of both the weights and the final model, based on the standard deviation with respect to the mean values. The details of how the MC is performed will be explained in the following sections.

Therefore, the analysis of the underlying stellar population, without considering the emission lines, comprises a non-linear and a linear problem. The non-linear problem includes the parameters that define the kinematic properties of the spectra and the dust attenuation, and the linear one involves the decomposition of the underlying stellar population in its components, either synthesized stellar populations (SSPs), or a library of stellar templates. The non-linear problem is addressed first and independently from the linear one.

The algorithms included in FIT3D were originally developed in , making an extensive use of the Perl Data Language (PDL 10). We have developed a version for which we have adopted most of the required algorithms included in PDL. The particularities of the version will be described elsewhere (Ibarra et al., in prep.).

2.1. Non-Linear parameters of the stellar populations

Unlike other tools (e.g., PXF Cappellari & Em-sellem 2004 b), FIT3D was not developed to provide a detailed kinematic analysis of the line-of-sight velocity distribution. However, it estimates the main properties of the kinematic structure, i.e., the stellar velocity (v sys) and velocity dispersion ( σ v ). These parameters are derived under the assumption that all stellar populations move at the same velocity, with a velocity dispersion following a Gaussian function. This is clearly an approximation, since it is well known that the velocity distribution of stellar populations deviates strongly from this simple approach (e.g. Rix & White 1992). However, most cases require a higher spectral resolution than the one presented here (R <2000) to make a clear distinction between this Gaussian approximation and more complex velocity distributions. The stellar kinematics is derived prior to the linear combination required to analyze the stellar population, based on an additional pseudo-random exploration within a pre-defined range of values. The pseudo-random exploration varies the parameter in question within a pre-defined range of values, but not following a fixed step. Instead it uses a random step, generated within +50% of the original step (by default, 1/30th of the parameter range). This loop is repeated several times (three by de fault). We consider that this non-regular exploration of the parameter space provides a better sampling within the given range of values.

For the velocity dispersion we adopted two different implementations: (i) measuring it in units of wavelength since for most of the spectra analyzed the dispersion is dominated by the instrumental one, and (ii) a different version that considers both a fixed instrumental dispersion along the spectral range in units of wavelength, and a physical velocity dispersion (in km/s). That second version requires a better knowledge of the instrumental dispersion and it is somewhat slower to run since it performs two different convolutions over the stellar templates.

In the first case, the stellar templates are convolved with a Gaussian function in wavelength space, and then shifted to the corresponding velocity. Although this is a crude approximation, for spectra dominated by instrumental resolution (e.g., PINGS and F-CALIFA data Rosales-Ortega et al. 2010; Mármol-Queraltó et al. 2011), there is no significant difference. In the second case, the stellar templates are first convolved in linear wavelength space with a fixed instrumental resolution (as in the first case), and then they are convolved in log(λ) space with a Gaussian that accounts for both the velocity and the velocity dispersion.

The fitting sequence is as follows: First, the velocity dispersion is fixed to an input guess value, and the velocity is changed incrementally from a minimum to a maximum value, using a semi-random step around a certain fixed velocity step. For each velocity a best linear combination of SSPs that reproduces the continuum is derived, based on the reduced X 2. Then a range of velocity dispersions is explored, following a similar procedure, now fixing the velocity to the best result found in the previous step. The procedure is repeated twice, with the second iteration focused on the best parameters found in the first one, and reducing both the range and the size of the step to 10% of the initial values. Figure 1 illustrates this process. It shows a detail of the central spectrum (5" aperture), extracted from the V500 setup data cube of NGC 2916 from the CALIFA survey, together with the best SSP templates derived for each of the scanned velocities and velocity dispersions within the ranges considered.

Fig. 1 Left Panel: Detail of the central spectrum (5" aperture, black solid-line) of NGC 2916 extracted from the V500 data cube of the CALIFA survey, together with the best fitted SSP template for the different scanned velocities (color-coded solid lines). The blueish solid lines correspond to SSP models blue-shifted with respect to the velocity, and the reddish solid lines correspond to red-shifted ones. Right Panel: Same spectrum shown in the left panel (black solid-line), together with the best fitted SSP template for the different scanned velocity dispersions (color-coded solid lines). The blueish solid lines correspond to SSP models with smaller velocity dispersion, while the reddish solid lines correspond to those with larger velocity dispersion. In both panels the residuals from the fits are plotted using the same color scheme, together with the residual of the best fitting model (grey solid-line). The two panels show two consecutive steps in the fitting procedure. The color figure can be viewed online. 

Once the velocity and velocity dispersion are derived, then the range of dust attenuations that better reproduces the data is explored. Again, a range of values is explored between a minimum and a maximum dust attenuation using a pseudo-random exploration with a variable step between consecutive values. Following the same procedure described before, the best linear combination of SSPs is derived for each particular value of the dust attenuation. Figure 2 illustrates this procedure. It shows the full wavelength range of the central spectrum (5״ aperture), extracted from the V500 setup data cube of NGC 2916 from the CALIFA survey, together with the best SSP templates derived for each of the scanned values of dust attenuation within a range between 0 and 1.6 mag.

Fig. 2 Central spectrum (5" aperture, black solid line) of NGC 2916 extracted from the V500 data cube of the CALIFA survey, together with the best fitted SSP template for the different scanned dust attenuations (Ay) as color-coded solid lines ranging from 0 to 1.6 mag. The blueish solid lines correspond to SSP models with low dust attenuation, while the reddish solid lines correspond to those with larger attenuation. The residuals from the fits are plotted using the same color scheme, together with the residual of the best fitting model (grey solid line). The color figure can be viewed online. 

Obviously, as described, this procedure would be computationally very intensive, and impractical in its current form. Thus we have implemented some modifications to speed up the process: (i) The MC iterations on the linear combination of SSPs are not applied in the steps required to derive the velocity, velocity dispersion, and dust attenuation; thus, a single linear combination is derived for each step. The fitting is performed using a standard matrix inversion method to do a least squares minimization weighted by the errors, as implemented in the PDL routine 1 11; (ii) a library with a smaller number of single SSPs was adopted in these initial steps (although it is not mandatory or restricted to do so); (iii) a restricted wavelength range could be considered for the derivation of the kinematic components, taking into account those regions with stronger stellar absorption features (e.g., 3700-4700 Ǻ in the case of CALIFA and MaNGA data). This approach is new, since in general it was considered that using just the regions with stronger stellar absorptions would create a bias towards the velocity of the younger stars. Indeed, the H-K Ca lines can be broadened due to the rotation of hot stars (e.g. Gerhard 2003; Walcher et al. 2005). However, for the range of velocity dispersion expected in galaxies this seems not to have a strong effect, based on the tests that we have performed; (iv) a well defined range for the velocity and velocity dispersion is derived by analyzing the integrated spectrum and the central/peak spectrum for each galaxy. In this initial exploration a wider range of parameters is considered. Extensive tests have shown that these approaches do not affect the results within the expected uncertainties, as we will describe later.

In all the three pseudo-random explorations, for each parameter we consider an input guess, a range of values, and an input step for which the parameter space is explored. The input guess is used as a fixed tracer of one of the parameters when the other two are explored. Therefore, the accuracy and precision of the derivation is limited to the input parameters injected. In order to perform the pseudo-random exploration, instead of using a constant step between the minimum and maximum, the step is variable within 80% of the input step. This ensures that the exploration is not uniform, and, due to the double iteration, it allows a better estimation of the parameter. For each iteration the best reduced X 2 derived from the linear analysis of the underlying stellar population is stored. Then, the minimum X 2 together with the two adjacent values are used to derive the parameter value using a parabolic minimum derivation (Sánchez 2006a), together with a rough estimation of the error in the parameter as the range in which the reduced X 2 changes up to X 2 + 1.

2.2. Linear parameters of the stellar populations

Once the non-linear parameters are derived, they are fixed in the derivation of the properties of the underlying stellar population, reducing the problem to the solution of a linear combination of templates (SSPs).

The implementation of the multi stellar-population fitting technique used in this article was already described in Sánchez et al. (2006 b, 2011). The basic steps of the fitting algorithm, spectrum by spectrum, are the following:

  • 1. Read the input and noise spectra and determine the areas to be masked. Let us define 0 i as the observed flux at a certain wavelength λ i, σi as the noise level at the same wavelength, and N as the number of elements of the masked spectrum.

  • 2. Create a Monte-Carlo realization of the data, adopting a Gaussian random distribution for the noise. Let us define G i as the MC-realization of the observed flux, defined as

  • Gik=Oi+ Riσi

  • where R i is a random value following a (-1,1) Gaussian distribution, and k is the running index of the current MC realization.

  • 3. Read the set of SSP template spectra, shift them to the velocity of the spectrum, convolve them with the velocity dispersion, and re-sample to the wavelength solution. Let us define F ji as the flux of the i th wavelength of the j th template (once shifted, convolved, and re-sampled), where M0 is the total number of templates used.

  • 4. Apply the derived dust attenuation to the templates. The attenuation law of Cardelli et al. (1989) was adopted, with a ratio of total to selective attenuation of Rv=3. 1. Let us define F AV ji as the flux of the i th wavelength of the j th template, after applying the dust attenuation corresponding to a certain attenuation of Av magnitudes.

  • 5. Perform a linear least-square fit of the input spectrum with the set of redshifted, convolved, and dust attenuated SSP templates. The fitting is performed using a standard matrix inversion method, weighted by the errors, as described before, based on the PDL routine 1 . Then a modified X 2 is adopted as the merit function to be minimized, with the form:

  • X2= 1N-M i=1NDi2,

  • where

  • Di= wi GiK- j=1Maj FjiAV,

  • In the above expression, w i is the weight of the i th pixel, defined as

  • wi= 1σi2 ,

  • a j is the coefficient of the j th template in the final modeled spectrum, and M is the number of templates considered in the fitting procedure. The actual number of templates included in the SSP-library depends on the science goals and the quality of the explored data. In our case for a typical N ≈ 2000, the number of templates if of the order of M ≈ 100. There is no general recommendation of how large should be Mto derive optimal results. It is therefore required to perform simulations to test it (e.g. Cid Fernandes et al. 2014). However, our experience indicates that it should not be larger than 2-3 times the typical signal-to-noise of the spectra, as a practical rule.

  • 6. Determine for which templates within the library a negative coefficient is derived in the linear combination (i.e. a j < 0). These templates will be excluded from the next iteration of the fitting procedure, that will be resumed in step (5). At each iteration, M is decreased by the amount of templates excluded. This loop ends once all the coefficients are positive. It is important to note here that excluding those templates within the library that produce negative values for the coefficients does not produce in general an increase of the reduced X 2 of the fitting. Therefore, to remove them does not decreases the goodness of the fitting in a systematic way. This procedure is adopted both in the linear and non-linear exploration of the parameters.

  • 7. The modeled spectrum flux at λ t for the kth MC realization is given by

  • Sik= j=1Majk FjiAV ,

  • where now M considers only the templates with positive coefficients a k j .

  • 8. The procedure is iterated for a fixed number of MC realizations (K) from step (ii). The final modeled spectrum, with its corresponding uncertainties at each λ i is given by the mean and standard deviation of the individual models for each realization of the MC iteration, at the given wavelength:

  • Si=mean (Si1..k)

  • σSi=stddev ( Si1..k)

  • The final reduced X 2 is then recomputed using this final average modeled spectrum. It is found that in general this X 2 is not significantly different from the individual ones for each of the MC realizations (the same is valid for the standard deviation, in a consistent way). This confirms the assumption that each individual fitting reproduces the analyzed spectrum equally well, and the final distribution is a more realistic description of the properties of the stellar population.

  • The final coefficients of the SSP decomposition and the corresponding uncertainties are derived using a similar procedure:

  • aj=mean (aj1..k)

  • σaj=stddev ( aj1..k)

  • where σ par is the uncertainty associated with the corresponding parameter.

For practical reasons, the algorithm provides the fraction of the total flux contributed by each particular SSP within the library at a certain wavelength and therefore a normalization wavelength is recommended: Let us define c j as these relative flux fractions. Along this article we have adopted the normalization at 5500 Ǻ; however, this is not mandatory..

In the derivation of the non-linear parameters a reduced version of this algorithm was adopted, excluding steps (2), (6), and (8). Thus, only a single linear least-square fit is applied, without cleaning the negative coefficients (that are normally absent when the template adopted is limited to a few stellar populations), and not iterating over a sequence of MC realizations.

This algorithm shares a basis with other procedures described in the literature (e.g., Cappellari & Emsellem 2004b; Cid Fernandes et al. 2005; Ocvirk et al. 2006; Sarzi et al. 2006a; Walcher et al. 2006; Sánchez et al. 2007 a; Koleva et al. 2009; MacArthur et al. 2009). If required, the low frequencies in the spectrum can be fitted by adding a polynomial function to the fitted templates, or by multiplying them by that function. There are different origins for these low frequencies in the spectrum, mostly related to defects in the spectrophotometric calibration of the data or the stellar libraries adopted for the SSP templates (e.g. Cid Fernandes et al. 2013) They may also be related to the way that the synthetic populations are generated, like differences/errors in the adopted IMF or the stellar evolution models being different for different stellar templates (e.g. Cid Fernandes et al. 2014). In the particular implementation shown along this article we have performed this correction only to subtract the stellar component to analyze the emission lines. However, for the full analysis of the stellar populations it is not implemented. This is evident in the residuals shown in Figure 3.

Fig. 3 Left Panel: Detail of the central spectrum (5" aperture, black solid line) of NGC 2916 extracted from the V500 data cube of the CALIFA survey, together with the different MC-realizations of this spectrum created along the fitting procedure (color-coded points). The real noise has been increased by a factor of 5 to visualize the different realizations more clearly. The residuals from the best model fitted in each MC-realization are included using the same color code as the model, together with the average of the different residuals (solid black line). Right Panel: Same spectrum shown in the left panel (black solid line) together with the best fitted SSP template derived using the MC linear regression technique (red solid line). The residual from the best fitting model is shown as a solid black line. The color figure can be viewed online. 

By construction, the fitting algorithm requires that certain regions in the spectrum be masked: (i) strong and variable night sky emission line residuals, (ii) regions affected by instrumental signatures, like defects in the CCD (e.g. dead columns) whose effect was not completely removed during the data reduction process, (iii) regions affected by telluric absorptions, not completely corrected during the flux calibration process, and (iv) regions containing strong emission lines from the ionized gas.

In the case of FIT3D a final iteration is implemented that makes this latter masking not critical. It is possible to combine the analysis of the stellar population and emission lines fitting in a single iterative process in which, for each iteration, the best model of the emission lines is removed prior to the analysis of the stellar population. We will describe later this possibility.

2.3. Average properties of the stellar populations

Once the best model for the underlying stellar population is derived based on a linear combination of SSPs, it is possible to derive the average quantities that characterize this population. The most common parameters derived are:

  • • The luminosity-weighted log age (Age L) and log metallicity (Z L) of the underlying stellar population,

  • logAgeL= j=lMcjlogAgej,

  • and

  • logZL= j=lMcjlogZj,

  • where Age j and Z j are the corresponding age and metallicity of the jth SSP template, and c, is the weight in light of each SSP. In this context we define Z as the mass fraction of metals, i.e., elements heavier than hydrogen and helium.

  • Hereafter, for the sake of simplicity, we will refer to the previous values as the 'luminosity-weighted age' (instead of the more cumbersome luminosity-weighted logarithmic age). The reader should take into account that a different definition would be:

  • AgeL= j=lMcjAgej,

  • and

  • ZL= j=lMcjZj

  • Note that the logarithm of the latter ones does not correspond to the former ones. Indeed the former ones are the geometrical average, weighted by the fraction of light corresponding to each single stellar population, while the latter ones correspond to the arithmetic average. Since the time sampling of the stellar templates is logarithmic, we consider the former one a better representation of the average ages and metallicities in galaxies, although it is possible to derive any of them using the FIT3D output.

  • • The mass-weighted log age (Agen M) and log metallicity (Z M ) of the underlying stellar population,

  • logAgeM= j=lMcjMLjlogAgej,

  • and

  • logZ= j=lMcjMLjlogZj,

  • where Age j and Z j are the same parameters described before, and ML j is the mass-to-light ratio of the jth SSP template (i.e., ML ).

  • Hereafter we will refer to them as the mass-weighted age and metallicity. Similar caveats should be applied to this definition as the one expressed regarding the luminosity-weighted ones.

  • •The average mass-to-light ratio,

  • ML= j=lMcj MLj

  • that in general is not the same as the mass-to-light ratio of the integrated spectra, that would be given by:

  • ML= j=lMMjj=lMLj

  • where

  • Mj = c j ML jL

  • and L is the total luminosity of the spectrum.

  • • We define the star formation history (SFH) as the evolution of the star-formation rate along the time. To derive it, we define the mass of stars of a given age:

  • Mage= j=lMMj, if agej=age.

  • Now, we estimate the cumulative amount of mass up to a time τ age where τ age = τUniverse - age and

  • τUniverse is the age of the Universe at redshift zero. To do so, we integrate the amount of stellar mass since the beginning of the cosmological time until time age):

  • Mage,c= j=lMMj, if agej>age.

  • From these cumulative masses it is possible to derive the star-formation rate (SFR) at any particular time.

  • SFRτage= Mageage ;

  • the distribution of these SFR τ age conforms the SFH as defined before.

  • At any time the mass remaining in stars is obtained by:

  • Mage,cor= j=lMfcor,jMj, if agej>age,

  • where f corj is the correction factor that takes into account the mass-loss and the mass locked into remnants, for each stellar population at the given age and metallicity (e.g. Courteau et al. 2014).

  • For each spectrum, both the SFH and the SFR are not single parameters, but arrays of a length equal to the number of ages included in the SSP template library.

The luminosity-weighted ages and metallicities should be considered as the first momentum of the distribution of weights, or, in other words, the first momentum of the star formation and chemical enrichment histories (in a logarithm scale as defined here). They differ considerably from the effective ages and metallicities, since in general they do not match the corresponding values if the population were composed of a single SSP (Serra & Trager 2007).

The luminosity-weighted ages and metallicities highlight the contribution of the young stellar populations (with a strong color effect), which are much more luminous for a certain mass than their older relatives. However, the mass-weighted ages and metallicities are less sensitive to the young stellar populations, and therefore trace better the bulk of the star formation history, rather than more recent star-forming events.

2.4. Characterization of the emission lines

One of the goals of this fitting procedure is to provide an accurate characterization of the underlying stellar population, to subtract this contribution from the original spectrum, deriving an emission-line only spectrum. This gas emission spectrum is derived by

Ci = Oi - Si.

To measure the intensity of each emission line detected, each emission line in the clean spectrum (C¡) is fitted with a single Gaussian function, plus a low order polynomial. Following the new philosophy described below, the emission lines are fitted splitting the non-linear and the linear components. As in the case of the stellar population, the Gaussian function comprises two non linear parameters, the velocity and the velocity dispersion, and a linear one, the intensity. On the other hand, the polynomial function included to fit the continuum only comprises a combination of linear components. The procedure is very similar to the one outlined before. First, we perform an pseudo-random exploration of the nonlinear parameters within a range of values, starting with the velocity, and then the velocity dispersion. An initial guess of the velocity dispersion is used as a fix parameter in the exploration of the velocity, and later the velocity derived is fixed in the exploration of the velocity dispersion. In each step of the pseudo-random exploration a least-squares linear regression is performed to derive the linear parameters (i.e., the intensity of each emission line and the coefficients of the polynomial function) using the same standard matrix inversion method, weighted by the errors, described before, which is based on the PDL routine 1 . The reduced X 2 is stored in each iteration as a proxy of the the goodness of the fitting.. The final best combination of non-linear and linear parameters is derived on the basis of the minimization of the reduced X 2 .

We estimate the errors of the optimized parameters by performing a MC simulation, where the original emission line only spectrum is perturbed by a noise spectrum that includes both the original estimated error and the uncertainties in the best fit SSP model. Therefore, the intensity fitted at each wavelength λ i in each MC realization is

where i is the spectral pixel, R i is a random value following a (-1,1) Gaussian distribution, k is the running index of the current MC realization, and

where σ i is the original error (the one used in the analysis of the SSP), and σ si. is the standard deviation around the best SSP model, S i (i.e., the uncertainty in this model), as defined in point 8 of § 2.2. Using this approach, we propagate the uncertainties of the subtraction of the underlying stellar population to the parameters derived for the emission lines.

Instead of fitting the entire wavelength range at once, for each spectrum we extract shorter wavelength ranges that sample one or a few emission lines. This allows us to fit the residual continuum with the simplest polynomial function (in general, of 1st or 2nd order), and to simplify the fitting procedure. When more than one emission line is fitted simultaneously, their velocities and dispersions may be coupled. Thus, they are forced to be equal in order to decrease the number of free parameters and increase the accuracy of the de-blending process (if required). This is particularly useful when the velocity dispersion is dominated by the spectral resolution, and for the recovery of the intensity of weak emission lines, with poorly defined kinematic properties. Figure 4 illustrates this procedure. It shows two spectral ranges centered on Hβ and Hα, of the central spectrum of NGC 2916 shown in Figure 1-3, once we subtracted the best fitted SSP model shown in Figure 3, together with the best fitted Gaussians to Hβ and [OIII]λ4959,5007, and to Hα and [NII]λ6548,6583, respectively. The different MC realizations adopted during the fitting procedure are shown to illustrate how the procedure evaluates the parameters and errors at the same time.

Fig. 4 Left Panel: Detail of central spectrum (5" aperture, black solid line) of NGC 2916 extracted from the V500 datacube of the CALIFA survey around the wavelength range of Hβ and [O III ] λ4959,5007, together with the different MC-realizations of this spectrum created along the fitting procedure (color-coded points). The best fitting model to the emission lines is shown as a green dashed line. The residuals from the different realizations are shown following the same color scheme (offset by -2.5 flux units for the sake of clarity), together with the residual of the best fitting model. Right Panel: Same spectrum shown in the left panel, but now around the wavelength range of Ha and [Nil ] λ6548,6583 (black solid line), together with the different MC realizations of this spectrum created along the fitting procedure (color-coded points). The best fitting model to the emission lines is shown as a green dashed line. In both panels the level of noise has been increased by a factor of 5 to visualize the different realizations more clearly. The color figure can be viewed online. 

A modeled emission line spectrum is created, based on the results of the last fit, using only the combination of Gaussian functions. This spectrum is given by

where L is the number of emission lines in the model, B k is the integrated flux of the K th emission line, and Gauss ik is the corresponding normalized kk) Gaussian evaluated at the ith pixel.

Finally, this modeled emission line spectrum is subtracted from the original spectrum (G i), deriving an emission line free spectrum given by

GFi = Gi - Ei.

This spectrum may be used to model again the stellar population, by applying the procedure described before, but without masking the emission line spectral regions, as outlined in § 2.2. Actually, our experience is that adopting a simple SSP template library for the initial iteration, and a more detailed one for the second iteration produces very reliable results for both the emission lines and the stellar populations. In general, a third iteration does not produce any significant change in the derived properties, although the implemented algorithm allows the user to select the desired number of iterations. It is important to remark that the final reduced X 2 adopted includes both models and the propagation of the uncertainties.

The full data-flow described in the previous sections is summarized in Figure 5. There, different steps of the fitting sequence are described, illustrating the products derived in each of them and the interdependence between those steps. The algorithm implemented in FIT3D allows the user to fit the stellar population without fitting any emission line, and at the same time it includes routines to fit the emission lines irrespectively of the underlying stellar population, in case this is required.

Fig. 5. Chart-flow summarizing the fitting sequence of FIT3D as described in the text. The flow indicates the main blocks of the fitting process, showing which operations are performed, and the main data-products derived. The MC-loop is indicated with a bracket, and the final iteration over the emission-line fitting with black arrows. In each box a label indicates either the section in the text where the procedure is described or the point in § 2.2 in which is explained. 

3. ACCURACY OF THE DERIVED PARAMETERS

The procedure outlined above to fit the stellar population and the emission lines implements a somewhat new concept. It does not rely on minimization algorithms frequently used to fit non-linear functions, like the LM algorithm implemented in other tools (e.g., Gandalf, the previous version of FIT3D), but on a sequential combination of pseudo-random exploration of the range of non-linear parameters, and a least-square inversion of the linear ones. The MC realizations and pseudo-random explorations share some ideas with other tools that explore the properties of the stellar populations, like S (Cid Fernandes et al. 2005). However, their minimization scheme is based on a sequence of Markov chains, and not on a sequence of linear regressions.

Therefore, we need to demonstrate that the new procedure yields reliable results. To do so, we perform a set of simulations to verify that (i) the adopted procedure creates accurate models of the data, (ii) the main parameters of the stellar populations and emission lines correspond at least to the expected ones with some well controlled inputs, and (iii) the procedure provides a reliable estimation of the uncertainties of the derived parameters.

3.1. SSP template library

As already noted by different authors (e.g. MacArthur et al. 2004; Cid Fernandes et al. 2014), this kind of analysis is always limited by the template library adopted, which comprises a discrete sampling of the SSP ages and metallicities. It is desirable that the stellar library be as complete as possible, and non-redundant. However, this would require an exact match between the models and the data, which is generally not possible to achieve, in particular if the stellar population includes more than one SSP.

The results should be as independent as possible of the currently adopted SSP template library. Therefore, we have tested the fitting routines using different ones, following Cid Fernandes et al. (2014):

  • gsd156 template library: described in detail by Cid Fernandes et al. (2013), it includes 156 templates that cover 39 stellar ages (1 Myr to 13 Gyr), and 4 metallicities (Z/Z ʘ = 0.2, 0.4, 1, and 1.5). These templates were extracted from a combination of the synthetic stellar spectra from the GRANADA library (Martins et al. 2005) and the SSP library provided by the MILES project (Sánchez-Blázquez et al. 2006; Vazdekis et al. 2010; Falcón-Barroso et al. 2011). This library has been extensively used within the CALIFA collaboration in different studies (e.g. Pérez et al. 2013; Cid Fernandes et al. 2013; González Delgado et al. 2014). Therefore, it is very interesting to know how accurate are our results using this particular library. Its spectral resolution has been fixed to the spectral resolution of the CALIFA V500 setup data (FWHM≈6 Ǻ).

  • miles72 template library: It includes 72 templates that cover 12 stellar ages (Myr to 12.6 Gyr), and 6 metallicities (Z/Z ʘ = 0.015, 0.045, 0.2, 0.4, 1, and 1.5). These templates were extracted from the SSP library provided by the MILES project (Sánchez-Blázquez et al. 2006; Vazdekis et al. 2010; Falcón-Barroso et al. 2011). They do not cover stellar populations as young as the previous ones, but they cover a wider range of metallicities. This template has a higher spectral resolution than the previous one, with FWHM≈2.5Ǻ (Falcón-Barroso et al. 2011).

  • bc 138 template library: It includes 138 templates that cover 23 stellar ages (1 Myr to 13 Gyr), and 6 metallicities (Z/Z ʘ = 0.005, 0.02, 0.2, 0.4, 1, and 2.5). The SSP models were created using the GISSEL code (Bruzual & Chariot 2003), assuming a Chabrier IMF (Chabrier 2003). They cover a range of metallicities similar to the previous one, but with a wider range of ages. However, they do not benefit from the spectrophotometric accuracy of the MILES stellar templates (Sánchez-Blazquez et al. 2006). The spectral resolution is slightly higher than that of gsd156, with a FWHM≈5.5Ǻ.

  • mar136 template library: It includes 136 templates that cover 34 stellar ages (0.1 Myr to 15 Gyr), and 4 metallicities (Z/Z ʘ = 0.05, 0.4, 1, and 2.23). The SSP models were created using the code by Maras-ton (2005), assuming a Salpeter IMF (Salpeter 1955). They cover a range of metallicities slightly smaller than the previous one, but with a similar range of ages. It is the template with the lowest spectral resolution, FWHM≈40Ǻ. This template was included to test the ability to recover the stellar population parameters when the spectral resolution is not very accurate.

Figure 6 shows the ranges of shapes covered by the different SSPs and how they compare with each other. For each SSP library we present the individual spectra within the wavelength range corresponding to the optical range between 3700-7000Ǻ, normalized at 5500Ǻ. We note that the sampling of all of these libraries differs in the distributions of age and metallicity. This is intrinsic to the templates, and imposes a bias in the sampling of the properties that is not intrinsic to FIT3D. To our knowledge, none of these libraries include the nebular continuum, which may produce a significant effect for very young stellar populations.

Fig. 6 Stellar population templates included in the gsd156 (first panel), miles72 (second panel), bc138 (third panel) and mar 136 (fourth panel) libraries. Each spectrum within each library is represented with a different color, with the blueish colors representing the younger stellar populations, increasing in age towards more reddish/orange/grey colors. The color figure can be viewed online. 

There are several caveats when applying model SSPs to the integrated light of an entire galaxy (or a portion of it), which have been clearly identified by previous authors (e.g. MacArthur et al. 2004; Cid Fernandes et al. 2014). The most important one is the assumption that the parameter space covered by the empirical library represents well that of the real data. However, in general, libraries are based on stars in the solar neighborhood (e.g. Sánchez-Blázquez et al. 2006) or on synthetic models (e.g. González Delgado et al. 2005; Maraston 2005), and therefore they may not represent well the stellar populations of other galaxies (or even of other regions of our Galaxy).

Most of these problems are not particularly important in the context of our tests, since our primary goal is to determine the accuracy and precision of the parameters recovered for a set of simulated data; however, they do affect the interpretation of the results.

In addition, it is important to note that the treatment of the dust attenuation may affect the resulting parameters (i.e. the luminosity-weighted age and metallicity of the stellar population). In this particular implementation of the analysis we adopted the Cardelli et al. (1989) attenuation law, which may not be the optimal solution to study the dust attenuation in star forming galaxies (e.g. Calzetti 2001). Other authors (MacArthur et al. 2009) assumed a completely different attenuation law, based on the two-component dust model of Chariot & Fall (2000), which is particularly developed to model the dust attenuation in star forming galaxies. The differences between the existing extinction laws in the optical wavelength range (3600-10000 Ǻ) are too small to affect significantly the results12. Indeed, no major differences are found when adopting the Cardelli et al. (1989), Calzetti (2001), or even a simple λ-l 3 extinction law. FIT3D adopted a fixed specific dust attenuation of Rv =3.1 (i.e., Milky-Way like), and the only fitted parameter is the dust attenuation at the V-band (Av), in magnitudes.

3.2. Accuracy of the properties of the stellar population

As mentioned before, the main parameters derived from the analysis of the stellar component are (i) the coefficients of the SSP templates in which the stellar continuum is decoupled, (ii) the luminosity- and mass-weighted ages and metallicities, and dust attenuation, (iii) the stel lar mass-to-light and mass included in the spectrum, and (iv) the star formation history. In order to assess the accuracy and precision of these parameters we have performed a set of simulations. The simulations should be performed trying to match as closely as possible some real data, especially the noise pattern. Since the primary goal of the currently developed package is the analysis of IFU data extracted from the CALIFA data, we will limit our tests to the wavelength range and spectral resolution of the CALIFA-V500 setup (Sánchez et al. 2012 a): a wavelength range between 3700-7500 Ǻ, and an instrumental resolution of inst =2.6 Ǻ. However, the simulations can be easily adapted to any other spectroscopic configuration.

The noise pattern is composed, in general, of both white noise (corresponding to the photon-noise of the source and the background, and electronic noise from the detector), and non-white noise (corresponding to defects/inaccuracies in the sky-subtraction, uncorrected defects in the CCD, errors in the spectrophotometric calibration, etc.). These noise patterns are different spectrum-to-spectrum, and wavelength-to-wavelength and are clearly difficult to simulate on a simple analytical basis. However, for a general test like the one performed here, we will limit our simulations to the use of white noise. In forthcoming articles we will explore the accuracy and precision of the parameters recovered using a more ad hoc set of simulations, following Cid Fernandes etal. (2014).

We perform two different sets of simulations. In the first one, a SFH dominated by a single burst of star formation is assumed, and therefore the stellar populations are dominated by stars of a certain age. The burst follows a Gaussian shape, with a width 0.5 dex in age (i.e., the length of the burst is proportional to the age selected). Under this assumption the burst lasts for a longer time at earlier times, and it is much shorter at later times. For the metallicity, we considered a similar probability distribution. Therefore, the simulated stellar population is dominated by stars of a certain metallicity with a Gaussian distribution of-width of 0.5 dex at each metallicity. This does not attempt to model any real physical enrichment mechanism in particular. Each simulated stellar spectrum is created by randomly selecting (within the range of the values considered in the SSP template) a particular age and metallicity as the center of the 2D Gaussian distribution.

In general, this simulated SFH may not be realistic. In most SFHs analyzed within the CALIFA survey most of the stars are formed following an exponential law with more than 80% of the mass formed in the first few Gyr (Perez et al. 2013). If we assumed this more realistic SFH we would not be able to test the ability of the fitting tool to decouple different stellar populations, and we would just test the very particular case that is more frequently observed in the Local Universe.

In the second set of simulations we try to create a more realistic SFH, in which the contribution of a stellar population within the template is proportional to a power of its age in Gyr (a j ∞ Agej 1/2). This is basically an exponential star formation rate along the evolution of the stellar population. The metallicity distribution follows a Gaussian probability function similar to the one described before. However, in this particular case the input metallicity is forced to follow a linear anti-correlation with age; i.e., older stellar populations are more metal poor than younger ones, following the recent results by González Delgado et al. (2014).

3.2.1. Results from the single burst simulation

For the first single burst SFH we created four dif ferent sets of simulations. In each one, 1000 simu lated spectra were created with the same normalized flux at 5000 Ǻ, corresponding to 100, 50, 10 and 5 10-16 erg s-1 cm-2 Ǻ-1, respectively. The velocity was randomly selected within the redshift range of the CALIFA footprint (0.005 < z < 0.03), and the velocity dispersion was selected randomly from σ = 2.5 - 7.5 Ǻ v 100 - 400 km/s). Then, the simulated spectrum was reddened by a dust attenuation randomly selected between Av = 0 - 1.1 mag. Finally, a purely Poissonian noise was added, to a level of 1 10-16 erg s-1 cm-2 Ǻ-1. As indicated before, this is not a totally realistic noise pattern. However, to perform a more detailed simulation would require to adopt ad hoc noise patterns, attached to a particular data set (e.g. Sánchez et al. 2011).

The fitting procedure was then applied to the simulated spectra using the same steps adopted to analyze the observational data, deriving, for each one, the different data products described in § 2.3. Figure 7 illustrates the result of this exercise for a particular spectrum with a signal-to-noise ratio of ≈ 100. It shows the input coefficients for a particular simulated spectrum created using the miles 72 template library (for both the simulation and the analysis), in comparison with the output distribution. This kind of distribution comprised the star-formation history and metal enrichment pattern, which are projections of this distribution. It is clear that the luminosity-weighted values derived from both distributions are very similar (log(age/yr)in = 9.76 vs. log(age/yr)out = 9.78 and [ Z/H ] in = -1.68 vs. [ Z/H ]out = -1.56), despite the fact that there are clear fluctuations between the input and recovered weighted/contributions of each particular stellar component. The largest discrepancies were found for the lower metallicity bin and the log(age/yr)=9.50 component, both of which differ by about 10%. Furthermore, the distribution of coefficients in the age-metallicity grid was remarkably similar. Previous analyses were performed by Ocvirk et al. (2006), Koleva et al. (2009), and Sánchez et al. (2011), although they covered a more reduced parameter space, and in general the kinematics and/or the dust attenuation were fixed in their simulations.

Fig. 7 Left panel: Distribution of input coefficients along the age-metallicity grid for a single simulated spectrum corresponding to the first set of simulations, i.e., those resembling a single burst of star-formation, with a broad extension in time. This particular example corresponds to a simulation based on the miles72 template library. For each age-metallicity pair the weight is shown as a solid square, with the size and color being proportional to the weight of the corresponding SSPage,met. For those SSPs with stronger contribution (above a 5%), the actual fraction is indicated on top of each square. The horizontal and vertical histograms represent the co-added coefficients for a particular age and metallicity, respectively. The green solid-star indicates the location of the luminosity-weighted age and metallicity corresponding to the distribution of stellar populations: log(age/yr)= 9.76, log(Z/0.02)= -1.68. Right panel :Similar distribution corresponding to the coefficients recovered for the given input simulation of a spectrum with a signal-to-noise ratio of ≈100, once the fitting procedure is applied. The green solid-star indicates the location of the luminosity-weighted age and metallicity corresponding to the recovered distribution of stellar populations: log(age/yr)= 9.78, log(Z/0.02)= -1.56. The color figure can be viewed online. 

Table 1 lists, for each set of simulated spectra, their average signal-to-noise ratio and the difference between the input and the recovered values for the following parameters: luminosity-weighted age, metallicity, and dust attenuation, together with the velocity and velocity dispersion, with their corresponding standard deviation. Here we are measuring two different effects: (i) the systematic deviation from the input parameters, characterized by the average offsets, (i.e., the accuracy of the measurement) and (ii) the precision in the measurement of the parameter (with or without a systematic bias or offset), characterized by the standard deviation. This corresponds to the expected the precision of each individual realization. However, the error of the systematic deviation described before is ≈30 times smaller, taking into account that each row corresponds to 1000 independent simulations. Only when the average offset is clearly larger than this error it is possible to claim that there is a systematic offset for the simulated dataset. This is more evident for the spectra with high S/N, where the systematic effects dominate over the precision of the measurement. However, in most of the simulated datasets there is a clear and well defined bias with a similar trend for each different template, in spite of the estimated precision. Nevertheless, if the standard deviation is clearly larger than the derived offset, the expected bias for an individual realization is too small compared with the systematic bias/offset, and thereforen it is sub-dominant. This is the case for most of the values in the table, and in particular for the spectra with lower S/N.

Table 1 Results of the simulations: single burst 

Each row shows the results of 1000 simulations, corresponding to a different S/N level, for each particular stellar template. Therefore, the errors in the derived offsets are a factor √1000 smaller than the values listed in the table. The fourth column lists the units for the dust attenuation Av in magnitudes, while Colums 5 and 6 list the velocity and velocity dispersion in km/s. The last column contains correlation coefficient between ΔAge/Age and ΔZ/0.02.

There are large differences between the accuracy and precision of the recovered parameters depending on the adopted stellar templates. For those templates with higher spectral resolution, miles72 and gsd156, the parameters are recovered with a better accuracy and precision for the spectra with higher signal-to-noise ratio. In those cases the standard deviation and the offsets be tween the input and output values are smaller, irrespectively of the possible systematic offsets. Those offsets are in general smaller than the standard deviations, and therefore they are sub-dominant for individual realizations, although in general they are statistically significant.

In summary, the simulations indicate that it is possible to recover the properties of the stellar populations for a S/N of at least ≈ 30-50 per pixel to derive accurate and precise results in this particular case. However, for those templates with coarser spectral resolution the parameters are recovered equally well/badly for a wide range of S/N. Figure 8 illustrates these results, showing for a particular set of simulations the comparison between the input and the recovered parameters for each individual spectrum in the dataset. In some cases, like in the age-age comparison, there are clear structures that indicates a possible bias in the recovery of this parameter.

Fig. 8 Results from the simulations. From left to right and from top to bottom, each panel shows the comparison between the input luminosity-weighted age, metallicity, dust attenuation, redshift, and velocity dispersion with the values recovered by our fitting technique. In addition, the degeneracy between different parameters is explored, showing a comparison of the differences between the input and output values for age and metallicity, age and velocity dispersion, metallicity and velocity dispersion, and age and dust attenuation. The results are presented for the gsd156 template and the simulations for a S/N≈60. 

A comparison of the offset between the output and input values for different parameters is also explored to understand the possible degeneracies among them. In particular, we explore the possible bias introduced by the age-metallicity degeneracy: In general older stellar populations are redder, like metal rich ones, while younger and metal poor ones are bluer (Worthey 1994). This in troduces a degeneracy in the derivation of both parameters, as discussed by previous authors (e.g. Sánchez-Blázquez et al. 2011). The strength of this degeneracy is parametrized by the correlation coefficient, r cor, between both parameters, listed in Table 1. The correlation shows a similar behavior for the different stellar templates when S/N¿10. The offset in age and metallicity is in most cases of the order of r cor ≈0.4-0.5.

The age-metallicity degeneracy shows the second strongest correlation. The weakest ones have not been included in the figure. They correspond to the possible degeneracy between the metallicity and both the dust attenuation and the velocity dispersion, and are of the order of r cor ≈0.1. Thus, there is no degeneracy between those parameters. The offset in age present a weak trend with the offset in velocity dispersion, with a typical r cor ≈0.3־. Finally, the strongest correlation is found between the offset of age and the offset of dust attenuation, for which the correlation is strong (r cor ≈0.8) in all the cases explored. This degeneracy is a well known one and indicates how important it is to have a good determination of the dust attenuation to derive the ages. The effect of this degeneracy is stronger for the younger than for the older stellar populations.

It is beyond the goal of this study to conduct a comparison between the properties of the different SSP libraries discussed here. We used different ones to illustrate that the procedure is able to recover the input parameters independently of a particular library. However, it is clear that there are differences. This is the case of the already quoted miles72 and gsd156 templates, when compared with bc138. The accuracy of the derived parameters is very similar in the three cases (as estimated by the bias/offset). However, the precision is much worse for the latter one. The parameters are recovered with far less precision when using the libraries be 138 and mar 136. This may reflect subtle effects of how the parameter space is covered, differences in the stellar templates adopted, or details in the code that generate the templates. For example, the miles72 template does not include SSPs younger than 60 Myr, while the gsd156 template does not include SSPs more metal poor than 0.2 Z ʘ, while the range of parameters covered by bc 138 includes both. Finally, the be 138 template uses the STELIB stellar library (Le Borgne et al. 2003), whose spectrophotometric calibration along the full wavelength range is worse than the one of MILES (Sánchez-Blázquez et al. 2006). We have repeated the analysis based on the be 138 library restrict ing its age/metallicity coverage to that of the miles72 and gsd156. Only in the case of the metal-rich version of the library was there an improvement on the accuracy and precision of the recovered parameters. However it never reached the values found for the two latter libraries. Therefore, although the coverage in the age/metallicity space has an effect on the results of the simulations, biasing towards those libraries with lower metallicity range, it cannot totally explain the reported differences.

Finally, the mar 136 library had the lowest spectral resolution, which strongly affected the results. A detailed discussion of the expected differences when using different stellar templates can be found in (González Delgado et al. 2015). For the other libraries, the stellar age, dust attenuation, and stellar velocity properties were recovered well, although there was a wide range in the accuracy and precision of the properties derived, as summarized in Table 1. For the SSP libraries with the better spectral resolution, and for S/N>50, the precision was of the order of ≈0.1-0.2 dex, similar to the one reported by other fitting procedures (e.g., , Cid Fernandes et al. 2014). The reported uncertainties were slightly larger than the ones estimated for other codes, like S (Ocvirk et al. 2006; Sánchez-Blázquez et al. 2011). Ocvirk et al. (2006) found a typical error Δage = 0.06 0.09 dex and Δmet = 0.10 - 0.14 dex, for signal-to-noise ratios of 50 and 20 respectively, and a spectral resolution R = 1000. However, in the latter simulations the effect of the precision of the velocity dispersion and the dust attenuation was not taken into account. Using a similar code, Sánchez-Blázquez et al. (2011) found that there was an intrinsic uncertainty in the derivation of the parameters associated with the velocity dispersion, a trend already reported by Koleva et al. (2009).

The velocity dispersion was better recovered for those stellar templates with higher spectral resolution (e.g., miles72), and worse for those with lower spectral resolution (e.g., mar136). Unexpectedly, the velocity did not present the same trend. There is an updated version of the Maraston (2005) SSP models with higher spectral resolution, based on different empirical and theoretical spectral libraries (Maraston & Strömbäck 2011). They cover a similar range of ages as the templates outlined here; however, they do not cover the same range of metallicities for the younger stellar ages. In general, they produce results similar to the ones provided by the miles72 and gsd156 models. Therefore, the main differences described here were introduced by the coarse spectral resolution of the Maraston (2005) models. When the updated/higher-resolution version was used, we found no significant differences.

The age-metallicity degeneracy is weaker for the high signal-to-noise simulations for the gsd156 and mar 136 libraries, slightly larger for be 13 8, and clearly stronger for miles 72. However, in all cases the correlation is significant, with a probability of not being produced by a random distribution of the data larger than a 99.99%. In the case of the be 13 8 library the same reason outlined before for the larger error in the parameters recovered could explain this small increase in the degeneracy. In the case of miles72, the fact that this library does not include stellar populations younger than 63 Myr could produce that effect, since the algorithm tries to compensate for the absence of these young stellar populations by giving more weight to the more metal poor ones. Finally, we notice that for some templates the strength of the degeneracy decreases with increasing noise, which indicates that for noisy spectra the random error dominates over this systematic bias (e.g., gsd156 and miles72). However, for other templates, the strength is enhanced by the noise (e.g. mar 136), highlighting the effect of a coarse spectral resolution in the degeneracy between both parameters.

A possible source of error in the interpretation of these simulations is the fact that the same set of templates was used to create the simulated data and to model them. In principle, we assumed that the grid of templates is representative of the real stellar population of the spectra analyzed, but this library is likely to be incomplete. This is not the case in our simulations, by construction. To study the effects of this possible incomplete representation of the observational spectra by the template library, we created a new set of simulated data using the gsd 156 library, that was then fitted using a much reduced version of the template library, called gsd32. This template includes a grid of 32 SSPs, corresponding to 8 ages covering the same range as gsd156, and with the same 4 metallicities. The results from this simulation are listed in Table 1, which shows that in the case of an incomplete sampling of the spectroscopic parameters by the model template, these parameters are still well reconstructed, albeit with a lower precision (as expected). We repeated the experiment with an even more reduced set of stellar libraries, with only 12 SSPs, with a similar result. It is important to note that in both cases the more extreme stellar populations from the original template were retained in the restricted one (i.e., the oldest/youngest and more metal rich/poor). Finally, we remark that although the average parameters that characterize the stellar populations are still well recovered, the details of the SFH and chemical enrichment are lost, obviously.

3.2.2. Results from the simulated SFH

The second set of simulations was created to follow a more realistic SFH. We created three different sets of simulations, each one comprising 1000 simulated spectra with the same normalized flux at 5000 Ǻ, corresponding to 100, 50 and 10 10-16 erg s-1 cm-2 Ǻ-1, respectively. The velocity, velocity dispersion, and dust attenuation were taken into account in a similar way as in the previous set of simulations described in § 3.2.1. For this study we show only the results using the miles72 template library, although we have repeated all the analysis using the other libraries, with similar results.

The fitting procedure was then applied to the simulated spectra using the same steps adopted to analyze the observational data, deriving, for each spectrum, the different data products described in § 2.3. Figure 9 illustrates the result of this exercise for three particular spectra, each one corresponding to a different S/N level. It shows the input coefficients for a particular simulated spectrum in comparison with the output distribution. It is clear that both the SFH, represented by the histogram in the bottom of each panel, was well reproduced for the different simulated spectra, although as expected, the input shape was more and more diluted as the S/N decreased. The metal enrichment pattern was well reproduced in all of the simulations, being more accurate for the older stellar populations than for the younger ones. The distribution of metallicities at different ages and ages at different metallicities, represented in the figure as yellow and orange solid circles, are a relatively good representation of the input correlation between the metals at different ages; thus, they are good tracers of the metal enrichment pattern.

Fig. 9 Top-Left panel: Distribution of input coefficients along the age-metallicity grid for a simulated spectrum corresponding to the second set of simulations, i.e, those resembling a star formation history with a continuous declining star formation rate, using the miles72 template library. For each age-metallicity pair the weight is shown as a solid square, with size and color proportional to the weight of the corresponding SSP age,met. For those SSPs with a stronger contribution, the actual fraction is indicated on top of each square. The horizontal and vertical histograms represent the co-added coefficients for a particular age and metallicity respectively. Top-Right panel: Similar distribution corresponding to the coefficients recovered for a spectrum with a signal-to-noise ratio of ≈100. The yellow filled circles show the average metallicity for a given age, while the orange filled circles show the average age for a given metallicity. Bottom-left panel: Similar distribution as shown in the previous panel but for a simulated spectrum with S/N≈50. Bottom-right panel: Similar distribution as shown in the previous panel but for a simulated spectrum with S/ N≈10. The color figure can be viewed online. 

Table 2 lists, for each set of simulated spectra, their average signal-to-noise ratio and the difference between the input and the recovered values for two of the derived parameters: luminosity-weighted age and metallicity. The kinematics and dust attenuation values are not listed since they did not present significant differences compared to the previous simulations, due to the decoupling of the derivation of these properties and the derivation of the SFH, as described in § 2. We note that the precision of the parameters recovered is better for this more realistic SFH than for the single-burst one, described in § 3.2.1, even for low S/N spectra.

Table 2 Results of the simulations: simulated SFH 

In general, the two luminosity-weighted parameters are recovered with a precision better, or of the order of, 0.1 dex.

3.2.3. Accuracy on the estimated errors

The two main goals of the current modification of FIT3D with respect to previous versions, regarding the analysis of the stellar populations are: (i) to provide a reliable estimation of the star formation history and the average properties describing the stellar population, and (ii) to obtain a useful estimation of the uncertainties of the derived parameters. In previous sections we have il lustrated how well are recovered the parameters that describe the underlying stellar population. In this section we explore how representative are the estimated errors of the real uncertainties of these parameters. For this purpose we used the first set of simulations, described in § 3.2.1, that cover a wider range of parameters, and we compared the errors estimated by the fitting procedure (σPAR) with the absolute difference between the recovered and the input parameters (ΔPAR), where PAR is the parameter considered: age, metallicity, or dust attenuation (Av).

FIT3D estimates the error of each parameter based on the individual values derived along the MC-chain of realizations described in the previous sections. The best fitting model was derived for each MC realization, and therefore, the best set of parameters that describe the input spectrum was obtained. It is generally assumed that the standard deviation of the values derived for each parameter is a good representation of the real errors (σPAR). Indeed, this is the basic assumption of any MC or bootstrapping procedure aimed to estimate the errors in any analysis. However, this assumption holds only if the different parameters are independent, and if the errors affect equally each of them. Therefore, there is no guarantee that σPAR. is indeed a good representation of the real uncertainty in the derivation of the parameter.

Figure 10 shows the distribution of σPAR. along ΔPAR, for the three parameters, and for each of the stel lar template libraries described before. The patterns in the comparison between the estimated and real errors are very similar for the different stellar templates, with most of the differences due to the ranges covered in age and metallicity. It seems that the errors derived by the fitting procedure are a good trace of the real differences between the input and recovered values, for both the luminosity-weighted age and metallicity, with a global underestima tion of the real errors by a factor ≈4, ΔAge ≈ 4σAge, and Δmet ≈ 4σmet. For those libraries with a smaller range of age or metallicity the trend between both errors is less obvious, as expected since the dynamical range of parameters is less well covered.

For the dust attenuation, there is also a bias between the recovered and the real error, but it is not obvious whether the estimated error is a good estimation of the precision with which this parameter is recovered. Indeed, in general it seems that the dust attenuation is recovered with an average uncertainty of ≈0.1-0.3 dex, without a clear correspondence with σ Av. However, we must recall that the error in the dust attenuation was not estimated on the basis of a MC-scheme, but corresponds to the value estimated from the X2 curve. In future versions of FIT3D we will try to implement a MC-scheme also for this critical parameter.

Fig. 10 Distribution of the estimated errors for the different parameters derived for the stellar populations, luminosity-weighted age (top panels), metallicity (middle panels), and dust attenuation (bottom panels), along with the absolute difference between the output and input parameter for the first set of simulations described in § 3.2.1. Each panel, from left to right, corresponds to a different stellar template library considered in the simulations. For each panel the contours correspond to the density of points, with the first contour encircling 95% of the points, and each consecutive contour encircling 15% fewer points. The color figure can be viewed online. 

In summary, the basic assumption that σPAR is a good representation of the real uncertainty does not hold. For the age and metallicity the estimated errors are a factor 4 lower than the real uncertainties, although there is a correlation between both of them and, therefore, it is possible to derive the latter from the former ones in a simple way (in general). However, for the dust attenuation, the errors estimated by FIT3D are unrealiable and there is no simple way to transform them into realistic ones.

3.3. Accuracy of the properties of the emission lines

Following a similar philosophy, we created a set of simulations in order to derive the accuracy and precision of the recovery of the properties of the emission lines. The first set of simulations was created assuming purely Poissonian noise, and without taking into account the possible effects of the imperfect subtraction of the underlying stellar population, and other sources of non Poissonian noise. We simulated three different systems of emission lines within this set: (i) Hβ and [ OIII ] λ4959,5007, with a flux intensity of 100, 66, and 22 10-16 erg s-1 cm-2 Ǻ-1, respectively, and covering the wavelength range 4800-5050 Ǻ, in the rest-frame; (ii) Hα, [ NII ] λ6548,6583, and [ SII ] λ6717,6731, with a flux intensity of 100, 66, 22, 25, and 25 10-16 erg s-1 cm-2 Ǻ-1, respectively, covering a wavelength range 6400-6900Ǻ; and (iii) the same emission lines included in the previous simulation plus a broad Hα component, covering the same wavelength range.

All the emission lines, except the broad component included in the last set, were simulated assuming the same range of velocity dispersion used in the simulations of the stellar populations, σ = 2.5Ǻ and 7.5Ǻ v 100 - 400 km/s). The velocity dispersion was then selected randomly within this range. Therefore, while in Case (i) the emission lines were always deblended, in Case (ii) there were instances in which the emission lines were deblended, but in many other cases they showed different levels of blending. Finally, in Case (iii) the narrow emission lines presented the same degree of blending than in Case (ii), but in all the cases they were blended with the broad emission line. This affects the accuracy and precision of the recovered emission line properties. For the broad component included in Case (iii) we considered a range of velocity dispersion σ =32.5 - 37.5 Ǻ v ≈ 1500 - 1700 km/s).

The velocity was randomly selected around 6000 ± 500 km/s, for the narrow emission lines. The broad component was randomly redshifted by 200 km/s from this velocity, which is a typical situation in the case of broad components generated by either outflows or type one AGNs. The emission line spectra were sampled with a pixel scale of 2Ǻ/pix, similar to the one provided by the CALIFA V500 data set. The results can be extrapolated to other instrumental configurations taking into account the ratio between the velocity dispersion and the size of the sampling pixel. The effect of the instrumental configuration has not been taken into account. It affects mostly the recovering of the physical velocity dispersion, which is a quadratic subtraction of the instrumental one from the measured dispersion. Therefore, it basically depends on this later parameter, although not in a linear way.

Figure 11 shows the comparison between the recovered and the simulated parameters for three sets of simulated emission lines, as a function of the signal-to-noise ratio of the emission lines. The S/N was estimated as the ratio between the peak intensity of the emission line and the noise level included in the simulated spectrum. As expected, there is a clear dependence with S/N (for S/N<10), with the worst recovered parameters corresponding to the lower S/N emission lines. For large S/N values (>20) there is almost no dependence of the uncertainty with the S/N. There is also a clear difference between the three sets of simulations, with the better precision corresponding to Case (i), and the worst to Case (iii). This is not an effect of the signal-to-noise, since in Case (i) and Case (ii) all the emission lines have the same integrated fluxes and a similar flux density per spectral pixel for a fixed dispersion.

Fig. 11 Relative difference between the recovered and the simulated parameters of the emission lines (top: flux intensity, middle: velocity dispersion, and bottom: velocity) for the three sets of simulated emission lines corresponding to [OIII ]+Hβ (left panels). [NII ]+Hα (central panels), and this latter one with an additional broad component for Hα (right panels), as a function of the signal-to-noise ratio of the emission lines. For each panel the contours correspond to the density of points, with the first contour encircling 95% of the points, and each consecutive contour encircling 15% fewer points. The color figure can be viewed online. 

A summary of these results is included in Table 3. For each set of simulated emission lines, it lists the mean and standard deviation of the relative difference between the recovered and simulated parameters for different ranges of signal-to-noise ratio. The program is able to recover the intensity of the emission lines with a precision better than ≈10% for weak emission lines (S/N≈3-10), in all the cases. The precision becomes worse, ≈20%, for emission lines swimming in the noise, although the accuracy is equally good. However, at this low S/N values the systematic effects (in the subtraction of the underlying stellar population and other sources of non Poissonian noise not considered in this simulation) will dominate, and therefore this precision should be considered as an upper limit. It is interesting to notice that although the precision increases towards higher signal-to-noise ratios, it never gets better than ≈3%, at least for the highest S/N ratios explored in this simulations. We consider that this threshold is related to the fact that we fit together the intensity, velocity, and velocity dispersion, and their derivation is somehow correlated. Therefore, the S/N of the integrated flux alone may not be a good parameter to characterize how well the emission lines are detected.

Table 3 Results of the simulations: emission lines 

1∆vsys is multiplied by ten to show it in the same scale as the other two parameters.

The velocity and velocity dispersion are also well recovered. In the case of the velocity dispersion it has a precision similar to that of the emission line intensity; it is recovered with a precision better than 10% in almost all the simulated spectra. The best precision for this parameter is ≈2%, even for high S/N spectra. Finally, the velocity is recovered with a very high precision, ≈10 times better than the other two parameters. The precision ranges from ≈12-54 km/s at the average redshift of the simulated spectra. This precision is between two and ten times better than the simulated velocity dispersion. These simulations confirm a widely accepted result indicating that the velocity derived using emission lines can be recovered with a precision as good as a few per cent of the velocity dispersion of the emission lines (typically 10%).

As already noted, the precision of the emission line parameters recovered is different for each set of simulations, being, in general, worst for the simulations with blended emission lines. The degradation of the precision could be as bad as two times the precision when the emission lines are not blended.

3.3.1. Interdependence of the accuracy between parameters

As indicated before, the precision in the estimation of one of the parameters may affect the accuracy with which others are recovered. In particular, the velocity dispersion could affect the recovered flux intensity, since this parameter is an integral that depends on the former one. Figure 12 shows the relative accuracy of the recovered flux intensity as a function of the accuracy of the recovered velocity dispersion for the three sets of simulations described before. In all three cases the distributions show a cloud of points with a dispersion of the order of ±0.1 Δ par/par. However, it is clear that the distribution is broader and more roundish for the simulations including blended emission lines than for the simulations including only well resolved emission lines (left panel). This is a consequence of the larger errors/inaccuracies in the derivation of both the flux intensity and the velocity dispersion for the blended emission lines, as shown in § 3.3. In this analysis we used all the simulations described before irrespectively of the S/N. We repeated the analysis for different S/N ranges without significant differences. Only for the very low S/N (<10) regime we found a broadening of the distribution, but the main trends held.

Fig. 12 Relative difference between the recovered and the simulated intensity of the emission lines for the three sets of simulated emission lines corresponding to [O III ] +Hβ (left panel), [NII ] +Hα (central panel), and this latter one with an additional broad component for Hα (right panels), as a function of the relative difference between the recovered and the simulated velocity dispersion of the corresponding emission lines. For each panel the contours correspond to the density of points, with the first contour encircling 95% of the points, and each consecutive contour encircling 15% fewer points. The color figure can be viewed online. 

There is a clear trend between both relative differences, stronger for the simulations with unblended emission lines, with a correlation coefficient r =0.41 (p > 99.99% of not being produced by a random distribution of the data). This indicates that the bias introduced in the estimation of the flux intensity by the inaccuracy in the derivation of the emission lines is stronger for clearly blended emission lines. However it is not a dominant factor in the error budget, for the signal-to-noise range of our simulations.

3.3.2. Accuracy of the estimated errors for the emission lines

As in the analysis of the stellar population, it is important to determine whether the estimated errors are representative of the real ones. Following the previous analysis, we compare for each emission line in each of the three simulated datasets the estimated error with the real ones, i.e., the absolute difference between the value recovered and the simulated (input) one. Figure 13 shows the results of this comparison.

Fig. 13 Relative error estimated for the recovered parameters of the emission lines (top: flux intensity, middle: velocity dispersion, and bottom: velocity) for the three sets of simulated emission lines corresponding to [OIII ] +Hβ (left panel), [N II ]+Hα (central panel), and this latter one with an additional broad component for Hα (right panels), as a function of the relative difference between the recovered and the simulated value of this parameter for the corresponding emission line. For each panel the contours correspond to the density of points, with the first contour encircling 95% of the points, and each consecutive contour encircling 15% fewer points. The red dashed line in each panel corresponds to the one-to-one relation. The color figure can be viewed online. 

As shown in previous sections, the real error in the derivation of the parameters is, in general, larger for the simulations including blended emission lines than in those including only un-blended ones. For the flux intensity, the estimated errors seem to be a good representation of the real ones for real errors larger than ≈4-5%. For errors smaller than this value the program seems to overestimate the errors, and derives a minimum value of ≈5%, regardless of the real difference between the input and output parameters. However, since in most simulations the emission intensities cannot be recovered with a precision better than a few percent, this lower limit seems still reasonable.

For the velocity dispersion the estimated errors seem to be of the same order of the real ones, at least for the simulations that did not include broad emission lines. In general, the errors seem to be slightly overestimated, with a similar trend towards a minimum value of the error of around ≈5%. For the simulations including broad emission lines the estimation of the error shows a wide scatter, with a clear underestimation of the error in many cases. A detailed investigation shows that the larger variations in the estimation of the error correspond to the broad emission line component included in this simulation. It is interesting to note that, in any case, the dispersion seems to be well recovered (see Figure 11), but for broad emission lines the program tends to underestimate the errors.

Finally, the estimated errors for the velocities are an order of magnitude lower than the errors derived for the other two parameters. They are just slightly overestimated compared to the real difference between the input and recovered velocities, being in general a good representation of the real errors. As in previous cases, the estimated errors are more accurate for the simulations including only narrow emission lines.

3.4. Simulations using real data

So far we have estimated the accuracy and precision of the parameters using simulations that only take into account the Poissonian noise. To investigate further the precision in the derivation of the parameters we need to take into account all the different sources of error that are included in real data.

For these particular simulations we used the data sets provided by the CALIFA survey. In particular, we adopted the V500 data cube corresponding to NGC 2916, an Sbc galaxy at z ≈0.01244, which is an average galaxy in this survey. This data cube includes 78x72 individual spectra, of which ≈2500 contain data and the rest are outside of the original hexagonal pattern of the instrument. We note that the simulations do not depend on the particular galaxy selected. We used this one because it presents both strong and weak emission lines within the FoV of the data cube, thus covering a wide range of signal-to-noise ratios.

The simulations were created using as input the results from the fitting using P 3D on this particular dataset. The analysis with P 3D will be detailed in forthcoming articles of this series, although some examples of the data products have been already presented (Sánchez et al. 2013). For the purpose of this study we should know that the pipeline performed a spatial binning over the data cube, aggregating adjacent spaxels based on their difference in relative intensity (at a reference wavelength); it provided final row-stacked spectra (RSS) comprising ≈500-1000 individual spectra with an increased signal-to-noise (in the continuum) with respect to the individual spaxels. Finally, each of the individual spectra was analyzed using the current version of FIT3D to provide all the data products regarding the stellar population and the emission lines described in previous sections. In this particular analysis we adopted the results derived using the gsd156 stellar library template.

The resulting star formation history, stellar kinematics, stellar dust attenuation, and the corresponding data products of the emission lines (intensity, velocity, and velocity dispersion) can be used as input parameters for a simulation with more realistic parameters. For each spatial bin we created a spectrum with a model stellar continuum and a model set of emission lines, using as input the output of the fit over the real data. Then, we renormalized the flux in each individual spaxel within each bin following the relative contribution of the flux to the bin by this spaxel in the original data. In this way we obtained a realistic model data cube with the same spatial distribution of the stellar population and ionized gas as the real data.

Then, we used the residuals from the fitting process over the observed data cube to generate a realistic distribution of the noise, which includes not only the Poissonian noise, but also any other source of error. To do so, we used the residual data cube obtained from the analysis with FIT3D, and, for each spaxel at each wavelength, we co-added a random value within a range plus-minus the median of the absolute value of the noise distribution within a box of ±3". In this way we respected the amplitude of the real noise at each wavelength and at each position within the data cube, excluding the peaks in the noise distribution (like the ones produced by cosmic-rays or inaccuracies in the sky subtraction). We should clarify that this derivation of the noise does not follow either a Gaussian or a square distribution.

We then applied the same spatial segmentation to extract the RSS-spectra at the same spatial locations as in the case of the original data cube. After that, we processed the spectra using the same fitting routines to recover both the properties of the stellar population and emission lines. We adopted the same stellar population library for the analysis, for a better comparison be־ tween the input and output parameters. We analyzed four sets of emission lines, with coupled kinematics, as described in § 2.4. Those emission lines were analyzed in four wavelength ranges: (1) [OII ]λ3727 (3700-3750Ǻ); (2) Hβ and [O III ]λ4959,5007 (4800-5050Ǻ); (3) Hα and [NII ]λ6548,6583 (6680-6770); and(4) [SII ]λ6717,6731. Altogether, they comprised more than 5000 individual emission lines, with different line intensities, in a wide range of wavelengths. Therefore, they were affected in a different way by the imperfect subtraction of the underlying stellar population and the wavelength dependence depth of the data. All the emission lines considered here were deblended at the spectral resolution of the V500 CALIFA data, however, some of them could be affected by the wing of the adjacent emission lines, in particular the doublet [SII ]λ6717,6731.

3.4.1. Accuracy of the properties recovered for the stellar populations

Figure 14 shows the results of comparing the data products derived for the stellar populations as a function of the corresponding input/simulated values, for spectra with S/N higher than 20. As in the case of Figure 8, it demonstrates that there is a good correspondence between the input and the recovered parameters, with a clear one-to-one correlation in all cases. A more quantitative comparison is included in Table 4, which shows the average S/N and the difference between the input and the recovered values, for different ranges of S/N. Despite the different procedure adopted, the results are very similar to the ones shown in Tables 1 and 2. For spectra with a S/N higher than ≈20, most parameters were recovered with a precision of ≈10%. In the case of the velocity the precision is an order of magnitude better. It is interesting to note that the velocity was recovered with a better precision when using a more realistic SFH than a purely simulated one.

Fig. 14 Results from the simulations based on the CALIFA V500 data cube of NGC2916. Each panel shows, from left to right and from top to bottom, the comparison between the input luminosity-weighted age, metallicity, dust attenuation, redshift, and velocity dispersion with the values recovered by our fitting technique. The results are presented for the gsd15 6 template. The final bottom right panel illustrates the degeneracy between the recovered age and metallicity, showing a comparison between the input-output differences for both parameters. 

Table 4 Results of the simulations: simulated NGC2916 data cube 

As in the previous cases, we explored the possible effect of the age-metallicity degeneracy in these more realistic simulations. The last panel of Figure 14 shows a comparison of the offset between the output and input ages, and the offset between the output and input metallicities of the stellar populations. There is a weak trend between the two offsets, in the expected direction induced by this degeneracy. The strength of this degeneracy is parametrized by the correlation coefficient between the two parameters, listed in Table 4. The trend is weak, and affects relatively more the metallicity than the age. This is expected since the relative error in the derivation of age is always smaller than the error in the derivation of metallicity. As expected, this degeneracy, a systematic effect, is enhanced by the S/N of the data. The effect is stronger for the more metal poor stellar populations, since the offset range is ±0.5 Z/Zʘ, i.e., ΔZ =±0.01. For the metal-rich population the effect is marginal, but for the stellar populations with sub-solar abundances it could bias the result by a factor of 2 or more.

We note here that the degeneracy and the range of offsets between the input and output parameters are larger than the values reported in previous studies (e.g. Ocvirk et al. 2006; Sánchez-Blázquez et al. 2011), as we indicated before, but they are of same order as the ones reported by more recent simulations (e.g. Cid Fernandes et al. 2014). Maybe the reason for this disagreement is that, although in the previous cases a full SFH and reasonable chemical enrichment processes were considered, they did not explore a wide range of possible uncertainties in the dust attenuation and the velocity dispersion. They do show a hint of its effect though (e.g., Figure B3 of Sánchez-Blázquez et al. 2011).

Figure 15 shows a comparison between the input and output coefficients of the decomposition of the stellar population with the gsd156 template library forthe spectra extracted using a 5" diameter aperture on the real and simulated CALIFA data cube. Both distributions have a very similar pattern, with the oldest stellar population being metal-poor (sub-solar metallicity), and a rapid chemical enrichment in the early epochs of the Universe. The star formation history presents several bursts, with most of the stars formed in the first 1-2 Gyr of the Universe. This pattern is very similar to the one described by Perez et al. (2013). The younger stellar populations (<1 Gyr) are more metal-poor, which could be interpreted as a result of infall of more pristine gas into the inner regions, thus, explaining the remanent star formation events in the very center of this galaxy. However, it could well be that these low intensity coefficients are just a product of the noise. In any case, the code is able to recreate realistic SFHs and chemical enrichment patterns, whatever their nature.

Fig. 15 Top panel: Distribution of coefficients along the age-metallicity grid derived from the analysis using FIT3D of the central 5" diameter aperture spectrum extracted from the CALIFA V500 data cube of NGC2916, used as input coefficients in the simulated spectra described in the text. Bottom panel: Similar distribution corresponding to the recovered coefficients for the simulated spectrum, resulting from the fitting procedure. All the symbols are similar to those described in Figure 9

Fig. 16 Comparison between the derived parameters of the stellar population for 448 central spectra extracted from the CALIFA V500 data cubes analyzed using FIT3D, S , S , and absorption line indices using the same template of stellar populations. The top panels show the comparison between the luminosity-weighted ages, and the bottom panels show the comparison between the luminosity-weighted metallicities. In all the panels the size of the symbols is inversely proportional to the S/N of the spectrum analyzed. The error bars indicate the largest and average errors in each of the parameters, as derived by FIT3D. In each panel the inset box shows the normalized histograms of the differences between the two parameters compared. 

3.4.2. Accuracy of the properties recovered for emission lines

Figure 17 shows the results of the comparison between the properties of the simulated emission lines and the values recovered for the simulation based on the CALIFA V500 data cube for NGC 2916. A more quantitive summary is presented in Table 4, where the average and standard deviation of the differences between the input and output values of the different parameters are presented for different ranges of the signal-to-noise ratio of the emission lines. Being a realistic simulation based on a galaxy with mild emission lines, the dynamical range for the signal-to-noise ratio is smaller than in the case of the simulations presented in § 3.3. The results are consistent in general. As expected, the accuracy and precision of the derivation of the different parameters is slightly lower for this simulation, which is more realistic in terms of the signal-to-noise ratio and the propagation of the imperfect subtraction of the underlying stellar population and other sources of non-Poissonian noise. For the well detected emission lines, S/N> 10, the flux intensity is recov ered within 10-15%, the velocity dispersion is recovered within ±10 km/s, and the velocity is recovered within 30-40 km/s. The precision decreases for the emission lines detected with low S/N (< 10), and the results are unreliable for barely detected emission lines.

Fig. 17 Relative difference between recovered and input values for the parameter analyzed for the emission lines (left: flux intensity, center: velocity dispersion and right: velocity) for the simulation based on the CALIFA V500 data of NGC2916 as a function of the emission lines flux intensities. For each panel the colored images and contours correspond to the density of points. with the first contour encircling 95% of the points. and each consecutive contour encircling 15% fewer points. The color figure can be viewed online. 

4. COMPARISON WITH OTHER FITTING ROUTINES

As indicated in previous sections, FIT3D is obviously not the first package developed with the aim of recovering the properties of the stellar population from optical spectra of galaxies. So far, we have illustrated how well it recovers the main parameters that define the stellar populations, based on different kinds of simulations. In this section we explore how the recovered parameters compare to similar ones derived using other fitting tools.

We require a set of spectra with enough S/N in the continuum to minimize the effects of noise in the comparison. Therefore, we selected the central spectra extracted from the 448 CALIFA galaxies observed using the V500 setup up to September 2014, co-adding the spaxels within a diameter of 5" around the peak intensity in the V-band. All those spectra have a S/N>50 per spectral pixel, as shown in Sánchez et al. (2012 b). The wide range of galaxy masses and morphologies covered by CALIFA guarantees that with those spectra we explore a wide range of stellar populations. We analyzed them with FIT3D, using the gsd 156 template library, and the procedures described in the previous sections.

Table 5 Results of the simulations: emission lines in the NGC2916 data cube 

1Kinematic parameters are given in km/s.

Using the same data set and the same library we applied, in addition to FIT3D, (i) (Cid Fernandes et al. 2005), following the prescriptions of Cid Fernandes et al. (2013), (ii) S (Ocvirk et al. 2006), following the prescriptions of Sánchez-Blázquez et al. (2011), Sánchez-Blázquez et al. (2014a), and Sánchez-Blázquez et al. (2014b), and (iii) the analysis based on the study of the absorption line indices described in Sánchez et al. (2007b) and Sánchez et al. (2011). In summary, after subtracting the emission lines, we derive for each spectrum a set of stellar indices (Hσ, Hβ, Mgb, Fe5270, Fe5335, and D4000), and we estimate the ages and metallicities based on a comparison with the expected values for each pair of indices with the corresponding ones for the same stellar template, using the program (Cardiel et al. 2003). It is beyond the scope of this article to describe in detail the different codes. However, we should note that they rely on very different assumptions and inversion algorithms. In essence, is more similar to FIT3D since it does a pseudo-random exploration of the parameter space based on a Markov chain method. However, as indicated before FIT3D does not use a Markov chain, but performs a MC realization process and a linear inversion to derive the weights of the stellar decomposition. Neither of these two codes makes any a priori assumption on the shape of the star-formation history. On the contrary, S performs a derivation of the star formation history using a peak-to-peak fitting algorithm (Ocvirk et al. 2006), that allows to introduce a Bayesian prior on the shape of the star formation history (e.g. Sánchez-Blázquez et al. 2014a). Finally, the analysis of the stellar indices presents larger uncertainties for a similar S/N level, as demonstrated by Sánchez-Blázquez et al. (2011). However, it is the method that requires fewer assumptions and computationally it is the simplest: It is a derivation, not a fitting, and performs a comparison with a single SSP model, not with a combination of them. An additional difference is that both the indices and S do not take into account the effects of dust attenuation, since they remove the continuum (which may have additional advantages if the spectrophotometric calibration has problems of any kind). In the case of the line indices any dust attenuation affects the age and metallicity derivation, since a priori it is assumed that it may affect D4000. Despite all these differences, the main aim of the different methods is to provide with a realistic estimation of the parameters that define the properties of the stellar population.

Figure 16 illustrates the results of this comparison. It shows the derived luminosity-weighted ages and metallicities for the different spectra as obtained by FIT3D, , S ,and the analysis of the stellar indices. The correspondence between the luminosity-weighted ages is very good: (Δage = -0.06±0.15 dex), S (Δage =0.03±0.22 dex), and stellar indices (Δage =0.06±0.15 dex). However, the distribution is compatible with a one-to-one relation taking into account the estimated errors only for the comparison with , with X 2 /v =1.14. The largest discrepancies correspond to the spectra with the youngest stellar populations, in particular for S. These populations are also the ones that present the strongest dust attenuations, and maybe this is the reason for the discrepancy. In more detail, the agreement between the ages seems to be better with than with S, which presents a systematic difference towards younger stellar populations campared with FIT3D.

The agreement between the metallicities is clearly better for (Δmet = -0.03±0.10 dex), than for the other two methods, S (Δmet =0.03±0.22 dex), and stellar indices (Δmet =0.10±0.27 dex), when compared with FIT3D. Indeed, the distribution is compatible with a one-to-one relation taking into account the estimated errors only for the comparison with with X2/v =1.08. In particular this happens in the range of more metal rich stellar populations, for which S collapses the values towards the highest metallicities in the SSP libraries. This may reflect the simi larities in the procedure between FIT3D and and the differences between both of them and S In the case of the stellar indices the agreement is better for the range of more metal rich stellar populations, for which the absorptions sensible to the metallicity are more prominent.

As a sanity check we perform the same comparison for the values derived between both and S , and we find similar results. The dispersion in the differences between the parameters derived is Δage =0.10±0.26 dex, and Δmet = -0.10±0.14 dex, respectively, confirming the trend described when comparing both procedures with FIT3D. Curiously, the values derived with FIT3D seem to correspond somehow to the average/intermediate values derived by both methods. If the discrepancies were dominated by systematic effects or by differences between the three methods, this should not be the case (appart from the effect in the metallicity when comparing with S ). For the stellar indices we found that the differences in ages and metallicities correlated with dust attenuation, a clear indication that this method tries to compensate the effects of dust by as suming an older and more metal-rich stellar population.

4.1. Comparison of the stellar velocity dispersion

FIT3D provides the main parameters of the stellar populations, including the stellar kinematics. In general, the stellar velocity is very well recovered by most of the stellar population analysis procedures with a good accuracy and precision. However, the stellar velocity dispersion is not recovered equally well, mostly due to the degeneracy between velocity dispersion and metallicity (e.g. Sánchez-Blázquez et al. 2011).

It is widely accepted that one of the best procedures to analyze the stellar kinematics is PXF (Cappellari & Emsellem 2004b). This routine allows to recover with high accuracy the line-of-sight velocity distribution (LOSVD) from absorption line spectra, permitting to model it not only with a simple Gaussian, characterized by its velocity and velocity dispersion, but with more complex profiles for the LOSVD, including asymmetries described by Hermite functions and/or multiple Gaussian-like velocity components. This is the reason why PXF is widely used in stellar population analyses, although additional tools are needed to derive other main properties of those stellar populations (e.g. Sánchez-Blázquez et al. 2014 a,b). The current capabilities of PXF regarding the analysis of stellar kinematics are far beyond the procedures implemented in the improved version of FIT3D. However, it is worth comparing the derived parameters, in particular the stellar velocity dispersion, to illustrate the precision of our current estimations.

Figure 18 shows the comparison between the stellar velocity dispersion derived using FIT3D and the one derived using PXF for the set of spectra described in § 4. There is a very good agreement between the estimations done by both procedures σ = -9.8±36.7 km/s), with a difference compatible with the expected precision for the given spectral resolution (R ≈850, Sánchez et al. 2012b). Therefore, despite the fact that the current version of FIT3D is not optimized for the analysis of the kinematic properties of stellar populations, it provides an accurate estimation of the stellar velocity dispersion.

Fig. 18 Comparison between the stellar velocity dispersion derived for 448 central spectra extracted from the CALIFA V500 data cubes analyzed using FIT3D and PFX. The inset box shows the normalized histograms of the differences between the two velocity dispersion measurements. 

In summary, we can conclude that the parameters of the stellar populations derived with FIT3D are consistent with those derived by other algorithms, with the difference that the new algorithm provides a reliable estimation of the errors of those parameters.

5. SUMMARY AND CONCLUSIONS

In this article we present a new version of FIT3D, a package developed to analyze the optical spectra of galaxies, in order to derive the main properties of their stellar population and the ionized gas. FIT3D implements a new MC iteration method that provides not only the numerical values of the parameters, but also an estimation of their errors. In this regard the procedure is a substantial modification of the previous version of the package.

The main results of this study are the following:

  • The fitting routine was confronted with simulations that demonstrate its ability to recover the properties of the underlying stellar population and the emission lines. We found that the algorithm is able to recover the luminosity-weighted ages and metallicities with a precision of ≈0.2 dex for spectra with a S/N per pixel better than ≈50, for most of the stellar templates analyzed in the simulations. We found systematic differences/offsets in the recovered pa rameters. However, they were sub-dominant for the individual realizations, and only significant in a statistical way. We recovered the well known age -metallicity - dust attenuation degeneracies. In general, the degeneracies were of the order of the precision of the parameters recovered. The use of different templates with different spectral resolutions and coverage of the parameter space affected the accuracy and precision of the recovered parameters in a similar way as a decrease of the S/N. In this regard, it is recommended to use ad hoc simulations in order to determine which stellar templates are optimal for each particular data set. In general we cannot recover the luminosity-weighted age, metallicity, and dust attenuation with a precision better than ≈0.1 dex, even for the highest signal-to-noise.

  • The algorithm is not tuned to derive accurate/precise stellar kinematics. However, it provides a good estimation of the stellar velocity and velocity dispersion, which, at a typical spectral resolution of R ≈1000, is never better than ≈25 km/s, for the highest S/N explored. The currently implemented approach to derive the kinematics introduces a systematic effect in the profile of the subtracted absorption lines. This translates into an unrealistic emission that is at most of the order of EW(Hα)≈0.4Ǻ, far below the typical detection limit for most of spectroscopic surveys (e.g., Gomes et al., in prep., Papaderos et al. 2013).

  • The program is able to recover the star formation and chemical enrichment histories for the galaxies when the signal-to-noise ratio per pixel is large enough (S/N>50). The precision of the recovered ages and metallicities is better when the SFH includes correlations between the distribution of the stellar populations than in the case of a simple dis tortion of a single burst of fixed age and metallicity.

  • The simulations indicate that the nominal errors derived by FIT3D for the parameters that characterize the stellar populations are a factor of four smaller than the estimated errors derived from the comparison of the input and estimated values for the stellar ages, metallicities, and dust attenuation.

  • A detailed set of simulations allowed us to determine the ability of the algorithm to recover the intensity, velocity dispersion, and velocity of the emission lines in three main cases: isolated lines, blended lines, and blended lines with an underlying broad emission line component. In general, it was possible to estimate the parameters with a precision of 0.05-0.10 dex for S/N above 3σ. As expected, the precision of the derived parameters was worse for the more blended emission lines and for the lower S/N.

  • A more realistic set of simulations was created based on real spectra of a typical galaxy extracted from the CALIFA survey data cubes, by using as input for the simulations the output of the analysis of this galaxy, based on the same algorithm, and perturbed by a realistic realization of the noise, based on the residuals of the analysis described. These simulations confirmed that the properties of the stellar populations can be recovered with a precision of the order of ≈0.1 dex only for S/N≈50.

  • The stellar parameters derived using FIT3D were compared with similar ones derived using widely used fitting routines, like and S and with those derived using classical stellar indices, for a set of high S/N spectra extracted from the CALIFA data cubes. This comparison showed that the parameters derived by FIT3D were compatible with those derived using the other tools. In particular, there was a very good agreement in the luminosity-weighted ages for the three cases, with a typical difference of ≈0.15-0.22 dex. The agreement was slightly worse for the luminosity-weighted metallicities in the case of S and the stellar indices, being better for the case of , with a typical difference of ≈0.1 dex.

In summary, we demonstrated that the new algorithm implemented in FIT3D is able to recover the main properties of the stellar populations and emission lines for moderate resolution spectra of galaxies. The tool is fully prepared to analyze IFS data, and, in particular, it is well suited for the analysis of IFU data from the major on-going surveys, like CALIFA (Sánchez et al. 2012b), MaNGA (Bundy et al. 2015), and SAMI (Croom et al. 2012). In forthcoming articles we will present a practical implementation to analyze spatially resolved spectroscopic data and a study of the possible uncertainties.

Acknowledgments

We acknowledge the referee for his/her comments and suggestions that helped to improve the manuscript. SFS thanks the director of CEFCA, M. Moles, for his sincere support to this project. We thanks C.J. Walcher and R. González Delgador for their valuables comments and suggestions that improved this manuscripts in many senses. SFS thanks the CONACYT-125180 and DGAPA-IA100815 projects for providing support in this study. We acknowledge support from the Spanish Ministerio de Economía y Competitividad, through projects AYA2010-15081 and AYA2010-10904E, and Junta de Andalucía FQ1580. EP acknowledges support from the IA-UNAM and from the Guillermo Haro program at INAOE. This study used data provided by the Calar Alto Legacy Integral Field Area (CALIFA) survey (http://califa.caha.es/). CALIFA is the first legacy survey performed at Calar Alto. The CALIFA collaboration would like to thank the IAA-CSIC and MPIA-MPG as major partners of the observatory, and CAHA itself, for the unique access to telescope time and support in manpower and infrastructure. The CALIFA collaboration also thanks the CAHA staff for the dedication to this project. This research is based on observations collected at the Centro Astronómico Hispano Alemán (CAHA) at Calar Alto, operated jointly by the Max-Planck-Institut für Astronomie and the Instituto de Astrofísica de Andalucía (CSIC).

REFERENCES

Aaronson, M., Cohen, J. G., Mould, I, & Malkan, M. 1978, ApJ, 223, 824 [ Links ]

Barrera-Ballesteros, J. K., Falcón-Barroso, J., García-Lorenzo, B., et al. 2014, A&A, 568, A70 [ Links ]

Bruzual, G. & Chariot, S. 2003, MNRAS, 344, 1000 [ Links ]

Bundy, K., Bershady, M. A., Law, D. R., et al. 2015, ApJ , 798, 7 [ Links ]

Calzetti, D. 2001, FASP, 113, 1449 [ Links ]

Cappellari, M. & Emsellem, E. 2004a, PASP, 116, 138 [ Links ]

__________. 2004b, PASP, 116, 138 [ Links ]

Cardelli, J. A., Clayton, G. C, & Mathis, J. S. 1989, ApJ , 345, 245 [ Links ]

Cardiel, N., Gorgas, J., Sánchez-Blázquez, P., Cenarro, A. J., Pedraz, S., Bruzual, G, & Klement, J. 2003, A&A , 409, 511 [ Links ]

Chabrier, G. 2003, PASP , 115, 763 [ Links ]

Chariot, S. & Fall, S. M. 2000, ApJ , 539, 718 [ Links ]

Cid Fernandes, R., González Delgado, R. M., Garcia Benito, R., et al. 2014, A&A , 561, A130 [ Links ]

Cid Fernandes, R., Mateus, A., Sodré, L., Stasińska, G., & Gomes, J. M. 2005, MNRAS , 358, 363 [ Links ]

Cid Fernandes, R., Pérez, E., Garcia Benito, R., et al. 2013, A&A , 557, A86 [ Links ]

Courteau, S., Cappellari, M., de Jong, R. S., et al. 2014, RvMP, 86, 47 [ Links ]

Croom, S. M., Lawrence, J. S., Bland-Hawthorn, J., et al. 2012, MNRAS , 421, 872 [ Links ]

Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., Riccia rdelli, E., Cardiel, N., Cenarro, A. J., Gorgas, J., & Peletier, R. F. 2011, A&A , 532, A95 [ Links ]

Gerhard, O. 2001, The Mass of Galaxies at Low and High Red-shift, in Proceedings of the European Southern Observatory and Universitäts-Sternwarte München Workshop Held in Venice, ed. R. Bender & A. Renzini (Springer Berlin, Heidelberg: EAS), 62 [ Links ]

Gil de Paz, A. & Madore, B. F. 2002, AJ, 123, 1864 [ Links ]

González Delgado, R. M., Cerviño, M., Martins, L. P., Lei-therer, C, & Hauschildt, P. H. 2005, MNRAS , 357, 945 [ Links ]

González Delgado, R. M., García-Benito, R., Pérez, E., et al. 2015, A & A , 581, 103 [ Links ]

González Delgado, R. M., Pérez, E., Cid Fernandes, R. et al. 2014, A&A , 562, A47 [ Links ]

Koleva, M., Prugniel, P., Bouchard, A., & Wu, Y. 2009, A&A , 501,1269 [ Links ]

Le Borgne, J.-F., Bruzual, G., Pelló, R., Lançon, A., Rocca-Volmerange, B., Sanahuja, B., Schaerer, D., Soubiran, C., & Vílchez-Gomez, R. 2003, A&A , 402, 433 [ Links ]

MacArthur, L. A., Courteau, S., Bell, E., & Holtzman, J. A. 2004, ApJ S, 152, 175 [ Links ]

MacArthur, L. A., González, J. J., & Courteau, S. 2009, MNRAS , 395,28 [ Links ]

Maraston, C. 2005, MNRAS , 362, 799 [ Links ]

Maraston, C. & Strömbäck, G. 2011, MNRAS , 418, 2785 [ Links ]

Mármol-Queraltó, E., Sánchez, S. F., Marino, R. A., Mast, D., Viironen, K., Gil de Paz, A., Iglesias-Paramo, J., Rosales-Ortega, F. F., & Vilchez, J. M. 2011, A&A , 534, A8 [ Links ]

Martins, L. P., González Delgado, R. M., Leitherer, C., Cerviño, M., & Hauschildt, P. 2005, MNRAS , 358, 49 [ Links ]

Oconnell, R. W. 1976, ApJ , 206, 370 [ Links ]

Ocvirk, P., Pichon, C., Lançon, A., & Thiebaut, E. 2006, MNRAS , 365, 46 [ Links ]

Panter, B., Heavens, A. F., & Jimenez, R. 2003, MNRAS , 343, 1145 [ Links ]

Papaderos, P., Gomes, J. M., Vilchez, J. M., et al. 2013, A&A , 555, L1 [ Links ]

Pérez, E., Cid Fernandes, R., González Delgado, R. M., et al. 2013, ApJ , 764, L1 [ Links ]

Rix, H.-W. & White, S. D. M. 1992, MNRAS , 254, 389 [ Links ]

Rosales-Ortega, F. F., Kennicutt, R. C., Sánchez, S. F., Díaz, A. I., Pasquali, A., Johnson, B. D., & Hao, C. N. 2010, MNRAS , 405, 735 [ Links ]

Salpeter, E. E. 1955, ApJ , 121, 161 [ Links ]

Sánchez, S. F. 2006a, Astronomische Nachrichten, 327, 850 [ Links ]

Sánchez, S. F., Aceituno, J., Thiele, U., Perez-Ramirez, D., & Alves, J. 2007a, PASP , 119, 1186 [ Links ]

Sánchez, S. F., Cardiel, N., Verheijen, M. A. W., Pedraz, S., & Covone, G. 2007b, MNRAS , 376, 125 [ Links ]

Sánchez, S. F., García-Lorenzo, B., Jahnke, K., Mediavilla, E., González-Serrano, J. I., Christensen, L., & Wisotzki, L. 2006b, New Astr. Rev., 49, 501 [ Links ]

Sánchez, S. F., Kennicutt, R. C., Gil de Paz, A., et al. 2012a, A&A , 538, A8 [ Links ]

Sánchez, S. F., Rosales-Ortega, F. F., Jungwiert, B., et al. 2013, A&A , 554, A58 [ Links ]

Sánchez, S. F., Rosales-Ortega, F. F., Kennicutt, R. C., Johnson, B. D., Diaz, A. I., Pasquali, A., & Hao, C. N. 2011, MNRAS , 410, 313 [ Links ]

Sánchez, S. F., Rosales-Ortega, F. F., Marino, R. A., et al. 2012b, A&A , 546, A2 [ Links ]

Sánchez-Blázquez, P., Ocvirk, P., Gibson, B. K., Perez, I., & Peletier, R. F. 2011, MNRAS , 415, 709 [ Links ]

Sánchez-Blázquez, P., Peletier, R. F., Jimenez-Vicente, J., Cardiel, N., Cenarro, A. J., Falcón-Barroso, J., Gorgas, J., Selam, S., & Vazdekis, A. 2006, MNRAS , 371, 703 [ Links ]

Sánchez-Blázquez, P., Rosales-Ortega, F., Diaz, A., & Sánchez, S. F. 2014a, MNRAS , 437, 1534 [ Links ]

Sánchez-Blázquez, P., Rosales-Ortega, F. F., Mendez-Abreu, J., et al. 2014b, A&A , 570, A6 [ Links ]

Sarzi, M., Falcón-Barroso, J., Davies, R. L., et al. 2006a, MNRAS , 366, 1151 [ Links ]

__________. 2006b, MNRAS, 366, 1151 [ Links ]

Serra, P. & Trager, S. C. 2007, MNRAS , 374, 769 [ Links ]

Stoklasová, I., Ferruit, P., Emsellem, E., Jungwiert, B., Pecontal, E., & Sánchez, S. F. 2009, A&A , 500, 1287 [ Links ]

Tinsley, B. M. 1980, FCPh, 5, 287 [ Links ]

Tojeiro, R., Heavens, A. F., Jimenez, R., & Panter, B. 2007, MNRAS , 381, 1252 [ Links ]

Tuffs, R. J., Popescu, C. C., Volk, H. J., Kylafis, N. D., & Dopita, M. A. 2004, A&A , 419, 821 [ Links ]

Vazdekis, A., Sánchez-Blazquez, P., Falcón-Barroso, J., Ce-narro, A. J., Beasley, M. A., Cardiel, N., Gorgas, J., & Peletier, R. F. 2010, MNRAS , 404, 1639 [ Links ]

Walcher, C. J., Boker, T., Charlot, S., Ho, L. C., Rix, H.-W., Rossa, J., Shields, J. C., & van der Marel, R. P. 2006, ApJ , 649, 692 [ Links ]

Walcher, C. J. , van der Marel, R. P. , McLaughlin, D., Rix, H.-W. , Boker, T. , Haring, N., Ho, L. C. , Sarzi, M., & Shields, J. C. 2005, ApJ , 618, 237 [ Links ]

Walcher, J., Groves, B., Budavari, T., & Dale, D. 2011, Ap&SS, 331,1 [ Links ]

Wilkinson, D. M., Maraston, C., Thomas, D., et al. 2015, MNRAS , 449, 328 [ Links ]

Worthey, G. 1994, ApJ S, 95, 107 [ Links ]

Received: May 25, 2015; Accepted: October 15, 2015

M. Cano-Díaz, J.J. González, C. López-Cobá, & S. F. Sánchez: Instituto de Astronomía, Universidad Nacional Autónoma de México, A.P. 70-264, 04510, México, D.F., México.

E. Perez: Instituto de Astrofísica de Andalucía (CSIC), Glorieta de la Astronomía s/n, Aptdo. 3004, E18080-Granada, Spain.

Y. Ascasibar & P. Sánchez-Blazquez: Departamento de Física Teorica, Universidad Autónoma de Madrid, 28049 Madrid, Spain.

F. F. Rosalez-Ortega: Instituto Nacional de Astrofísica, (Óptica y Electronica, Luis E. Erro 1, 72840 Tonantzintla, Puebla, Mexico.

A. Gil de Paz & R. A. Marino: CEI Campus Moncloa, UCM-UPM, Departamento de Astrofísica y CC. de la Atmosfera, Facultad de CC. Físicas, Universidad Complutense de Madrid, Avda. Complutense s/n, 28040 Madrid, Spain.

M. Molla: Departamento de Investigacion Basica, CIEMAT, Avda. Complutense 40 E-28040 Madrid, Spain.

A. R. Lopez-Sánchez: Australian Astronomical Observatory, PO BOX 296, Epping, NSW 1710, Australia. J. Barrera-Ballesteros: Instituto de Astrofísica de Canarias (IAC), E-38205 La Laguna, Tenerife, Spain.

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License