© ESO, 2014

1. Introduction

Large optical spectroscopic surveys have dominated the panorama of extra-galactic astrophysics over the last decade. Stellar population analysis of these data has produced information on stellar masses, mean stellar ages, abundances, dust content, and diagnostics of their star formation histories (SFHs), whose analysis per se and in relation to other properties (for instance, morphology or environment) has advanced the understanding of the physics of galaxies and their evolution through cosmic time. Yet, all this information is derived from spatially unresolved spectra, and thus reflects global averages over properties that are known to vary accross the face of galaxies. Integral field spectroscopy (IFS) is a natural next step to move from global properties to the internal physics of galaxies. Dissecting the different structural components will probably lead to a much more complete picture of galaxies, and to a clearer understanding of the results previously obtained in non-spatially resolved surveys. IFS has been around for over two decades, but has only recently been promoted to the scale of modern galaxy surveys. CALIFA, the Calar Alto Integral Field Area survey (Sánchez et al. 2012; Husemann et al. 2013), is a pioneer in this blooming field.

To harvest the full potential of IFS data one needs to adapt techniques developed for spatially unresolved galaxy spectra, folding in the information on the spatial structure. In a companion article (Cid Fernandes et al. 2013, hereafter Paper I) we have presented our method for extracting stellar population properties out of CALIFA datacubes by means of the starlight spectral-fitting software (Cid Fernandes et al. 2005). The packages developed to pre-process and partition the reduced datacubes (QBICK) and to post-process the results of the spectral synthesis (PyCASSO, the Python CALIFA Starlight Synthesis Organizer) were described in detail. The products of this work flow were examplified using data for the nearby spiral NGC 2916 (CALIFA 277). Diagnostics such as 2D maps of physical properties (mean ages and metallicities, mass surface densities, star formation rates, etc.), 1D averages in the temporal and spatial dimensions, projections of the stellar light and mass growth onto radius-age diagrams, and cuts through the (x,y,age) cubes were presented to illustrate the potential of the combination of IFS data with spectral synthesis as a tool for studying galaxy evolution in time and space simultaneously.

Uncertainties in this whole process stem from (i) noise and calibration of the data; (ii) limitations of our spectral synthesis method; and (iii) evolutionary synthesis models, the key ingredient in the mapping of observed spectra to astrophysical information. These different sources of uncertainty, each one comprising a field of work by itself, intertwine and propagate in complex ways through our pipelines. Uncertainties in spectral synthesis products have been previously explored in the literature, such as the noise-effect simulations in Cid Fernandes et al. (2005), the model-impact study of Panter el al. (2007), and the suite of tests presented by Tojeiro et al. (2007), but never in the context of IFS.

This article complements Paper I by addressing the effects of these uncertainties. Our central goal is to provide quantitative and qualitative guidelines to allow an assessment of the robustness of results produced with the methods and tools explained in Paper I, which are being used in a series of publications by our collaboration. Pérez et al. (2013), for instance, analyzed the spatially resolved mass assembly history of 105 CALIFA galaxies, finding clear evidence for an inside-out growth for massive galaxies, and signs of outside-in growth for the less massive ones (see González Delgado et al. 2012 and 2013 for other results). Despite the focus on starlight and CALIFA, the results reported here are also useful in the broader context of stellar population analysis of IFS data; for instance, experimental datacubes recently obtained by MaNGA run through our QBICK and PyCASSO pipelines without any modification.

The first part of this paper (Sects. 2 and 3) presents extensive simulations that explore the effects of uncertainties related to the data and spectral synthesis method. Random noise and shape-changing perturbations are added to both observed and synthetic spectra of all zones of the same galaxy analyzed in Paper I. These simulations are then used to evaluate uncertainties in global properties (mean ages, masses, etc.) and higher-order products such as SFHs. This is done both at the level of individual spectra and for spatial regions cut out from the datacube.

The second part (Sect. 4) is dedicated to experiments carried out to evaluate how our results are affected by different choices of evolutionary synthesis models for the spectra of simple stellar populations (SSP). This study uses ~100 thousand spectra from 107 galaxies observed by CALIFA, all of which were processed with three sets of SSP models used in the current literature. An inspection of the quality of the spectral fits and a comparison with results obtained for SDSS galaxies is also presented.

Section 5 summarizes our main results. Throughout this paper the notation and definitions established in Paper I are followed.

2. Uncertainties associated to the data and the method: description of the simulations

A drawback of starlight is that it provides no error estimate on its output. The most straightforward (albeit computationally expensive) way to do so is to run the code several times perturbing the input data according to some realistic prescription of the errors involved. This section describes a set of extensive simulations taylored to match the characteristics of the CALIFA data. The simulations are designed to address the effects of

• 1.

  random noise,

• 2.

  shape-related calibration uncertainties, and

• 3.

  multiplicity of solutions,

upon the physical properties derived with starlight. Throughout these experiments we keept the set of SSPs fixed at the same set as used in Paper I, which is denoted “base GM” in Sect. 4, where the effects of this choice are investigated. This base comprises N_⋆ = 156 elements with 39 ages between t = 1 Myr and 14 Gyr, and 4 metallicities from Z = 0.2 to 1.5 solar. The simulations are based on the data and starlight fits of 1638 zone spectra¹ of CALIFA 277 studied in Paper I. Given the variety of properties found accross this galaxy and the fact that all other datacubes in CALIFA were processed in the same way, the uncertainty estimates presented here are representative of CALIFA data in general.

Each spectrum was perturbed ten times for each of the R, C, and E error prescriptions discussed below and refitted. This data set is then used to evaluate uncertainties in the derived properties, both at the level of individual zones and for spatially averaged properties like those shown in Figs. 10 and 11 of Paper I. Observed as well as synthetic spectra were perturbed in the simulations. The first represent real-life conditions, while the latter are essentially a theoretical exercise which explores the sensitivity of the starlight results to different types of perturbations.

2.1. Types of errors and simulations

The first type of error considered is that related to the variations on the derived properties induced by random fluctuations of the observed fluxes within the noise level. To study this, each zone spectrum was perturbed ten times with Gaussian noise with a λ-dependent amplitude given by the error spectrum ϵ_λ (hereafter OR1 runs, with O standing for observed, R for random, and 1 for the amplitude of the perturbation). Experiments were also conducted with pertubations of 2 and 3 × ϵ_λ (OR2 and OR3 runs, respectively) to emulate what would be obtained under conditions considerably worse than those in CALIFA².

A second set of simulations was carried to address the problem of the overall continuum shape calibration. starlight relies on flux-calibrated input spectra, so it is relevant to access how much its solutions vary for spectra whose shapes vary within expected uncertainty ranges. Husseman et al. (2013) found that CALIFA synthetic g − r colors typically differ by 0.05 mag from those obtained by SDSS photometry. We used this value to gauge the effects of continumm calibration uncertainty in each individual specrum of CALIFA 277. For simplicity, but with no loss of generality, these simulations keep the flux at λ_r = 6231 Å stable, such that the 0.05 mag variation in g − r translates into a 0.02 dex uncertainty in the flux at λ_g = 4770 Å. To emulate the effects of such variations we produced ten versions of each spectrum by adding A_g × (λ_r − λ)/(λ_r − λ_g) to the original log F_λ, with A_g = Δlog F_λg sampled from a Gaussian distribution of zero mean and dispersion σ(log F_λ4470) = 0.02 (OC002 runs, where C stands for color). This produces both redder and bluer versions of the original spectrum. Experiments were also carried out for σ(log F_λ4470) = 0.01, and 0.04 dex (OC001 and OC004, respectively). We take OC002 runs as default for CALIFA, even though the spaxel-to-spaxel variations in the continuum calibration for a single galaxy are probably smaller than the g − r dispersion quoted by Husemann et al. (2013), which was derived from the 100 galaxies in the CALIFA DR1 sample.

Finally, fits were also conducted without any perturbation at all (OE0 runs), the motivation for which is to quantify a qualitatively different source of uncertainty embedded in all other simulations. starlight provides a single best-fit set of parameters out of typically many millions tried out during its likelihood-guided sampling of the parameter space. This single solution, however, is not meant to be mathematically unique. In fact, given the pseudo-random nature of its Markov chains, starlight solutions depend on the input seed for the random-number generator. For large samples, and each CALIFA datacube is a large sample in itself, these variations have no relevant effect on the overall results of the analysis. For individual spectra, however, the multiplicity of solutions leads to uncertainties in the derived properties. The OE0 runs quantify these uncertainties, which ultimately reflect an intrinsic limitation of our spectral synthesis method.

2.2. Simulations based on synthetic spectra

The OE0, OR*, and OC* runs work upon the original observed spectra (hence the prefix O). A paralell set of simulations was carried out applying these same perturbation recipes to the synthetic spectra obtained in our starlight fits of CALIFA 277. We label these runs as SE0, SR1, SR2, SR3, SC001, SC002, and SC004 (S for synthetic).

The S runs use exactly the same sampling, masks, flags, and ϵ_λ errors as in the O runs, thus following the pattern of CALIFA data in a realistic fashion. At the same time, these are idealistic simulations in the sense that the perturbations act upon perfect input spectra instead of on observed ones. This theoretical exercise therefore overlooks inadequacies in the SSP models (including problems such as potential flux calibration problems in the libraries, limited coverage of the stellar parameter space, missing or incorrectly modeled evolutionary phases) and other assumptions involved in the modeling (namely, extinction and kinematics).

The purpose of these runs is to evaluate to which extent starlight recovers known input properties in the presence of perturbations. This kind of input versus output comparison is more commonly carried out by feeding the spectral fitting code with data generated with ad hoc parametric descriptions of the SFH (e.g., the so-called tau models), while our theroretical SFHs are inspired in those derived from actual observations. A valid criticism of these simulations is that the same base is used to generate and fit spectra. The alternative would be to generate models with some other set of t′s and Z′, but since our base is so large (39 ages and all four available metallicities), this would make very little difference in practice.

2.3. Notation

For clarity and quick reference, we summarize the notation and intended purpose of the simulations described above.

R runs (OR1, OR2, OR3, SR1, SR2, and SR3) explore the effect of noise by adding Gaussian perturbations scaled to have amplitudes 1, 2, or 3 times the nominal error ϵ_λ. C runs perturb the input spectra in their global shape by multiplying the original F_λ by a power law in λ, scaled such that the flux at λ_r remains fixed while that at λ_g differs by (on average over all realizations) 0.01 (OC001 and SC001), 0.02 (OC002 and SC002) and 0.04 (OC004 and SC004) dex with respect to the input flux. (In a log F_λ × λ plot, these shape-pertubations correspond to straight lines anchored at λ_r and different slopes.) These runs are meant to emulate the effect of a flux calibration uncertainty in the data. E0 runs do nothing to the input spectra. They only differ in the random number seed feed to starlight. The variance among these fits reflects the multiplicity of SSP combinations which lead to a similar total spectrum.

In all cases, a prefix O indicates that the perturbations act upon the observed spectra, while an S prefix means that synthetic spectra are perturbed. Table 1 summarizes the simulation-related nomenclature used thoughout the paper.

3. Results of the simulations

The simulations described in Sect. 2 provide the material to evaluate uncertainties associated with the data and spectral synthesis method. We split the presentation of this material into an assessement of the errors in global properties (Sect. 3.1), radial profiles (Sect. 3.2), and SFHs (Sect. 3.3). We conclude this part of our investigation with a summary of the results of the simulations (Sect. 3.4).

3.1. Uncertainties in global properties

The main scalar properties derived from the spectral synthesis are mean ages and metallicities (weighted by light or mass), extinction and stellar masses³. These first order products of the analysis convey enough information to draw a quite detailed sketch of the 2D distribution of stellar population properties in a galaxy. Before incorporating the spatial arrangement of the spectra, however, it is useful to assess uncertainties in the individual spectra as if they came from different galaxies.

Figure 1 shows histograms of the Δ ≡ simulation − original value (i.e., that obtained from the unperturbed spectrum) of each of these quantities, and for six different sets of simulations: E0, R1, and C002 runs for the O and S versions. Statistics for all runs are listed in Table 2.

3.1.1. E runs: noise-free simulations

The OE0 and SE0 panels show that variations due to the random number seed fed to starlight have a tiny effect on the derived properties, with no bias () and one σ uncertainty all in the second decimal place.

Intriguingly, SE0 runs have somewhat higher dispersions than OE0 runs. The mathematical reason for this behavior is subtle, but simple to grasp with the concepts developed by Pelat (1997, 1998). By construction, SE0 spectra lie within the synthetic domain (SD), i.e., the space generated by convex combinations of the base spectra⁴. A point inside this space can be exactly fitted by at least one set of parameters. Even if the solution is mathematically unique, a whole range of other solutions produce nearly identical spectra. Because of the limitations of evolutionary synthesis models, observed spectra cannot possibly be fitted exactly even in the ideal limit of perfect calibration and S/N = ∞ signal-to-noise ratio, so OE0 spectra are inevitably outside the SD. Fits of OE0 spectra therefore have access to a smaller subspace of acceptable solutions than fits of SE0 spectra, which results in the lower dispersions observed in Fig. 1. This same interpretation applies to other runs. Even when the synthetic spectra are perturbed (for instance, in the SR1 simulations), they remain closer to the SD than the corresponding O spectra, thus reaching a larger region of the solutions space.

3.1.2. R runs: random-noise simulations

Random noise at the nominal level (OR1 and SR1 runs, middle rows in Fig. 1) lead to Δ-distributions with σ around 0.1 dex in ages, metallicities and masses, and 0.05 mag in A_V. As expected (and as also seen in the E0 runs), mass-weighted ages and metallicities have broader Δ distributions than their luminosity-weighted counterparts. In any case, these dispersions are all relatively small. Indeed, they are smaller than those induced by the use of different SSP models (Sect. 4, Fig. 11).

Table 2 reports results for simulations with noise levels 2 and 3 times worse than ϵ_λ. For the OR2 runs, for instance, the uncertainties in the properties ploted in Fig. 1 increase by roughly a factor of 2 with respect to those found in the OR1 simulations. The Δ distributions (not plotted) remain symmetric, without significant biases. We recall, however, that there is no reason to believe that the ϵ_λ values produced by the CALIFA reduction pipeline are underestimated. On the contrary, Husemann et al. (2013) showed that, if anything, they are slightly overestimated.

3.1.3. C runs: shape perturbations

The bottom panels in Fig. 1 show the results for simulations where the continuum shape was perturbed to produce (on average) 0.02 dex variations in the flux at λ_g for a fixed flux at λ_r. The one σ variations in ages, metallicities, and masses for OC002 runs are larger than those obtained with the OR1 runs, but not as large as in OR2 (Table 2). Similar conclusions apply to the corresponding S runs. Predictably, A_V is more severely affected by these global shape perturbations, with σ(ΔA_V) ~ 2.5 times higher than those in the OR1 runs.

Unlike in other simulations, the Δ distributions are skewed, particularly those for A_V and ⟨log t⟩_L. This is a by-product of the A_V ≥ 0 limit imposed upon the fits. 582 of the 1638 zones of CALIFA 277 were fitted with A_V = 0 (see Fig. 4 in Paper I). About half of the shape-perturbed versions of these spectra are bluer than the original ones. Since starlight was not allowed to decrease A_V to fit them, it gets stuck at A_V = 0, compensating the extra blueness with younger populations, whicg in turn produces the tail toward negative Δ ⟨log t⟩_L seen in the bottom rows of Fig. 1. This is also the reason why there are more instances of ΔA_V > 0 than < 0, producing the positively skewed ΔA_V distributions.

Increasing the shape errors by a factor of 2 (C004 runs) leads to uncertainties roughly twice as large, biases of the order of − 0.1 dex in ages, and + 0.15 mag in A_V (Table 2), as well as larger skewness of the Δ distributions. Conversely, C001 runs result in uncertainties ~2 times smaller, with no bias nor much skewness.

3.1.4. Covariances

As usual in multivariate problems, variations in one property correlate with variations in others. In the context of spectral synthesis, such covariances are known as age-metallicity-extinction degeneracies: A_V affects colors while t and Z affect both absorption lines and colors. Previous work (e.g., Cid Fernandes et al. 2005; Sánchez-Blázquez et al. 2011) showed that λ-by-λ fits help in reducing these degeneracies, but do not eliminate them.

Figure 2 shows how the main scalar properies co-vary in our simulations. S runs were used in this plot-O runs produce similar plots, with slightly more compact contours. The first two columns (SR1 and SR2 runs) reveal the classical degeneracies of population synthesis, with mean ages increasing as mean metallicities decrease, and extinction variations also anti-correlated with variations in mean age. The bottom row shows how stellar masses increase as ages become older, and vice-versa.

In the R-runs, perturbations do not change the overall shape of the input spectrum. Hence, if a starlight solution turns out to be somewhat older than the original one, its extinction and/or metallicity will decrease to compensate for the redder colors and stronger absorption lines of the older stars. C-runs, on the other hand, do change colors-that is, in fact, exactly what they are designed to do. This slightly changes the nature of “degeneracies”.

If the shape perturbation followed the same functional form of the extinction law, it would suffice to change A_V, leaving the population mixture essentially unchanged, varying ts and Zs at the level they would naturally vary without color perturbations (i.e, due to noise and the multiplicity of solutions). To first order, this is indeed what happens in the C-runs. Compare, for instance, the (SR1, SR2) and (SC001, SC002) pairs of ΔA_V × Δ ⟨log t⟩_L panels in Fig. 2. In both pairs of plots the corresponding type of error changes by a factor of 2 in amplitude, yet one sees that the ΔA_V contours grow more from SC001 to SC002 than from SR1 to SR2. Furthermore, the ΔA_V contours are nearly vertically oriented in the C-runs, while in R-runs the anti-correlation with Δ ⟨log t⟩_L is clear (the age-extinction degeneracy). We thus confirm the expectation that, to first order, continuum shape uncertainties primarily affect the A_V estimates.

There are still two second-order features to explain about the ΔA_V × Δ ⟨log t⟩_L C-run maps in Fig. 2. First, the inner contours have a slightly positive slope, opposite to that expected from the age-extincion degeneracy (clearly seen in the SR1 and SR2 panels). This behaviour arises because our shape-perturbations (where Δlog F_λ changes linearly with λ) cannot be fully absorbed by the reddening curve used in the fits, because of the different functional forms. Since changing A_V alone is not quite enough to produce optimal spectral fits, starlight searches for other ways of making the synthetic spectrum redder or bluer (depending on the sign of the perturbation), and skewing the age mixture is the most effective one. In any case, the change in mean age is small: for SC002 runs (Table 2). Secondly, for SC002 runs, the outer contours, which contain about 20% of the points, are limited to ΔA_V ≤ 0 whenever simulated mean ages are more than 0.1 dex younger than the original ones (Δ ⟨log t⟩_L < − 0.1). Spectra along this odd-looking ridge were originally fitted with A_V = 0, and then perturbed to bluer colors. Since our fits do not allow A_V < 0, starlight fitted them with the same A_V (hence ΔA_V = 0), but younger ages, leading to the pattern observed in this panel.

3.1.5. Checking the effects of imposing A_V ≥ 0

As seen just above, the effects of a hypothetical systematic redness associated to spectroscopic calibration problems are almost entirely absorbed by biasing A_V upward, while an excessive blueness of the same nature stumbles upon the A_V ≥ 0 constraint, leading to biases in the age mixture. An artificial way of fixing this assymetry is to allow for A_V < 0 in the spectral fits. A trivial statistical justification to allow for this unphysical possibility is that when A_V is truly = 0, unbiased estimates should oscillate around 0, including both negative and positive values. We also note that completely independent studies (different methods, different data) also allow for or deduce negative extinctions (e.g. Kauffmann et al. 2003, their Fig. 11).

An extra set of OE0-like simulations, but imposing A_V > −1, was run to investigate this issue. In practice what happens is that only zones that were originally fitted with A_V ~ 0 are systematicaly refitted with negative A_V. The 591 zones that result in A_V < 0 in these new fits include all the 582 for which we had originally found A_V = 0, and the remaining 9 had A_V very close to 0. The mean new A_V of these 591 zones was − 0.12 mag, with a standard deviation of 0.08 mag. The lowest value reached was A_V = −0.4 mag, but in only 16% of the cases A_V went below − 0.2 mag. The variations in other global properties (for these same 591 zones) were , , , , and dex. Mean ages are the most affected property, but to an acceptably small degree, with biases of the same order as the one-σ variations produced by random and shape-related uncertainties at the nominal level (R1 and C002 simulations).

Zones originally fitted with A_V > 0 did not show systematic variations in their global properties in these extra runs.

3.1.6. The effect of S/N

The Δ histograms in Fig. 1 and the statistics in Table 2 represent the uncertainties of the CALIFA 277 data set as a whole, including zones of different S/N. The general expectation is that a lower S/N leads to larger uncertainties, such that the tails of the Δ distributions in Fig. 1 are dominated by the noisier spectra, while the cores (small Δ values) correspond to a higher S/N.

To verify this, Fig. 3 shows the Δ ⟨log t⟩_L variations against the S/N of the original spectra. The lines represent the mean ± 1σ range. Uncertainties do decrease with increasing S/N in the SE0, SR1, and SR2 runs, but not in the C-runs. This is in fact expected because the shape perturbations have the same amplitude for all spectra, independent of their S/N. The right-most panel illustrates the bias toward younger ages that results from the A_V ≥ 0 barrier.

We note in passing that, as hinted in Fig. 3, the distribution of S/N values in CALIFA 277 is skewed toward the target minimum S/N of 20 stipulated in the spatial-binning scheme explained in Paper I. This same threshold was applied to all other datacubes, so the uncertainty estimates presented here should be representative of the CALIFA data in general, even though they are based on a single galaxy.

3.2. Exploring the statistics in IFS data: uncertainties in radial profiles

So far our error analysis has treated the datacube as a collection of unrelated spectra. Spectral synthesis does not care where a spectrum comes from, so the positional information provided by IFS has no effect on the uncertainties for individual zones. However, averaging results obtained for different zones brings obvious statistical benefits, reducing uncertainties in properties derived from individual spectra.

Figure 10 in Paper I shows azimuthally averaged profiles for a series of global properties. Despite the scatter at any given radius, gradients in age and metallicity were clearly detected in CALIFA 277. Figure 4 in this paper repeats that plot, but now showing the radial profiles derived from the simulated spectra, so that we can evaluate the uncertainties in azymuthal averages. The original profiles are plotted in solid black, while colored shades represent the ranges spanned by OR1 (blue) and C002 (red) runs. Each panel also shows an (arbitrarily shifted) simulation minus the original difference profile at its bottom.

As expected, random noise does not change the radial profiles significantly. A small bias of about 0.02 mag is seen in for radii beyond ~1 half light radius (HLR). This is again a side effect of the A_V ≥ 0 constraint: when a spectrum originally fitted with A_V ~ 0 is perturbed, some will be refitted with A_V > 0, but none will be refitted with A_V < 0. The resulting bias propagates to minor offsets in other properties, like slightly younger ages (panels e and f). Overall, however, these are all tiny effects.

Color perturbations (red-shaded regions) produce larger effects. The main effect is on the profile (panel c), where the typical offset matches the expected mag reported in Table 2. Leaving aside this bias, the variations are much smaller than those in individuals zones. For instance, considering all radial bins and all 10 × 1638 versions of the OC002 spectra, the rms in is only 0.05 mag, while individual zones are subjected to uncertainties of σ(A_V) = 0.16 mag (Fig. 1).

Larger than typical variations occur in the nucleus for the trivial reason that it comprises a single spaxel, so that no actual statistics is performed. More robust results for the nucleus require extending its definition to include more spaxels. For instance, Pérez et al. (2013) define the nucleus as the inner 0.1 HLR, which in the case of CALIFA 277 comprises the central-most nine spaxels. Uncertainties also increase in the outer regions (R > 2 HLR), induced by the lower S/N of the spectra and by the decrease in the number of zones compared with intermediate radii.

Figure 4 confirms the expected gain in dealing with averaged results. The simulations show that the uncertainties in the radial profiles are just too small to question the reality of the age and metallicity gradients in this galaxy. This conclusion holds even for our most pessimistic simulations (C004 and R3). This gain comes at the expense of a loss of spatial information–compressing the x,y plane to a single dimension (R) in the example presented here. Other forms of spatial averaging can be envisaged, such as isolating structural components (bulge, disk, bar, arms, inter-arm regions). From the results above, it is easy to see that as long as enough zones are involved in the statistics, any averaging scheme should lead to small uncertainties.

3.3. Uncertainties in star formation histories

starlight outputs SFHs as population vectors: x and μ. The first gives the fractional contribution of each base element to the total flux at the chosen reference wavelength (5635 Å in our case), while the latter gives the corresponding mass fractions. The 39 ages in our base⁵ sample the 1 Myr to 14 Gyr interval with a mean log t separation of just 0.1 dex. Individually, the strength of each component is not reliable (see Fig. 1 in Cid Fernandes et al. 2004 for an example), given the many ways that spectrally similar components can rearrange their x_j values to produce a similar total spectrum. Yet, their combination in indices such as ⟨log t⟩_L produces robust quantities, as shown by the simulations.

None of the properties studied in the previous sections relies heavily on the details of the spectral decomposition performed by starlight. Mean ages, for instance, are defined as the first moment of the age distributions encoded in x and μ. By construction, however, SFHs do rely on higher-order descriptions of the age distribution. This section aims to evaluate the uncertainties in our SFHs. Before doing so, we digress on the actual meaning of the “H” in SFH.

There is no standard semantics for the meaning of “history” in spectral synthesis studies. On the contrary, the same term is used in very different ways in the literature, always qualitatively referring to some measure of the how star formation unfolded over time, but quantitatively varying from estimates of mean ages to measures of the “recent” star formation (where recent can mean anything between ~10 Myr and a few Gyr, depending on the tracers), and full t-by-t descriptions with varying age resolution. In addition to different timescales, different metrics can be used to quantify the SFH, such as luminosity or mass associated to stars in a given age range, in differential or cumulative forms, etc. This diversity is particularly natural in the case of non-parametric codes such as starlight. Because no pre-established functional forms for SFR(t) are imposed, one has plenty of freedom to design different indices and functions for characterizing the SFH. Paper I is itself an example of this diversity. We have presented there a series of SFH diagnostics obtained by manipulating the x and μ vectors in different ways, always with the underlying philosophy of compressing the overly-detailed output of starlight into coarser but more robust descriptions of the SFH.

In short, there are numerous alternative ways to describe an SFH. In what follows we evaluate the uncertainties in a paper I-inspired selection of SFH-descriptions. Following the sequence of our analysis of global properties, we first present results for uncertainties in individual zones (Sect. 3.3.1) and then for spatially averaged SFHs (Sect. 3.3.2).

3.3.1. SFHs for individual zones

One way to characterize SFHs is to compress the population vectors into age ranges, an approach that dates back to the early days of spectral synthesis (Bica 1988). In the context of starlight, Cid Fernandes et al. (2004 and 2005) showed that a robust summary of the SFH can be obtained by binning x into young (Y), intermediate (I) and old (O) groups. This x → (x_Y,x_I,x_O) compression scheme was employed to produce the maps shown in Figs. 8 and 9 of Paper I, where these age groups were separated by frontiers placed at t = 0.14 and 1.4 Gyr.

We have used our simulations to assess the uncertainties in these fractions. The top three rows in Table 3 list the O-run statistics of Δx_Y, Δx_I and Δx_O, where, as before, Δ ≡ simulated − original value. The table shows that random noise at the one-σ level (R1 runs) induces uncertainties in x_Y, x_I and x_O of ~3, 9 and 9% respectively⁶. Larger uncertainties (~15%) in x_I and x_O apply for shape-changing perturbations at the level expected for CALIFA (OC002 runs). Part of this increase is due to the aformentioned A_V ≥ 0 effect, which forces low A_V zones that are perturbed to bluer spectra to become younger to compensate for the inability to reduce A_V below 0. This is also reflected in the biases, with %. At any rate, comparing these uncertainties with %, the average values over all zones in CALIFA 277, one concludes that this compact description of the SFH is reasonably robust.

Figure 5 show how the variations in these light fractions correlate with each other. Δx_O anti-correlates very strongly with Δx_I, which is not suprising given the x_Y + x_I + x_O = 100% constraint and the fact that these are the age groups that carry most of the light in CALIFA 277. The Y group, on the other hand, is roughly independent of the others.

Figure 6 shows OR1-based results for a finer-graded description using five age groups defined by log t   [yr] = 6–7.5, 7.5–8.5, 8.5–9.3, 9.3–9.8, and 9.8–10.2. The anti-correlations among adjacent age groups are visible along the diagonal. The σ_Δx values are 3, 6, 9, 13, and 10% (Table 3) for mean fractions of , 3, 24, 34 and 26%, respectively. This more detailed description shortens the difference between σ_Δx values and the corresponding . The uncertainties are therefore larger, as expected, but these higher-resolution SFHs are still useful, especially if used in a statistical way.

As the name implies, starlight fits light, so x-based SFH descriptions are the natural way to analyze its results. Yet, it is the stellar mass which carries the astrophysical information on the galaxy history. One can obviously group μ elements just as we did for x, but while light fractions in different age groups can be directly intercompared, mass fractions cannot. With a factor of ~1000 decrease in light-to-mass ratio from the youngest to the oldest elements in the base, wildly different mass fractions, say, μ_Y ~ 1% and μ_O ~ 99%, can be equally relevant for the fits.

To evaluate uncertainties in the masses associated to different age groups it is better to use the masses themselves instead of mass fractions. Table 4 reports the uncertainties in log M′, where the prime denotes the initial stellar mass, uncorrected for the mass returned to the interstellar medium by stellar evolution. The table reports the results of the simulations for both the (Y, I, O) and the five age groups descriptions of the SFH discussed above. Again, the fact that we are now exploring higher moments of the age distribution implies larger uncertainties than for global properties. A column-by-column comparison with the Δlog M′ statistics in Table 2 confirms this. For instance, while the total masses have uncertainties of just 0.08 dex in OR1 runs, the masses in age-groups are typically three times more uncertain, an unavoidable but still reasonable price to pay for a time-resolved description of results of the spectral synthesis.

Additional subdivisions lead to a quasi-continuous regime that is better handled in ways other than age-binning, such as cumulating x or μ. Figure 7 illustrates this. The left panels show the cumulative light fraction, shown after a cosmetic smoothing by a FWHM = 0.2 dex Gaussian in log t. Six individual zones, including the nucleus (top) and the last one (bottom), were chosen to exemplify this alternative description of SFHs and its uncertainties. The thick solid line traces the original starlight solution, while thin gray lines trace the solutions obtained for each of the ten OR1 perturbations, whose median is drawn as a dashed line. The red-shaded bands mark the range spanned by the six SFHs closer to the median at each t. The right panels show the corresponding mass assembly histories, all scaled to 100% at t = 0. The six zones in Fig. 7 cover different distances from the nucleus, S/N, and SFHs, thus providing representative examples of starlight results for CALIFA data.

Except for the nucleus, where the perturbed spectra and the associated SFHs are nearly indistinguishable from the original ones because of the very high S/N, there is a clear spread in cumulative light and mass distributions. The dispersion is higher for ages older than a few Gyrs, which contain a substantial fraction of the light and most of the mass. The mild spectral evolution in this age range (Conroy 2013) explains the larger uncertainties, because components can be interchanged with relatively little variation in the total spectrum. Gauging the uncertainties in the cumulative SFHs by the width of the red band, we find it to be of the order of 10–20%, comparable to those reported in Tables 3 and 4 for age-grouped quantities.

Despite uncertainties, the example SFHs in Fig. 7 reveal a gradual shift toward younger ages for larger radii. This trend is the higher-order description of the clear mean age radial gradients previously shown in panels e and f of Fig. 4, where results for individual zones at same R were averaged.

3.3.2. Spatially averaged SFHs

In analogy with Sect. 3.2, uncertainties in the SFHs of individual zone spectra should be greatly reduced when averaging over spatial regions comprising many zones. This is confirmed in Fig. 8, which shows cumulative light and mass distributions (in both absolute and relative units) for four spatial regions: the nucleus (gray), R < 0.5 (green), R < 1 (red), and R > 1 HLR (blue). The smoothing here is by FWHM = 0.5 dex in log t for compatibility with Fig. 11 of Paper I. OR1 runs are shown in the left panels, OC002 in the right panels.

The statistical damping of the uncertainties is so strong that one can hardly notice that, as in Fig. 7, shaded bands are used to show the vertical ranges spanned by the simulations. The bands are so narrow that they look like thick lines. The exception is the nucleus (gray), where the uncertainties are visible for the simple reason that unlike for other regions, it does not involve any actual spatial averaging because it comprises a single spaxel.

Close inspection of the right panels reveals that the shape-changing simulations tend to have SFHs slightly skewed toward younger ages than SFHs inferred from spectral fits of the original data (marked by dashed lines). This is again the effect of some perturbations making low A_V zones too blue to be fitted with A_V ≥ 0, so starlight compensates it with lower ages.

In conclusion, spatial averaging effectively eliminates the dispersions seen in the SFHs of individual zones. This remains valid even for simulations where the level of perturbations is higher than that expected for CALIFA spectra.

3.4. Summary of the simulations

The overall message from these experiments is that degeneracies inherent to stellar populations and uncertainties in the data have a limited impact upon our analysis. For the level of uncertainties expected in CALIFA, here represented by the R1 and C002 runs, the simulations yield acceptably small uncertainties in global properties. More precisely, mean ages and metallicities have one sigma variations of ~0.1 dex when weighted in luminosity and ~0.15 weighting by mass, while stellar masses are good to ~0.1 dex (round numbers), with negligible biases. Uncertainties in A_V are ~0.07 mag if only random noise is considered, increasing to 0.16 mag when shape-related uncertainties are accounted for. In the latter case, a bias of mag appears as a side-effect of forbidding negative A_V. Because they are a higher-order product of the spectral synthesis, SFHs are naturally less strongly constrained than global properties, but the age-grouped and cumulative descriptions explored above were found to be relatively stable against random and spectral-shape perturbations.

It is a consensus in the field of spectral synthesis that the results are best used in a statistical sense (e.g., Panter et al. 2006), and the present study endorses this view. First, averaging reduces uncertainties, as confirmed by our experiments with radial profiles (Fig. 4) and spatially averaged SFHs (Fig. 8), which showed that the results are stable against the types of uncertainties studied here. Second, and more importantly, large samples allow for comparative studies, which can reveal astrophysically valuable trends while mitigating worries about data or method-related systematics affecting the absolute values. This point has been abundantly explored in spectral synthesis studies of SDSS galaxies, where unequivocal relations between, for instance, age, mass, and metallicity, have been identified in spite of the uncertainties in the properties derived for individual galaxies. The analogy with IFS data is clear. In IFS each galaxy is a large sample by itself, and the very goal of IFS is to compare properties in one place with those somewhere else.

This concludes our study of data- and method-related uncertainties. This whole study was based on a single set of SSP models, a fundamental prior in the analysis. We now turn to the uncertainties resulting from this choice.

4. Uncertainties associated to evolutionary synthesis models: an empirical evaluation

SSP spectra are the key astrophysical ingredient in any spectral synthesis analysis, the dictionary that translates stellar properties into spectroscopic information, which can then be compared with data through different methods to extract physical properties and SFHs. The computation of model SSP spectra combines an assumed initial mass function (IMF) with stellar evolution theory (isochrones) and a library to provide the spectra of stars with temperature, gravity, and chemical composition matching the values dictated by the evolutionary tracks. In practice, this algebric exercise stumbles upon a series of problems that directly affect the predicted SSP spectra.

First, despite its (well-deserved) status as a founding stone of modern astrophysics, stellar evolution is still an open business. Rotation in massive stars (Levesque et al. 2012), rapid phases during the post-main-sequence evolution of intermediate-mass stars (Maraston 2005), and binarity (Li et al. 2013) are examples of phenomena that affect the radiative output of stellar populations, but are still at the forefront of research, and thus prone to uncertainties. A practical example of how incomplete treatment of evolutionary phases affects spectral synthesis analysis is discussed by Koleva et al. (2008), Cid Fernandes & González Delgado (2010), and Ocvirk (2010). Using different codes, these studies found that spectral synthesis of old Milky Way and LMC globular clusters suggests the presence of very young stellar populations accounting for some 10–20% of the optical light, an astrophysically absurd result that also happens in some old, “red and dead” galaxies (Cid Fernandes et al. 2011). These “fake young bursts” (Ocvirk 2010) are in fact by-products of inadequate modeling of the horizontal branch (and/or blue straglers, Cenarro et al. 2008). Seeking for a best match, the spectral synthesis uses young hot stars to replace the old hot ones missing in the models.

Second, spectral libraries, whether empirical or theoretical, have their own limitations (Martins & Coelho 2007). The widely used Bruzual & Charlot (2003; BC03) models, for instance, are based on the STELIB library (Le Borgne et al. 2003), a major advance in its days, but limited to 249 stars, some of which have severe spectral gaps. Its relatively poor coverage of the stellar parameter space propagates to the predicted SSP spectra, with collateral effects on the results of a spectral synthesis analysis. Koleva et al. (2008), for instance, noticed that because of the lack of truly very metal-rich stars in the library, BC03 spectra for 2.5 Z_⊙ SSPs behave like ~Z_⊙ models of an older age. In short, this too is an evolving field, undergoing constant revisions.

In this unstable scenario, it is relevant to evaluate to which extent the results of our spectral synthesis analysis depend on the choice of SSP models. To this end, we performed spectral fits using three bases extracted from different sets of evolutionary synthesis models, listed in Table 5. We stress that this is by no means an exaustive exploration of this matter, but it suits our goal of estimating uncertainties related to this choice. Section 4.1 describes the sets of models to be compared. The CALIFA data used in this study are described in Sect. 4.2. Comparisons of the properties obtained with different sets of SSPs are presented in Sect. 4.3. Fit-quality indicators and spectral residuals are discussed in Sect. 4.4, while Sect. 4.5 compares residuals in CALIFA and SDSS spectra.

4.1. SSP spectral models

Base GM is a combination of publicly available SSP spectra from Vazdekis et al. (2010)⁷, which start at an age of 63 Myr, with the González Delgado et al. (2005)⁸ models for younger ages. The former are based on stars from the MILES library (Sánchez-Blázquez et al. 2006), while the latter relies on the synthetic stellar spectra from the granada library (Martins et al. 2005). Only minor adjustments were needed to match these two sets of models. First, we smoothed the granada models to the spectral resolution of MILES (2.5 Å FWHM, cf. Beifiori et al. 2011 and Falcón-Barroso et al. 2011), and then multiplied them by a factor of 1.05, estimated from a direct comparison of the predicted spectra in the 63 Myr to 1 Gyr range, where these models agree very well in both continuum shape and absorption features. These independently derived SSP spectra complement each other, providing a base suitable for spectral fitting analysis of galaxies of all types.

The evolutionary tracks in base GM are those of Girardi et al. (2000), except for the youngest ages (1 and 3 Myr), which are based on Geneva tracks (Schaller et al. 1992; Schaerer et al. 1993a,b; Charbonnel et al. 1993). The specific subset of SSPs used in the spectral fits comprises 39 ages between 1 Myr and 14 Gyr. Four metallicities are included: Z = 0.2, 0.4, 1, and 1.5 solar. The Salpeter IMF is adopted. Like the other bases discussed next, GM does not adequately cover the blue horizontal branch, both because it does not reach very low Z and because, as thoroughly discussed by Ocvirk (2010), these models perform poorly in this phase anyway.

Base CB is an updated version of the BC03 models (Charlot & Bruzual 2007, priv. comm.), replacing STELIB by spectra from the MILES and granada libraries. The IMF is that of Chabrier (2003), and the evolutionary tracks are those by Alongi et al. (1993), Bressan et al. (1993), Fagotto et al. (1994a,b), and Girardi et al. (1996), collectively labeled “Padova 1994” by BC03. N_⋆ = 160 SSPs comprising 40 ages in the 1 Myr to 14 Gyr range and four metallicities (0.2, 0.4, 1, and 2.5 Z_⊙) were chosen for the spectral fits.

Finally, base BC is built from the “standard” BC03 models, whose ingredients are identical to those of base CB except for the stellar library, which is STELIB. Base BC allows us to check whether differences between results obtained with bases GM and CB are associated to isochrones and/or IMF instead of stellar libraries.

These three bases cover the same age range, and all have the same number of metallicities. All bases start at Z = 0.004 (0.2 Z_⊙)⁹, but GM stops at 1.5 Z_⊙, while both CB and BC extend to 2.5 Z_⊙. Because they are both based on stars not far from the Sun, neither MILES nor STELIB actually contain many stars as metal-rich as 2.5 Z_⊙, hence the latter models should be interpreted with care. On the other hand, the Z ≤ 1.5 Z_⊙ limit of base GM may be too low to accomodate the metal-rich inner regions of galaxies, possibly leading to saturation effects. Base GM also differs from the others in IMF and isochrones. The differences in opacities and equation of state between Padova 2000 (GM) and 1994 (CB and BC) tracks produce somewhat warmer (hence bluer) stars in the red giant branch in the former than in the latter. In principle, one may thus expect to obtain older ages with GM than with CB and BC, but the diffences in color are very small in the spectral range of our data, only becoming noticeable at longer wavelengths (see Fig. 2 of BC03).

4.2. Data: 98291 spectra from 107 CALIFA datacubes

The experiments described below were based on data with the same characteristics and pre-processed in the exact same way as the data for CALIFA 277. The difference with respect to the simulations presented in the first part of this paper is that we now use 107 galaxies, out of which a total of 98291 spectra corresponding to individual spaxels (or Voronoi zones) were extracted following the methodology outlined in Paper I.

starlight fits of these spectra were carried out using the three bases described above. Our interest here is to identify and quantify differences in the results obtained with different model SSPs, so we overlook the spatial information and treat each spectrum as an individual. We now use this massive dataset to compare the results obtained with bases GM, CB, and BC.

4.3. Comparison of global properties

Figure 9 compares results obtained with bases CB (horizontal axis) and GM (vertical). The top panels compare values of ⟨log t⟩_L, log ⟨Z⟩_M, A_V and log M′ for all 98921 individual zones, grayscale-coding by the density of points. Note the wide ranges of values in all panels, illustrating the diversity of properties spanned by the data. Each panel quotes the mean Δ and its standard deviation, where Δ = property(GM ) − property(CB).

GM-based initial stellar masses are higher than CB ones by 0.14 dex on-average, reflecting their different IMFs¹⁰. Apart from this offset, the two masses agree to within 0.16 dex. Mean ages and extinctions also agree well, with dispersions of 0.14 dex in ⟨log t⟩_L and 0.06 mag in A_V, and insignificant offsets. Metallicities are obviously related, but with a substantial scatter due to a combination of the difference in evolutionary tracks and IMF, coarse Z grids, intrinsic difficulties in estimating metallicities, and the aforementioned difference in maximum Z.

To estimate the effect of this in-built difference we carried out experiments replacing the 2.5 Z_⊙ SSPs in CB by models interpolated to 1.5 Z_⊙, thus matching the Zs in GM. This modified CB base produces log ⟨Z⟩_M values that differ by σ_Δ ~ 0.15 dex from the GM-based ones, slightly better than the 0.18 dex obtained with the original CB. Leaving aside caveats associated to the interpolation of SSP spectra, this result indicates that the difference in maximum Z between GM and CB has a relatively modest impact on the difference between log ⟨Z⟩_M values obtained with these two bases, so that differences in evolutionary tracks and IMF are also relevant.

The bottom panels repeat the comparison of the same four properties, but now for galaxy-wide values, so there are only 107 points¹¹. These cleaner plots confirm the results from the zone-by-zone comparison. Naturally, galaxy-wide averages are much more robust than values for individual zones, so the σ_Δ dispersions are significantly lower than those quoted in the top panels. ⟨Z⟩_M estimates distribute around the identity line up to Z_⊙, where the relation bends down due to the differences in the maximum Z between the bases. In these cases all one can be sure of is that metalicities are over-solar.

Figure 10 compares results for bases BC (vertical axis) with CB (horizontal). Masses, ages, and extinctions agree to within 0.15 dex, 0.16 dex, and 0.09 mag, respectively, with a tendency of A_V to be 0.04 mag smaller in BC fits. This relatively good level of agreement does not apply to mass weighted metallicities, which differ substantially. Given that these two sets of models coincide in isochrones, IMF and evolutionary synthesis code, this disagrement must be related to the stellar libraries. Because of the several improvements in MILES over STELIB (Sánchez-Blázquez et al. 2006), metallicities from the CB and GM bases should be considered more reliable than those obtained with BC (see also Koleva et al. 2008). This is in line with previous empirical tests by González Delgado & Cid Fernandes (2010), who showed that MILES based SSP spectra perform better than STELIB-based ones in matching the metallicities of star clusters derived from spectroscopy of individual stars or CMD analysis of the same systems.

Figure 11 shows histograms of the pair-wise Δs: GM − CB (top), BC − CB (middle), GM − BC (bottom). The width of these histograms (the σ_Δ values listed in each panel) gives a measure of the uncertainties related to the choice of evolutionary synthesis models for the spectral fits. As usual, mass-weighted properties present larger uncertainties than luminosity-weighted ones. Metallicities, for instance, have dispersions in Δlog  ⟨Z⟩_L of 0.14 and 0.19 dex (GM versus CB, and BC versus CB, respectively), lower than the 0.18 and 0.27 dex found for Δlog  ⟨Z⟩_M.

The overal conclusion of these comparisons is that from a statistical perspective, physical properties obtained from a spectral synthesis analysis do not depend strongly on the choice of SSP models. Mean ages and extinctions retrieved with the three models investigated here all agree to within relatively small margins, and the same is true for stellar masses once in-built differences in IMFs are accounted for. Metallicities are generally harder to estimate, and our experiments confirm this. Unlike other properties, mean Zs are sensitive to the choice of models. Base BC gives the most discrepant results, highlighting the influence of stellar libraries.

Rounding up numbers, we find that uncertainties related to the SSP models adopted in the analysis are about 0.15 dex in ⟨log t⟩_L, 0.23 dex in ⟨log t⟩_M, 0.15 dex in log M′, and 0.1 mag in A_V. Excluding BC models, log ⟨Z⟩_L and log ⟨Z⟩_M have uncertainties of 0.14 and 0.18 dex, respectively.

4.3.1. Models versus data- and method- related uncertainties

How do these uncertainties compare with those related to the data and synthesis method derived in the first part of this article?

To answer this, we focus on the OR1 and OC002 simulations described in Sect. 2.1, because they represent random noise and shape-calibration errors in CALIFA. Comparing the σ_Δ values in Fig. 11 with those for the OR1 simulations in Fig. 1, we see that the former are roughly twice as high as the latter. In fact, the σ_Δ values in the top panels of Fig. 11 are closer to those obtained in the simulations where the noise was artificially amplified by a factor or 2 (OR2 runs in Table 2). Shape-related uncertainties also lead to uncertainties somewhat smaller than those found in this section, with the exception of A_V, whose uncertainty is 0.16 mag in the OC002 runs but 0.1 mag or less in ΔA_V histograms in Fig. 11. A visual comparison of Figs. 1 and 11 corroborates the conclusion that uncertainties related to differences among SSP models are larger than those related to data and method.

These comparisons, however, are not entirely fair, since model-related uncertainties were evaluated from ~10⁵ zones of 107 galaxies, while data- and method- related uncertainty estimates were based on 1638 zones of a single galaxy, perturbed ten times in the simulations. Repeating Fig. 11 for CALIFA 277 alone, we still find σ_Δ values higher than those found in Sect. 3, but by a smaller margin. The important message, therefore, is not so much that variance among current SSP spectral models leads to larger uncertainties than those related to data and method, but that both types of uncertainties are of the same order.

4.4. Quality of spectral fits

In principle, one could use the quality of the spectral fits as a criterion to favor one base over another. In practice, however, all bases studied here provide ~ equally good fits to the data. This is illustrated in Fig. 12, which shows the distributions of two figures of merit. Green, red, and blue correspond to bases GM, CB, and BC, respectively. The histograms are nearly identical. Considering all 98291 fits, the mean (median) values of the relative spectral deviation (Eq. (6) of Paper I) are (3.92), 4.34 (3.93) and 4.37 (3.94) % for bases GM, CB, and BC respectively, while the corresponding mean (median) χ² per fitted flux are 1.14 (0.96), 1.15 (0.97), and 1.18 (0.99). Formally, this puts base GM in first place, with CB in second and BC in third, but these statistics reveal that it is clearly not possible to favour one or another model in terms of fit quality. Because of their more complete spectral libraries, we favor bases GM and CB over BC.

Figure 13 provides a more direct visualization of the quality of the spectral fits. The top panel shows (in black) the average of all zone spectra (for examples of fits to individual zones see González Delgado et al., in prep.). Each spectrum is first divided by the median flux in the 5590–5680 Å normalization window, such that all zones weight equally in the average. The average synthetic spectra are shown in green, red, and blue for fits with bases GM, CB, and BC, respectively, vertically offset for clarity. The corresponding residual spectra are shown in the bottom of the same panel, but can hardly be distinguished on this scale. This similarity persists even after zooming-in by a full order of magnitude in the flux scale, as shown in Fig. 13b. Gaps in this zoomed plot correspond to our generic emission line mask, i.e., λ-windows that were ignored in the starlight fits. Note that, because of bad-pixel flags, augmented by windows around telluric features (Paper I), the number of specra involved in these averaged residuals varies for each λ, as shown in the bottom panel.

4.5. Spectral residuals and comparison with SDSS

Clearly, the amplitude of the average spectral residuals is very small. Whereas the typical relative deviation for fits of individual zones is % (Fig. 12), for the average spectra in Fig. 13 it is a tiny %. This four-fold reduction indicates that random noise is responsible for a large share of the residuals in spectral fits of individual zones. Yet, if the residuals were truly random, this reduction should be by a factor of the order of , not just 4. Residuals may be small, but they are definitely real.

The systematics of the spectral residuals are visible in Fig. 13b. These deviations must be related to problems in the method, data, and/or models.

The narrow positive peaks in the average residual spectrum are due to our fitting method. They are in fact emission lines not included in our generic emission line mask. [NeIII]λ3869, H8λ3889, [NI]λ5199, and [OI]λ6364 are the clearest ones in Fig. 13b, but weaker lines such as H9λ3835 and HeIλ6678 are identified upon closer inspection. The conservative clipping parameters adopted in the fits (see Paper I) prevent starlight from authomatically clipping these lines completely. Other starlight-based studies circumvented this problem by building emission line masks taylored for each spectrum individually (Mateus et al. 2006; Asari et al. 2007), a refinement that is yet to be implemented in our CALIFA pipelines.

The most noticeable features in Fig. 13b, however, are not the emission lines, but the broad O_λ < M_λ trough around Hβ and the O_λ > M_λ bumps in the 5800–6000 Å region. To understand whether these features are particular to CALIFA data, we compiled SDSS spectra for 50 galaxies in common with the current sample and analyzed them using the same masks, wavelength range, sampling, and starlight configuration as used for CALIFA. These spectra were then averaged as in Fig. 13, and compared with the spectral statistics obtained from the nuclear CALIFA spectra of the same 50 galaxies. Figure 14 shows the results of this comparison. The figure follows the same layout as Fig. 13, with spectral statistics for the 50 CALIFA nuclei and the corresponding SDSS spectra plotted in the left and right panels, respectively. The spikier appearance of the SDSS residuals is due to its higher spectral resolution (3 versus 6 Å FWHM), but this is irrelevant for the broad residual features discussed above.

The humps around the ~5800 Å region of CALIFA residuals are clearly absent in the SDSS. These features are seen as shallow slanted bands in Fig. 16 of Husseman et al. (2013), where residuals for 100 galaxies were sorted by redshift. We have tracked the origin of this feature to an imperfect removal of telluric NaI bands around 5900 Å in the observed frame. In analogy with what happened with the SDSS data releases, as CALIFA data are used imperfections like this will be identified and, when possible, corrected for in future updates of the data base.

The trough around Hβ, however, is present in both SDSS and CALIFA spectra. This problem has been noticed before (Cid Fernandes 2006; Panter et al. 2007; Walcher et al. 2009), but not satisfactorily explained to date. The fact that this feature appears in both data sets suggests that the models overpredict the continuum in a ~200 Å wide region centered around Hβ. On the other hand, the fact that the problem is present for both MILES- and STELIB-based models suggests that its origin is more likely related to data-calibration issues. This latter view is supported by studies such as the one by Ferré-Mateu et al. (2012), where starlight fits of William Herschel Telescope long-slit spectra of massive compact galaxies with the Vazdekis et al. (2010) models (included in our GM base) did not reveal the Hβ trough seen in Fig. 13b. While a full investigation of this matter is beyond the scope of this paper, we note that the trough tends to be more pronounced in younger systems, as previously reported by Cid Fernandes (2006, see his Fig. 5).

The negative residual seen in the region of Mg lines near 5175 Å is almost as strong as the Hβ trough. This deficit is due to differences in the abundances of α elements between galaxies (particular the central regions of massive ones) and stars in the stellar libraries. That this residual is stronger in Fig. 14 (nuclear spectra) than in Fig. 13 (all zones) is most likely due to α/Fe gradients. Walcher et al. (2009) reported steps toward incorporating α/Fe in spectral synthesis analysis.

4.6. Public distribution of base GM spectra

Given the good performance of the GM base in the tests discussed in this section, we make these spectra available to the community¹². We stress that these models were produced by González Delgado et al. (2005) and Vazdekis et al. (2010), to whom all credit should be given. Our distribution just puts these two sets of models together.

5. Summary

We have performed an extensive study of the uncertainties in the stellar population properties derived from starlight fits of optical integral field spectra from the CALIFA survey. In the first part of this investigation, Monte Carlo simulations were used to explore the effects of random noise, continuum-shape calibration uncertainties and degeneracies in the spectral synthesis. The simulations are based on the data and starlight fits of 1638 zone spectra of the Sb galaxy CALIFA 277 (NGC 2916). Uncertainties in mean age, metallicity, extincion, stellar mass, and SFHs were evaluated both at the level individual zones (treating the datacube as a collection of unrelated spectra) and exploring the statistical power of IFS by means of radial profiles and other spatial averages. Second, we explored the variations in physical properties caused by using different sets of models for SSP spectra. We compared physical properties and spectral residuals obtained from fits with three different SSP bases: (a) a combination of the Gonzalez Delgado et al. (2005) and Vazdekis et al. (2010) models (GM); (b) preliminary Charlot & Bruzual models using the granada and MILES libraries (CB); (c) the standard Bruzual & Charlot (2003) models with STELIB (BC). This comparison employed ~10⁵ spectra from 107 CALIFA galaxies.

Our main results can be summarized as follows:

• 1.

  The level of noise in CALIFA spectra leads to one-σ uncertainties in mean ages and metallicities of ~0.08 dex when weighting in luminosity and ~0.15 dex weighting by mass. Stellar masses are uncertain at the ~0.08 dex level, while σ(A_V) is about 0.06 mag.

• 2.

  Shape-related calibration uncertainties of 0.05 mag in g − r produce uncertainties on the same order as those induced by random noise, except for A_V, which becomes uncertain at the level of σ ~ 0.16 mag and biased by + 0.05 mag.

• 3.

  Time-resolved SFHs are higher-order products of the spectral synthesis, inevitably subjected to larger uncertainties than global properties. For instance, we found that, for individual zones, stellar masses formed in broad age ranges have uncertainties about three times larger than those of the mass formed over all times. While SFHs of individual zones do carry useful information, a statistical approach is more advisable.

• 4.

  Spectral synthesis results are best used in comparative studies involving large samples, and IFS data fit this requirement perfectly. Spatial averaging reduces uncertainties while preserving enough information on the history and structure of stellar populations. Radial profiles of mean ages, metallicities, mass density, and extincion are found to be very robust due to the large number statistics of datacubes. For instance, uncertainties in the data and limitations of the synthesis have a negligible impact upon the negative age and metallicity gradients detected in CALIFA 277. Similarly, SFHs averaged over different spatial regions are much more stable than those obtained for individual zones.

• 5.

  Mean ages, extinction, and stellar masses obtained with GM, CB, and BC models are reasonably consistent with one another, with one-σ differences in ⟨log t⟩_L, ⟨log t⟩_M, log M′, and A_V of ~ 0.15, 0.24, 0.15 dex, and 0.08 mag, respectively. GM and CB produce luminosity- (mass-) weighted metallicities consistent to within 0.14 (0.18) dex rms, while BC-based metallicities hardly correlate with the others. The variations induced by the use of different SSP models are larger than uncertainties due to data and method by roughly a factor of 2, highlighting the important role of this ingredient in the whole analysis.

• 6.

  The quality of the spectral fits is very nearly the same whichever set of SSP models is used. Spectral residuals are on the order of 4% for individual zones and 1% when averaged over all zones and galaxies. The strongest feature in the residual spectra is a broad but shallow (~3%) trough around Hβ, also present in SDSS spectra, and its origin remains unidentified. Other systematic features in the residual spectra were tracked to incomplete masking of weak emission lines, imperfect removal of telluric features, and the lack of α enhanced models in the base.

These results, along with the tables and figures in this paper, fulfill our goal of assessing the uncertainties on the several products of starlight applied to CALIFA datacubes presented in Paper I. The combination of these spatial and temporal diagnostics of stellar populations with the diversity of galaxy types in the CALIFA survey offers new ways of studying galaxy evolution, that we aim to explore in other publications by our collaboration, and this paper provides quantitative and qualitative guidelines on how to interpret them appropriately.

---

Acknowledgments

CALIFA is the first legacy survey being 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 infrastructures. We also thank the CAHA staff for the dedication to this project. R.C.F. thanks the hospitality of the IAA and the support of CAPES and CNPq. A.L.A. acknowledges support from INCT-A, Brazil. BH gratefully acknowledges the support by the DFG via grant Wi 1369/29-1. Support from the Spanish Ministerio de Economia y Competitividad, through projects AYA2010-15081 (PI RGD), AYA2010-22111-C03-03 and AYA2010-10904E (SFS), AYA2010-21322-C03-02 (PSB) and the Ramón y Cajal Program (S.F.S., P.S.B. and J.F.B.), is warmly acknowledged. We also thank the Viabilidad, Diseño, Acceso y Mejora funding program, ICTS-2009-10, for funding the data acquisition of this project. Funding for the SDSS and SDSS-II has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, the U.S. Department of Energy, the National Aeronautics and Space Administration, the Japanese Monbukagakusho, the Max Planck Society, and the Higher Education Funding Council for England. The SDSS Web Site is http://www.sdss.org/. The SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions. The Participating Institutions are the American Museum of Natural History, Astrophysical Institute Potsdam, University of Basel, University of Cambridge, Case Western Reserve University, University of Chicago, Drexel University, Fermilab, the Institute for Advanced Study, the Japan Participation Group, Johns Hopkins University, the Joint Institute for Nuclear Astrophysics, the Kavli Institute for Particle Astrophysics and Cosmology, the Korean Scientist Group, the Chinese Academy of Sciences (LAMOST), Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, Ohio State University, University of Pittsburgh, University of Portsmouth, Princeton University, the United States Naval Observatory, and the University of Washington.

References

All Tables

All Figures