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 ². 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 ¹¹; (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 ² derived from the linear analysis of the underlying stellar population is stored. Then, the minimum X ² 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 ² changes up to X ² + 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}, ²⁰¹¹⁾. The basic steps of the fitting algorithm, spectrum by spectrum, are the following:

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

C_i = O_i - S_i.

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

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-³, 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 (λ _k,σ_k) 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

GF_i = G_i - E_i.

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 ² 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.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)}:

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 results¹². 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 ∞ Age_{j ^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 SSP_age,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 

¹∆v_sys 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 ². 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.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.