Resolving galaxies in time and space: I:

Resolving galaxies in time and space: I:

Applying STARLIGHT to CALIFA data cubes
Key Words.:
galaxies: evolution – galaxies:stellar content – galaxies: general– techniques: imaging spectroscopy



Aims:Fossil record methods based on spectral synthesis techniques have matured during the past decade, and their application to integrated galaxy spectra fostered substantial advances on the understanding of galaxies and their evolution. Yet, because of the lack of spatial resolution, these studies are limited to a global view, providing no information about the internal physics of galaxies.

Methods:Motivated by the CALIFA survey, which is gathering Integral Field Spectroscopy (IFS) over the full optical extent of 600 galaxies, we have developed an end-to-end pipeline which: (i) partitions the observed data cube into Voronoi zones in order to, when necessary and taking due account of correlated errors, increase the signal-to-noise ratio, (ii) extracts rest-framed spectra, including propagated errors and bad-pixel flags, (iii) feeds the spectra into the starlight spectral synthesis code, (iv) packs the results for all galaxy zones into a single FITS or HDF5 file, (v) performs a series of post-processing operations, including zone-to-pixel image reconstruction and unpacking the spectral and stellar population properties derived by starlight into multi-dimensional time, metallicity, and spatial coordinates. This paper provides an illustrated description of this whole pipeline and its many products. Using data for the nearby spiral NGC 2916 as a show case, we go through each of the steps involved, and present a series of ways of visualizing and analyzing this manifold. These include 2D maps of properties such as the velocity field, stellar extinction, mean ages and metallicities, mass surface densities, star formation rates on different time scales and normalized in different ways, plus 1D averages in the temporal and spatial dimensions, leading to evolutionary curves and radial profiles of physical properties. Projections of the stellar light and mass growth onto radius-age diagrams are introduced as a means of visualizing galaxy evolution in time and space simultaneously, something which can also be achieved in 3D with snapshot cuts through the cubes.

Results:The results provide a vivid illustration of the richness of the combination of IFS data with spectral synthesis, of the insights on galaxy physics provided by the variety of diagnostics and semi-empirical constraints obtained, as well as a glimpse of what is to come from CALIFA and future IFS surveys.


1 Introduction

Last decade technology has lead to massive surveys either in imaging or spectroscopy modes that have been very successful in providing information of the spectral energy distribution or central/integrated spectra for a large number of galaxies (e.g. 2dF, Folkes et al. 1999; SDSS, York et al. 2000; COMBO-17, Bell et al. 2004; ALHAMBRA, Moles et al. 2008; COSMOS, Ilbert, et al 2009). These data have allowed significant insight on the integrated properties of galaxies, as well as on the evolution of global quantities such as the mass assembly (Pérez-Gonzalez et al. 2008) and the star formation history (SFH) of the universe (Lilly et al. 1996; Madau et al. 1996). However, these surveys are limited in either spectral or spatial resolution. Imaging provides spatial information, but at the expense of coarse spectral coverage (typically broad band filters), limiting the amount of information on the stellar populations, and providing little or no information on the gas (emission lines) nor the kinematics of the galaxy. Surveys based on integrated spectra, offer a richer list of diagnostics of the gaseous and stellar components, but lack spatial resolution and suffer from aperture effects (one spectrum per object, centered at the nucleus and not covering the full galaxy). Moreover, integrated galaxy spectra do not allow to, for example, isolate morphological components (bulge, disc, bars), map the effects of mergers, trace secular processes such as stellar migration, and other features of galaxy formation and evolution which can only be observationally tackled with a combination of imaging and spectroscopic capabilities.

The near future will see a proliferation of Integral Field Spectroscopy (IFS) surveys, which will allow a detailed look at physics within galaxies, as opposed to the global view offered by integrated light surveys. CALIFA, the Calar Alto Legacy Integral Field Area survey is a pioneer in this area (Sánchez et al. 2012), and others will follow soon, like SAMI (Croom et al 2012), VENGA (Blanc et al. 2009), and MaNGA (Bundy et al., in prep). CALIFA is obtaining IFS for 600 nearby galaxies () of all types, mapping stellar populations, ionized gas and their kinematics across the full optical extent of the sources. The number of spectra to be generated is of the same order of the whole SDSS ().

Clearly, however, these will not be “just another million galaxy spectra”. Processed through the machinery of fossil record methods, which recover the history of a stellar system out of the information encoded in its spectrum (Walcher et al. 2011 and references therein), CALIFA will provide valuable and hitherto unavailable information on the spatially resolved SFH of galaxies. An example of what can be done was presented by Pérez et al. (2012), where we have used data on the first 105 galaxies observed by CALIFA in conjunction with the spectral synthesis code starlight (Cid Fernandes et al. 2005) to study the spatially resolved mass assembly history of galaxies. We find that, for massive galaxies, the inner parts grow faster than the outer ones, a clear signature of inside-out growth, and that the relative growth rate depends on the stellar mass of the galaxy, with a maximum relative growth efficiency for intermediate masses ().

The application of spectral synthesis methods to IFS datacubes opens new ways of studying galaxies and their evolution. In practice, however, this apparently simple extension of work done for spatially integrated spectra involves a series technical details and methodological issues, as well as challenges in the visualization and representation of the results.

The purpose of this paper is to present a detailed account of all steps in the processing of CALIFA datacubes to the multidimensional products derived from the starlight analysis. We thus focus basically on methodology, leaving the exploration of results for subsequent studies. Paper II (Cid Fernandes et al. 2013) presents a thorough study of the uncertainties involved in these products, including the effects of random noise, spectral shape errors, and evolutionary synthesis models used in the analysis. Despite the emphasis on CALIFA and starlight, the issues discussed here and in Paper II are of direct interest to spectral synthesis analysis of IFS data in general.

The paper is organized as follows. Section 2 presents the data. Section 3 provides a detailed description of the pre-processing steps which prepare the reduced data for a full spectral fitting analysis. Section 4 reviews the starlight code and the ingredients used in our analysis. It also introduces the pycasso pipeline, used to pack and analyze the results of the synthesis. Our main results are presented in Section 5, which uses the nearby spiral NGC 2916 (galaxy 277 in the CALIFA mother sample) as a show case. This long section starts with considerations on how to evaluate the quality of the spectral fits, and then proceeds to the presentation of the synthesis products in different ways, like 2D maps of the stellar light and mass, extinction, mean ages and metallicities, and star formation rates maps on different time-scales, 1D maps on either the spatial or temporal dimensions, and attempts to visualize the evolution of galaxies in time and space simultaneously. The emphasis is always on methodological issues, but the examples naturally illustrate the richness of the multidimensional products of this analysis. Finally, our results are summarized in Section 6.

2 Data

The goals, observational strategy and overall properties of the survey were described in Sánchez et al. (2012). In summary, CALIFA’s mother sample comprises 939 galaxies in the 0.005–0.03 redshift range, extracted from the SDSS imaging survey to span the full color-magnitude diagram down to M mag, and with diameters selected to match the field of the IFU instrument ( in diameter). A thorough characterization of the sample will be presented in Walcher et al. (in prep). CALIFA will observe galaxies from this mother sample, applying purely visibility-based criteria at each observing night. The final sample will be representative of galaxies in the local universe, statistically significant and well selected.

The data analyzed in this paper were reduced using the CALIFA Pipeline version 1.3c, described in the Data Release 1 article (Husemann et al. 2013). Each galaxy is observed with PMAS (Roth et al. 2005, 2010) in the PPAK mode (Verheijen et al. 2004; Kelz et al. 2006) at the 3.5m telescope of Calar Alto Observatory and with two grating setups: the V500, covering from to 7500 Å at FWHM Å resolution, and the V1200 (–4600 Å, FWHM Å). The wide wavelength coverage of the V500 is ideal for our purposes, but the blue end of its spectra, which contains valuable stellar population tracers such as the 4000 Å break and high order Balmer lines, is affected by vignetting over a non-negligible part of the field of view. To circumvent this problem we worked with a combination of the two gratings, whereby the data for Å comes from the V1200 (which does not suffer from this problem at this wavelength range) and the rest from the V500. The spectra were homogenized to the resolution of the V500.

For each galaxy the pipeline provides a cube with a final sampling of . The cube comprises absolute calibrated flux densities at each spaxel (), corrected by Galactic extincion using the Schlegel, Finkbeiner & Davies (1998) maps and the extinction law of Cardelli, Clayton & Mathis (1989) law with . Besides fluxes, the reduced cubes contain carefully derived errors (), and a bad-pixel flag () which signals unreliable flux entries, due to, say, bad-columns, cosmic rays and other artefacts (see Husemann et al. 2013 for details).

3 Pre-processing steps: From datacubes to STARLIGHT input

A series of pre-processing steps are applied to the reduced datacubes, with the general goal of extracting good quality spectra for a stellar population analysis with starlight or any other similar code. These are:

  1. Definition of a spatial mask

  2. Outer mask

  3. Refinements of the (flag) and (error) spectra

  4. Rest-framing

  5. Spatial binning

  6. Resampling in

This section describes each of these operations.

Figure 1: (a) SDSS stamp for NGC 2916 (CALIFA 277). As in other images throughout this paper, the circles mark and 2 Half Light Radius (HLR). (b) CALIFA image in the rest-frame Å interval, after applying the spatial mask. (c) Map of the at 5635 Å. (d) map after the Voronoi binning. (e) Number of spaxels in each Voronoi zone. Note that no spatial binning () was necessary throughout the inner HLR. (f) Percentage of bad pixels in the 3800–6850 Å range. North is up and east is towards the left in all images.

3.1 Spatial masks

The first step is to mask out foreground stars, artefacts and low regions. As a starting point, candidate spurious sources are identified applying SExtractor (Bertin & Arnouts 1996) to the SDSS -band image. The masked image is then matched to the CALIFA field of view using the WCS information available in the headers to reproject the SDSS masks to the scale, orientation and pixel size of the CALIFA cubes. To improve upon this first approximation we perform a visual inspection of each galaxy, correcting its mask interactively whenever necessary. For instance, we manually fixed a few unmasked pixels which remained in the borders of each masked object (probably due to a combination of the slightly better PSF of the SDSS images and the accuracy of the CALIFA WCS information used in the matching). Also, in a number of cases we found that the initial mask included bright Hii regions, while in others very dim point-like sources can be distinguished in the SDSS image, but have no clear counterpart in the CALIFA image nor spectra. In all cases, the masks were corrected by hand.

The resulting masks are free from spurious spaxels. Unlike the other operations described below, this step could not be fully automated. For the size of CALIFA, this is still a feasible, albeit laborius, task. A last step in defining the spatial mask is to impose a spectral quality threshold based on the of the data. This step has been fully automated in the pipeline described next.

3.2 Flags, errors and segmentation: The qbick pipeline

The next steps (items ii to vi) were all packed in a fully automated package, named qbick, built specifically for CALIFA but sufficiently general to handle any FITS format to pre-process datacubes for a 3D spectral fitting analysis.

The flags produced by the reduction pipeline ( for good pixels and for flagged ones) do a good job in identifying problematic pixels. For uniformity, we enforce around the strongest sky lines seen in the CALIFA data: HgI 4358 Å, HgI 5461 Å, OI 5577 Å, NaI D (around 5890 Å), OI-OH (6300 Å), OI 6364 Å, and the B-band atmospheric absorption (Sánchez et al. 2007). We have further set when the value of the error is very large (specifically, when ), but this is merely a cosmetic choice, as these pixels would not influence the spectral fits anyway. Anomalously small values, on the other hand, can have a disproportionately large weight on the analysis. In the very few cases where this happens, the pixels were flagged. These refinements are useful to circumvent potential problems in the spectral fitting analysis, specially when the whole data set is processed in “pipeline mode”.

After these fixes all spectra are rest-framed using the redshift derived by the reduction pipeline for the central spectrum (Sánchez et al. 2012), included in the original headers. This step includes a correction for the energy, time and bandwidth changes with redshift.

We used the (rest wavelength) window between 5590 and 5680 Å to estimate a ratio, defining the signal as the mean flux and the noise as the detrended standard deviation1 of , both excluding flagged pixels.2 Regions outside the outer iso-contour were added to the spatial mask, and discarded from the analysis. In a few cases this threshold was lowered due to low surface brightness.

Fig. 1b shows the Å image of CALIFA 277, after the final spatial mask. The hole above the nucleus corresponds to the foreground source seen in the SDSS image (panel a). The hexagonal pattern of the fiber bundle is clearly visible, since in this example over most of the field of view (panel c).

The 5635 Å image is also used to compute the galaxy’s Half Light Radius (HLR), adding up fluxes in spaxels sorted by distance to the nucleus and picking the radius at which half of the total flux is reached. For CALIFA 277, HLR . This datacube-based HLR may not be the ideal measure of the size of a galaxy, but it is the most convenient metric for the present analysis for the simple reason that it is based on the exact same data. In Fig. 1 and all images throughout this paper circles are drawn at and 2 HLR to guide the eye.

The following step (item v in our list) consists of applying a segmentation structure to the data, grouping spaxels into spatial zones. All spaxels belonging to the same zone are assigned a common integer ID label, unique for each zone. qbick stores the zone-to- tensor needed to map zones to pixels and vice-versa. In this paper we use spatial binning to increase the for the spectral analysis. Other applications may prefer to group spaxels according to, say, morphology or emission line oriented criteria. qbick admits an externally provided segmentation map to allow for such choices.

The zoning scheme adopted here is based on the Voronoi tesselation technique, as implemented by Cappellari & Copin (2003), but with modifications to account for correlated errors (Section 3.3). The code is feed with the signal and noise images in the Å range discussed above. The target is set to 20. In practice, most individual spaxels inside 1 HLR satisfy this condition. Fig. 1d shows the map of CALIFA 277 after the spatial binning. For this galaxy, the 3649 non-masked spaxels were grouped into 1638 zones, the overwhelming majority (1527) of which are not binned at all. As seen in the map of the number of spaxels per zone (, Fig. 1e), spatial binning effectively starts at HLR. Discounting the single spaxel zones, the median is 10, but values up to 92 are reached in the outskirts.

Once the spatial binning map is defined, the flux, error and flag spectra are computed for each zone. This operation must account for the flagged pixels. The total flux in zone , containing spaxels whose fluxes are , is given by


where for flagged pixels () or 1 otherwise (), and is the number of useful spaxels at wavelength in the zone. The zone flux is thus times the average flux over all non-flagged pixels. When none of the fluxes are flagged, this equation simply sums all the fluxes in the zone, whereas when this is not the case, flagged entries are replaced by the average of non-flagged ones.

This scaling recipe is reasonable as long as the number of flagged pixels is not comparable to . For example, in a zone of spaxels, one would not like to take the resulting seriously if 8 of them are flagged. To deal with this issue we flag away entries whenever less than half of the spaxels in a zone have credible data at the given :


For our example galaxy, the fraction of flagged pixels in the 3800–6850 Å interval used in the spectral synthesis analysis ranges from 5 to 8%, with a median of 6% (see Fig. 1f).

The computation of the error in follows a recipe similar to eq. 1, but adding in quadrature


where accounts for correlated errors (see below).

qbick also produces the total spatially integrated spectrum, following these same prescriptions. This is useful to compare the analysis of the whole to the sum of the analysis of the parts (which we do in Section 5.9).

Finally, all the spectra are resampled to 2 Å. The resulting zone spectra are stored in ASCII and/or FITS files, with flux , error and flag columns necessary for a careful spectral analysis. All the quality control images produced in this process are packed in a single FITS file along with the technical parameters adopted in the different steps ( clipping, segmentation map, etc.)

3.3 Spatial binning and correlated errors

Figure 2: Ratio of the measured over the idealized noise for Voronoi zones comprising spaxels. The idealized noise is that derived under the assumption of uncorrelated errors, while the measured one is obtained from the detrended standard deviation of fluxes in a Å range around 5635 Å. Gray dots represent values for 99267 spaxels from 109 galaxies. Circles show the median for each , and the solid line shows equation 4. The top panel show the histogram of -values (notice the logarithmic scale).

As already mentioned, the data cubes have a spatial sampling of . Because of the three-fold dithering pattern used in the observations, each spaxel has a varying contribution from a number of fibers. Since one fiber can contribute to more than one spaxel, the noise in a given spaxel is partially correlated with that of its neighbors. This implies that when spectra from adjacent spaxels are added together the noise in the resulting spectrum has to be calculated taking spatial covariances into account. As an illustrative extreme example, consider completely correlated (ie., identical) spaxels in a zone. Clearly, co-adding their spectra would not result in the desired increase in the , as both signal and noise scale linearly with . Unless told otherwise, however, the zoning algorithm assumes that the spaxels are independent, and hence that the noise in the zone is given by , resulting in an improvement of the of the co-added spectrum by , whereas in fact it has not changed at all.3

An empirical correction was devised to account for the global behaviour of the spatial correlation of the signal in CALIFA data. The way we proceed is as follows: First we define initial Voronoi zones with the Cappellari & Copin (2003) algorithm. Because it assumes uncorrelated input, the code presumes the zone error is . We then measure the “real” noise () in the resulting zone summed spectra from the (detrended) rms in the Å range. The ratio of these two numbers, which we call , is represented versus the number of spaxels in each zone in Fig. 2. Small grey dots are values for the individual zones, circles represent the median value computed at each , while the shaded region represents the 16 to 84 percentiles. There are nearly points in this plot, obtained from 244196 spaxels in 109 galaxies (all those observed with the V500 and V1200 setups up to May 2012). The top panel shows the histogram of .

The rise and flattening of the - relation reflects the expected behavior. Zones with small group adjacent and highly correlated spaxels, producing the initial rise of . As farther apart spaxels are aggregated, the errors become less correlated, leading to the flattening observed. The scatter in the relation is likely related to the way the Voronoi algorithm assemble spaxels. Also, the three-fold dithering pattern is not completely uniform, so the amount of correlation is not exactly the same across the face of the fiber bundle. Still, the relation portrayed in Fig. 2 is tight enough to derive an excellent statistical correction.

The solid line represents the fit of the following function:


where the fit was carried out considering only the filled circles. We can now apply this correction directly in the zoning code, replacing its idealized computation of the noise by


which accounts for the effects of spatially correlated errors. This correction was applied by (cf. equation 3). An independently derived correction very similar to equation 4 is presented in Husemann et al. (2013).

This scheme naturally leads to larger zones than those obtained under the (erroneous) hypothesis of independent spaxels, but ensures that the target of 20 is actually achieved in all but the outermost zones (when the algorithm runs out of spaxels to bin). Typically, we obtain 1000 zones per galaxy.


starlight is a spectral synthesis method which combines the spectra from a base in order to match an observed spectrum . The search for an optimal model and the coefficients of the population vector associated with the base elements also allows for reddening, a velocity off-set and a stellar velocity dispersion . The code was first introduced in Cid Fernandes et al. (2004) and substantially optimized by Cid Fernandes et al. (2005). The code, didactic introductions and a user manual are available at We are actually using a new and yet to be documented version of starlight, but since most of the new features are not really used in our analysis, the publicly available manual remains an accurate reference for the purposes of this paper.

The virtues and caveats of full spectral fits are discussed in a number of articles (e.g., Panter et al. 2003; Ocvirk et al. 2006; Cid Fernandes 2007; Tojeiro et al. 2007; McArthur, González & Courteau 2009; Sánchez-Blázquez et al. 2011). A general consensus in the field is that uncertainties in the results for individual objects average out for large statistical samples (Panter et al. 2007). With CALIFA, each galaxy is a statistical sample per se! Hence, even if the results for single spaxels (or zones) are not iron-clad, the overal trends should be robust (see Paper II for a detailed evaluation of the uncertainties in our starlight-based analysis of CALIFA data.) It is also worth pointing out that despite the diversity of spectral synthesis methods, substantial changes on the results are more likely to come from revisions in the input data (as happened, for instance, with the recalibration of SDSS spectra between data releases 5 and 6) and from updates in the base models, the single most important ingredient in any spectral synthesis analysis.

The results reported in this paper rely on a spectral base of Simple Stellar Populations (SSPs) comprising 4 metallicities, , 0.4, 1 and , and 39 ages between and yr. This base (dubbed ”GM” in Paper II) combines the Granada models of González Delgado et al. (2005) for Myr with those of Vazdekis et al. (2010, as updated by Falcón-Barroso et al. 2011) for larger ages. They are based on the Salpeter Initial Mass Function and the evolutionary tracks by Girardi et al. (2000), except for the youngest ages ( Myr), which are based on Geneva tracks (Schaller et al. 1992; Schaerer et al. 1993a,b; Charbonnel et al. 1993). A comparison of results obtained with other bases is presented in Paper II, to which the reader is also referred for examples of the spectral fits.

The fits were carried out in the 3800–6850 Å interval. Reddening was modeled with the Cardelli, Clayton & Mathis (1989) curve, wit . As in Cid Fernandes et al. (2005), the main emission lines as well as the NaI D doublet (because of its sensitivity to ISM absorption) were masked. All runs were performed in the IAA-GRID, a network of computers which allows us to process spectra in less than a day.

As documented in the user manual, starlight outputs a large array of quantities. In broad terms, these can be categorized as

  1. Input data: File names, configuration options, and other information provided either by the user explicitly or derived from the user-provided spectrum. Base-related data are also reported, including ages and metallicities of each component, as well as the corresponding light-to-mass ratios (at the chosen normalization wavelength, ) and a returned-mass correction factor.

  2. Fit results: Figures of merit (, ), kinematic parameters (, ), the V-band extinction (), total stellar masses and luminosities, and population vectors, expressed in terms of light () and mass () fractions.

  3. Spectral data: The input (), best model fit ), and weight spectra (, which equals except for flagged, masked and clipped pixels).

All these data are stored in a plain ASCII file, one for each fitted spectrum. The population vectors can be post-processed in a number of ways to obtain SFHs, to compute the mass growth as a function of time, to condense the results into first moments like the mean and , etc. This scheme works fine for studies of independent objects, like SDSS galaxies or globular clusters.

IFS data, however, require a more structured level of organization, as it is cumbersome and inefficient to handle so many files separately. To handle datacubes we developed pycasso, the Python Califa Starlight Synthesis Organizer. pycasso comprises three main parts:

  1. A writer module, which packs the output of all starlight fits of the individual zones of a same galaxy into a single FITS (or HDF5) file, which also contains all information generated during the pre-processing steps (Section 3), as well as information propagated from the original datacube.

  2. A reader module, which reads this file and structures all the data in an easy to access and manipulate format.

  3. A post-processing module, which performs a series of common operations such as mapping any property from zones to pixels, resampling and smoothing the population vectors, computing growth functions, averaging in spatial, time or metallicity dimensions, etc. All these functionalities are conveniently wrapped as python packages, easily imported into the user’s own analysis code.

pycasso will be the main work-horse behind a series of articles by our collaboration, so we dedicate the remainder of this paper to illustrating its main products.

5 Results

The Sb galaxy NGC 2916 (CALIFA 277) was chosen as a show-case. For the adopted distance of 56 Mpc, ( CALIFA spaxel) corresponds to 270 pc. With an -band absolute magnitude of , and colors, and concentration index , NGC 2916 can be characterized as a blue-disk galaxy with an intermediate mass (). Rudnick, Rix & Kennicutt (2000) have shown that the galaxy show some degree of lopsidedness, indicating a weak interaction (possibly with its irregular satellite away, Gutiérrez, Azzaro & Prada 2002) or a minor merger. Moorthy & Holtzman (2006) have obtained long-slit spectra and measured emission lines and Lick indices along the position angle PA . They found that the nuclear spectrum is well within the AGN wing of the BPT diagram [Oiii]/H vs. [Nii]/H, in a location close to what is now accepted as borderline between Seyferts and LINERs (Kewley et al. 2006).

Our Voronoi binning of the CALIFA datacube produces zones. The corresponding spectra (the corresponding to the spatially integrated spectrum) were fitted with starlight and the results feed into pycasso. The following subsections describe how we handled the data and the products of the analysis.

5.1 Fit quality assessment

Figure 3: Maps of spectral-fit quality indicators. (a) Mean relative deviation , (b) per fitted flux, (c) number of fluxes used in the fit, and (d) number of clipped fluxes.

As with 1D spectra, the quality of the spectral fits should be examined before proceeding to interpretation of the results. This subsection describes some standard methods to do this, the caveats involved and strategies to circumvent them.

Fig. 3a shows the map of the fit quality indicator


where and are the observed and model spectra, respectively, and the sum is carried over the wavelengths actually used in the fit, ie, discarding masked, flagged and clipped pixels.4 For this galaxy spans the range between 1.0 and 5.6%, with a median value of 3.0%. increases towards the outer regions, where the is smaller (Fig. 1c). The map (Fig. 3b) shows the opposite behavior, since the errors also increase outwards. Because of its non-explicit reliance on the (often hard to compute, if not altogether unavailable) spectrum, as well as because of its easily grasped meaning, is a more useful figure of merit to assess fit quality.

Inspection of the highest spectra often reveal non-masked emission lines or artifacts, like imperfect masking of foreground sources, slightly misplaced flags, or values which should have been clipped by starlight but were not because of large . We emphasize that such bugs are very rare in CALIFA. The median for the zone spectra analysed so far is just 4% (corresponding to an equivalent of 25), and in less than 2% of the cases exceeds 10%. This high rate of success is mainly due to the carefully derived errors and flags in the reduction pipeline and further refined in our pre-processing steps (Section 3).

One can use such maps to define spatial masks over which the spectral fits satisfy some quality threshold. Experiments with the whole data set suggest that a reasonable quality-control limit should be in the to 10% range. Not surprisingly, the large residuals tend to be located in the outer regions of the galaxies. Generally speaking, results beyond 2–3 HLR should be interpreted with greater care. No -based cut is needed for our example galaxy. For the benefit of starlight users, it is nevertheless worth making some further remarks on quality control.

Subtleties and caveats on quality control

In full spectral fitting methods, seemingly trivial operations like imposing a fit-quality threshold often hide subtleties which can easily go unnoticed, particularly when processing tons of data. The paragraphs below (focused on starlight but extensible to other codes) review some of these.

First, one should distinguish cases where is large because of bad data from those where large residuals arise because of the user’s failure to properly inform, through spectral masks () or flags (), spectral regions which should be ignored in the fit. For instance, one can have an excellent spectrum of an HII region with several strong emission lines besides those included in a generic emission line mask feed into starlight, causing an artificially large . Large velocity offsets can have a similar effect, as emission lines get shifted out of fixed windows. Similarly, sky residuals and other artefacts missed out by the flags lead to large even when the overall spectrum is good. In short, not everything which fails a blind quality control is actually bad.

The clipping options implemented in starlight help spotting pixels which are too hard to fit and thus probably represent non-stellar or spurious features. Our fits use the “NSIGMA” clipping method and a conservative threshold, meaning that we only clip pixels when . Fig. 3c maps the number of clipped points in CALIFA 277. One sees a peak in the central pixels, where the is so high that our clipping does not forgive even small deviations (most often associated with problems with the base models rather than with the data). Elsewhere, very few pixels were clipped. The point to highlight here is that clipping only worked because we have a reliable . Other clipping methods can be used when this is not the case, but in all cases the user should always check carefully the output, as it may easily happen that too many points are clipped.

Custom-made spectral masks also help. For instance, since Mateus et al. (2006), the starlight analysis of SDSS spectra employs individual masks constructed by searching for emission lines in the residual spectrum obtained from a first fit, and taking into consideration the local noise level. is then refitted with this tailor-made , circumventing some of the issues above. This refinement has not yet been implemented for CALIFA.

Errors, flags and masks are, of course, secondary actors in any spectral synthesis analysis, but these general remarks illustrate that it pays off to dedicate them some attention.

5.2 2D maps: Stellar light, mass and extinction

Figure 4: (a) Synthesitc surface brightness at 5635 Å, (b) stellar extinction map (), (c) dereddened image, (d) stellar mass surface density.

One of the main products of spectral synthesis is to convert light to mass. Fig. 4 illustrates the results for CALIFA 277. Its top-left panel shows the surface density of the luminosity at Å, while the top-right panel show the derived extinction map. The dust corrected image and the mass surface density are also shown.

The effects of the Voronoi binning are present in all panels, but are much more salient in the map. This is because the light and mass images were “dezonified” by scaling the value at each spaxel by its fractional contribution to the total flux in a zone (). For instance, for the mass surface density () this operation reads


where () denotes the area in a spaxel (zone), and


with as the mean flux in the Å region, the same used in the Voronoi zoning (Section 3). This operation was applied to luminosity and mass related quantities, producing somewhat smoother images than obtained with . Intensive properties (like , mean ages, and etc.), however, cannot be dezonified.

The stellar extinction map shows low values of (of order 0.1–0.2 mag), with slight enhancements in the nuclear region and the arms. Overall, however, there is relatively little variation across the face of the galaxy (see also Fig. 10). The light and mass images have a similar structure, both showing the bulge at HLR, and the disc beyond that. The spiral arms are less prominent in terms of mass than in light because of the higher of stars in the arms.

The total stellar mass obtained from the sum of the zone masses is . This is the mass locked in stars nowadays. Counting also the mass returned by stars to the interstellar medium, were involved in star formation.These values ignore the mass in the masked region around the foreground star northeast of the nucleus. pycasso can fill in such holes with values estimated from (circular or elliptical) radial profiles. For CALIFA 277, this correction increases the masses quote above by 4%.

5.3 2D maps: Kinematics

Figure 5: Kinematical products of the spectral fitting. The values are not corrected for instrumental broadening, which dominates below 140 km/s. The horizontal stripes in the position-velocity and panels correspond to Voronoi zones.

Fig. 5 shows the and fields estimated by starlight, as well as a position-velocity diagram and a radially averaged profile. The field indicates a projected rotation velocity of km/s. The V500 resolution prevents the determination of values below km/s, and the flat profile beyond HLR reflects this resolution limit (notice that we have not corrected the values for the instrumental resolution in this plot). Only in the inner regions the values are trustable, reaching 200 km/s in the nucleus (equivalent to 140 km/s after correcting for instrumental broadening).

The kinematical information derived from our analysis will be superseded by studies based on the higher spectral resolution V1200 datacubes (Falcón-Barroso et al, in prep). Eventually, one can envisage feeding the parameters derived from these more precise analysis back into the starlight fits (using its fixed-kinematics mode). In fact, this feedback may well turn out to improve our stellar population analysis, given the potential degeneracy between and (decreasing the former while increasing the latter have the same global effect of making metal lines deeper; Koleva et al 2009; Sánchez-Blázquez et al. 2011).

5.4 2D maps: Mean ages and metallicities

Figure 6: Luminosity (top) and mass (bottom) weighted mean age (left) and metallicity (right) maps for CALIFA 277.

The crudest way to quantify the SFH of a system is to compress the age and metallicity distributions encoded in the population vectors to their first moments. For this purpose we will use the following definitions:


where is the fraction of light at due to the base population with age and metallicity . The mass weighted versions of these indicators, and , are obtained replacing by the mass fraction .

Fig. 6 shows the light and mass weighted mean (log) age and metallicity maps. The image shows a steady increase towards the center. Outside 1 HLR, traces of the spiral arms are noticeable as regions of lower age (as in the SDSS color image in Fig. 1a, the arms are brighter the western half of the image). Because of the large weight of old populations, spans a smaller dynamical range than , hence producing lower contrast maps, but the outwardly decreasing age is still visible. Negative gradients are also clearly present in metallicity, with indications of flattening within the bulge region. As in other 2D maps, the small scale fluctuations towards the edges are at least in part due to the lower .

Figure 7: Mean age versus mean metallicity versus extinction, coding zones by their distance to the nucleus. Green diamonds: HLR. Red triangles: 0.5–1 HLR. Blue stars: 1–2 HLR. Black dots: HLR.

The fact that and grow in the same sense suggests that the results are not badly affected by the age-metallicity degeneracy. As shown by Sánchez-Blázquez et al. (2011), full spectral synthesis is much less sensitive to this than conventional index-based approaches. Fig. 7 plots extinction versus mean age versus metallicity, the main properties involved in spectral synthesis. Points are color and symbol coded by their distance from the nucleus. One sees that the highest values of are found for old central ( HLR) and young outer ( HLR) populations, and that bears no correlation with . The middle panel shows the positive - relation inferred from the 2D maps. Within 0.5 HLR (green diamonds), however, anticorrelates with , most likely due to the age-metallicity degeneracy. The mean values in this central region straddle 1 and 1.5 , the two highest metallicities in our base.

5.5 2D maps: , ,

Figure 8: Maps of the percentage contribution of Young ( Gyr), Intermediate age (0.14–1.4 Gyr), and Old ( Gyr) stars to the observed light at 5635 Å.

Population synthesis studies in the past found that a useful way to summarize the SFH is to condense the age distribution encoded in the population vector into age ranges. This strategy comes from a time when the analysis was based on equivalent widths and colors (Bica 1988; Bica et al. 1994; Cid Fernandes et al 2001, 2003), but it was also applied to full spectral fits (Gonzalez Delgado et al. 2004).

Fig. 8 presents images of the light fraction due to Young, Intermediate and Old populations (, and ), defined as those with , , and Gyr, respectively. As usual, the choice of borderlines is somewhat subjective, constrained only by the underlying idea of grouping base elements covering relatively wide age ranges. For the base used in our starlight runs the Y, I and O bins contain , and SSPs, respectively (where the comes from the 4 metallicities).

The plots show that the spiral arms stand out more clearly in the map, specially in the western half of the galaxy, as also seen in the SDSS color composite (Fig. 1a). Very little of the light from HLR comes from young stars. Intermediate age populations do contribute more, but the central HLR is completely dominated by old stars. Overall, however, despite the added informational content, these maps do not visibly add much to the image (Fig. 6).

5.6 2D maps: Star formation rates and IFS-based variations over Scalo’s parameter

Figure 9: Spatially resolved SFR surface density (eq. 11) on the last (top), 142 (middle) and 1420 (bottom) Myr, in units of different reference values. The left panels compare to , i.e, the all-times average at , thus providing a local version of Scalo’s parameter (eq. 12). Middle panels compares the local to the all-times average over the whole galaxy (eq. 13). Panels on the right compare to that of the galaxy as a whole over the same time scale (eq. 14). All images are on log scale from to , such that only values above the corresponding reference value are visible.

The base used in our fits comprises instantaneous bursts (ie., SSPs), so, despite the large number of ages considered, our SFHs are not continuous and hence not derivable. There are, nonetheless, ways of defining star formation rates (SFR).

The simplest way to estimate a SFR is to cumulate all the stellar mass formed since a lookback time of and perform a mass-over-time average:


where is the mass (at and per unit area) of stars formed at lookback time .5 Eq. 11 gives the mean SFR surface density since . One can tune to reach different depths in the past, but as increases this estimator becomes increasingly useless, converging to the mass density divided by Gyr (the largest age in the base).

It is often more useful to consider SFRs in relation to some fiducial value, instead of absolute units. The classical example is Scalo’s birthrate parameter, , which measures the SFR in the recent past () with respect to its average over the whole lifetime of the system.6 This is trivially obtained dividing by its asymptotic value


The superscript “loc” is to emphasize that the reference lifetime average SFR is that of same spatial location . One should recall that while young stars in a spaxel were probably born there (or near it), old ones may have migrated from different regions. Thus, despite the identical definitions, the physical meaning of is not really the same as for galaxies as a whole.

IFS data allow for other definitions of . For instance, one might prefer to compare to the of the galaxy as a whole, the bulge, the disc or any general region. This variation measures the “present” () “here” (at ) to the “past” () in a spatial region :


A formally similar but conceptually different definition is obtained using for the reference value in the denominator the SFR surface density within the same time-scale but evaluated in a different region:


i.e., to compare not “now versus then”, but “now here versus now there”. Together with equations 12 and 13, this definition cover the different combinations of time and space enabled by the application of fossil methods to IFS data.

Fig. 9 shows maps of (left panels), (middle), and (right) for three values of : 14 (top), 142 (middle) and 1420 Myr (bottom). In all panels the color scale is deliberately saturated to highlight regions where is larger then chosen reference value—the dynamical range of the images goes from 1 to 5 (0 to 0.7 in log) in the corresponding relative units.

Panel a shows that spaxels which have formed stars over the past 14 Myr at a larger rate than their respective past average are all located in the outer regions. Nowhere within HLR one finds . Considering the past 142 Myr (panel d), the inner ”deficit” covers as much as 1.2 HLR. On the longest time scale considered (g) one again sees that the inner spaxels have been less active than their life-long average.

Maps look considerably different in the middle column of panels, where the local is normalized to the past average over the whole galaxy (eq. 13). The spiral arms of CALIFA 277 appear better delineated in (panel b) than in any map. Interestingly, the map is practically featureless for Myr (e). This suggests that star-formation is not continuous over this time scale, but instead happens in a bursty mode. Over the past 1.4 Gyr (h), one again sees traces of the arms, but with a lower amplitude than in panel b, qualitatively consistent with the cumulative effect of an intermittent sequence of short duration bursts. (Due to the logarithmic age resolution of fossil methods, such short bursts can only be recognized as such at the very young ages sampled in panel b.)

The right column panels show our novel relative SFR index (eq. 14), using the whole galaxy as reference region. Unlike Scalo’s , this index does not compare present to past, but present to present elsewhere. Its 2D maps can be understood as “snapshots” of the SFR in the galaxy taken with an “exposure time” . Keeping the “diaphragm” open for only the last 14 Myr (panel c) highlights the ongoing star-formation in the galaxy. The spiral arms are clearly visible, and the structures are more focused than in the map (panel b). Integrating for 142 Myr (f), the inner parts of the arms fade, but their outer ( HLR) parts brighten up in terms of . Over the past 142 Myr, these regions have formed more stars than anywhere else in the galaxy, even though comparing with panel e we find that only parts of the regions are forming stars at a rate per unit area larger than that of the galaxy as a whole through its entire life (ie, ). For Gyr (panel i) the image becomes more blurred. Because of the long time scale, the map is nearly indistinguishable from that obtained with for the same (panel h).

One thus sees that, despite some degree of redundancy (the numerator is always the same), these definitions offer different and complementary views of the star formation in a galaxy.

5.7 1D spatial maps: radial profiles

Figure 10: Radial profiles of several properties. Grey points correspond to values in individual spaxels. The green circles and solid lines represent the average of the plotted quantity, while the error bars represent the dispersion in the -bin. The dashed magenta horizontal line (only in the plots not involving surface densities) marks the value derived from the starlight analysis of the integrated spectrum, i.e. collapsing the dimensions of the datacube. The stars in the bottom panels mark the 4 metallicities in the base models.

The information contained in 2D maps like the ones shown above can sometimes be hard to absorb. Azimuthal averaging is a useful way to compress and help digesting 2D data. pycasso provides both circular and elliptical -to- conversion tensors. The examples below are based on a circular mapping.

For quantities like , one can think of two types of radial averaging. The first is to add up all the stellar mass in spaxels within a given -bin and then divide by the bin area (counting only non-empty spaxels). The second is to average the surface density values for each spaxel directly. The same applies to, for instance, (eq. 9): One can either add up the value of product for all spaxels at a given and then divide by the corresponding , or else simply average the values for each spaxel in the -bin. The first method, which we may call area averaging, effectively collapses the galaxy to 1D, but this kind of averaging cannot be applied to quantities which do not involve light or mass, like and . For the sake of uniformity, in this subsection we chose the second type of radial averaging (henceforth spaxel averaging), noting that the two methods almost always give nearly identical results.

Fig. 10 shows several pycasso products as a function of . The grey points represent values for individual spaxels. The solid line and green circles show the mean value of the plotted quantity among the spaxels in the same bin, and the error bars map the corresponding dispersion. (As discussed in Paper II, actual error bars in radial profiles are much smaller due to the large number statistics.) Several of the features discussed while describing the 2D images are clearly seen in these plots. The and gradients, in particular, are cleanly depicted in these plots. As pointed out in Section 5.4, means ages and metallicities increase simultaneously towards the nucleus, while does not change much, indicating that the classical degeneracies in stellar population studies are not playing a strong role in shaping the results, with the possible exception of the central regions. As in Fig. 5, horizontal stripes of gray points in panels c, e, f, g and h of Fig. 10 correspond to spaxels in the same Voronoi zones. The stripes disappear in panels a, b and d because these extensive quantities were ”dezonified” (cf. eq. 7).

5.8 1D temporal maps: evolutionary curves

Figure 11: Time evolution of light and mass representations of the SFH, as derived from the starlight fits. The discrete and ”population vectors” involved in these curves were smoothed in , marginalized in and integrated over the spatial regions indicated. Dotted grey lines indicate the nucleus (central pixel). Green, red, and blue lines correspond to HLR, HLR and HLR regions, respectively. The black line represents the whole galaxy evolution as reconstructed from the sum of its parts, while (as in Fig. 10) the magenta dashed line is used to represent the results obtained for the spatially integrated spectrum.

The main power of fossil record methods is to open the time domain, allowing evolutionary studies. The caveat is that it does so with the typical logarithmic resolution inherent to stellar evolution, but several studies show that much can be learned about galaxy evolution despite this fundamental constraint. Clearly, though, our 39 ages base is highly overdimensioned, so some sort of SFH compression is necessary. As discussed by Asari et al. (2007), age resolutions between 0.5 and 1 dex are reasonable. We thus smooth the light and mass population vectors with a gaussian of FWHM dex in .

Fig. 11 plots a number of our synthesis products against age. To keep an at least partial representation of the spatial information, we plot evolutionary curves derived for different radial regions: The nucleus, defined as the central pixel (plotted in grey dashed lines), (dashed green), (solid red), and HLR (solid blue).

The top panels show the growth of luminosity (Fig. 11a) and mass (b).7 Among other things, this plot illustrates the difference between light and mass, in the sense that regions inside and outside HLR do not have the same mass, despite having (by definition) the same light. While HLR kpc, the radius containing half of the mass is 3.4 kpc. This is a direct consequence of the stellar population gradients in this galaxy. González Delgado et al. (in prep) investigates the relation between light and mass based radii for different types of galaxies.

It is also seen that both light and mass grow at different speeds for different regions. This is better appreciated in panels c and d, were each growth curve is plotted on a 0 to 1 scale, with 1 representing the present values (which tantamounts to cumulating the equivalent and vectors for each region). The nucleus reached 80% of its mass at Gyr, while the HLR region did so later, at Gyr. As shown by Pérez et al. (2013), this inside-out ordering of the mass assembly history applies to essentially all massive galaxies.

Fig. 11e shows . The SFR per unit area decreases from the nucleus outwards through most of the time, but the trend is reversed in the last Myr. Notice that, because of the smoothing, is now a continuous function of time, so that this panel represent the instantaneous SFR, as opposed to the mass-over-time definition in eq. 11. Fig. 11f shows Scalos’s parameter for each region, which does use the running mean of eq. 11. This plot helps interpreting Fig. 9, which opens up the radial regions into full 2D maps, but compress the time axis by stipulating fixed time scales.

5.9 The whole versus the sum of the parts

One of the applications of CALIFA is to assess the effects of the lack of spatially resolved information on the derivation of physical properties and SFHs from integrated-spectra surveys. The horizontal dashed magenta lines in Fig. 10 illustrate this point. They represent the results obtained from the analysis of the spectrum of datacube as a whole, ie., adding upp all spaxels to emulate an integrated spectrum. Qualitatively, one expects that this should produce properties typical of the galaxy zones as a whole. This expectation is borne out in Fig. 10, where one sees that the extinction, mean ages and metallicities marked by the dashed magenta line do represent an overall average of the spatially resolved values. In this particular example, they all match quite well the values at 1 HLR. The stellar masses obtained from the whole and the sum of the parts are also in excellent agreement, differening by just 1%. Similarly, galaxy-wide luminosity and mass weighted mean ages and metallicities computed in these two ways agree to within 0.05 dex.

But can one derive the SFH of a galaxy out of an integrated spectrum? Fig. 11 says that, at least in CALIFA 277, the answer is yes. The dashed magenta and the black solid lines are very similar in all panels of this figure. As before, the former represents the results obtained from the analysis of the spatially compressed datacube, while the black curves are computed adding the starlight results for each spaxel. We are thus comparing the whole against the sum of the parts, and they match.

At first sight, this may seem a trivial result, but it is not. Formally, one only expects this to happen if is the same at all positions, a condition which is approximately meet in CALIFA 277. In objects where dust has strong spatial variations, the spectral fitting of the global spectrum with a single will inevitably operate compensations, like increasing the age to compensate for an understimated and vice-versa. Also, the combined effects of population gradients and kinematics (e.g., disc population contributing more to the wings of absorption lines than the ones from the bulge) are hard to predict. Since we are presenting results for a single example galaxy, it remains to be seen how general this “coincidence” is.

5.10 Space time diagrams: SFHs in 2D

Figure 12: - diagrams showing the radially averaged distribution of light, mass and SFR as a function distance from the nucleus and age. Left: Luminosity at Å per unit area. . Middle: Stellar mass formed per unit area . Right: Time dependent SFR per unit area . The solid (white) lines represent the sum over all spaxels for a given age.

As is evident at this point of the paper, one of the main challenges involved in the analysis of the multi dimensional data built from the combination of the spatial dimensions with the and dimensions opened by population synthesis is how to visualize the results. All examples shown so far project (or average over) two or more of these axes.

Fig. 12 shows an attempt to visualize galaxy evolution as a function of both space and time. The trick is to compress into , and collapse the axis, producing a radially averaged SFH map. The left panel shows the luminosity density at each radial position and for each age . As for the other panels in this figure, the original array from which this map was derived was smoothed in , marginalized over and collapsed into one spatial dimension. Unlike in Fig. 10, we now use the area averaging method discussed in Section 5.7, but the results do not depend strongly on this choice.

The diagram shows that the light from old stars is concentrated in the bulge, while the youngest stars shine equally at all . Also, Gyr stars are more smoothly distributed in ages in the disc than in the bulge, where they are concentrated in two populations, one at Gyr and the other older than 10 Gyr. Fig. 12b shows the formed stellar mass surface density as a function of and . As usual, due to the highly non-linear relation between light and mass in stars, the young populations which show up so well in light practically disappear when seen in mass. Young stars reappear in the right panel, which shows the evolution in time and space of the SFR surface density. As in Fig. 11g, these are obtained from the instantaneous SFR at each spaxel, computed as in Asari et al. (2007):


The age sampling is logarithmic ( constant), so that the SFR is proportional to the stellar mass formed at divided by . The and diagrams therefore carry the same information; it is the factor which makes them look so different.

When interpreting these or any other plot involving age, one should always keep in mind the logarithmic age resolution. For instance, taken at face value, the spatially integrated SFR of CALIFA 277 in the last few Myr is almost equal to its historical peak, Gyr ago (solid line in Fig. 12c). However, these two peaks span vastly different time intervals, roughly by a factor of Gyr/Myr = a thousand. The peak SFR around 2.5 Gyr ago, and indeed even well before that (at Gyr, when most of the mass was formed), were surely higher than the one we are now seeing at a few Myr, but cannot be resolved in time. No fossil record method will ever be able to distinguish bursts much shorter than their current age. These limitations are well known in the field, but it is fit to recall them to avoid missinterpretation of the results.

Notwithstanding such age resolution issues, Fig. 12 represents a novel way of looking at galaxies both in time and space, and it shows clear evidence in favor of an inside out growth scenario, in which the outer regions assemble their mass at a slower pace and a over more extended period. This behavior is not unique to this galaxy. Pérez et al (2012) find just the same in a study of the first 105 galaxies observed by CALIFA.

5.11 Space vs. time ”snapshots”

Figure 13: ”Pseudo snapshots” of the cumulative stellar mass surface density of CALIFA 277 at different ages: From [yr] to 8 (as labeled in each panel).

There is so much one can plot in a single image. To go beyond the - diagrams of Fig. 12 one needs more dimensions than a sheet of paper can accommodate. An alternative is shown in Fig. 13, which shows -slices through the cumulative cube. Unlike the previous diagrams, which compressed the SFH in one way or another, these plots reveal the full 2D nature of the mass assembly history of CALIFA 277.

Age sequences like those in Fig. 13 can be constructed for several, but not all properties inferred though starlight or any other fossil method. One can, for instance, replace mass, by intrinsic luminosity, SFR or stellar metallicity, in cumulative or differential form, but it is obviously impossible to reconstruct maps of , , as a function of . It is also worth pointing out that the panels in Fig. 13 are not truly snapshots of a movie, ie., they are not pictures of the galaxy as it appeared at different look-back times. Instead, these are maps of where stars of a given age are located nowadays.

Simulators should take notice of these simple facts. IFS data plus fossil methods provide a rich, yet inevitably limited form of time-travel. Illustrative and beautiful as they are, movies of stars and gas particles moving through time and space will never be directly compared to anything observational. In other words, simulations should be convolved through this “reality filter”. The observationally relevant predictions are the distribution of stellar ages and metallicities as a function of .

5.12 Emission line maps

Figure 14: Residual spectral cube in a 3D rendering.

Independent of the stellar population applications described throughout this paper, the starlight fits provide a residual spectrum which is, inasmuch as possible, free of stellar light. This is of great aid in the measurement of emission lines, and hence in the derivation of a series of diagnostics of the warm ISM properties, such as its kinematics, nebular extinction and chemical abundance. The complementarity of nebular and stellar analysis has been amply explored in studies of integrated spectra (eg, Cid Fernandes et al. 2007; Asari et al. 2007; Pérez-Montero et al. 2010; García-Benito & Pérez-Montero 2012), and the addition of two spatial dimensions to this analysis holds great promise. Kehrig et al. (2012), for instance, combine starlight analysis and nebular emission to study the impact of post-AGB populations on the ionization of the ISM in early type galaxies.

Fig. 14 shows a residuals cube, with some of the main emission lines marked. The plot shows that emission lines are brighter towards the outside, in particular over the spiral arms, as could be guessed from the distribution of actively star-forming regions traced by our SFH analysis (e.g., Fig. 9). Despite its active nucleus, the nuclear line fluxes are not particularly bright, implying a low luminosity AGN.

6 Summary

This paper described a complete ”analysis pipeline” which processes reduced IFS data through the starlight spectral synthesis code. Data for the spiral galaxy NGC 2916 (CALIFA 277) were used to illustrate the journey from datacubes to multi-dimensional maps of the synthesis products, focusing on the technical and methodological aspects. The three main steps in the analysis can be summarized as follows:

  1. Pre-processing: This step comprises (a) construction of a spatial mask; (b) adjustments in the flag and error spectra; (c) rest-framing and resampling; (d) spatial binning. The last three stages were coded in a generic python package (qbick) which can run in a fully automated fashion. Spatial binning was accomplished with an Voronoi tesselation scheme, tuned to produce zones with spectra. In practice, because of the good quality of the data, no binning was necessary within HLR. The issue of correlated errors in spatial binning was discussed and an empirical recipe to correct for it was presented. The output of these pre-processing steps are datacubes with fluxes, errors and flags ready for a detailed -by- spectral analysis, as well as spatial masks, a zone map and other ancillary products.

  2. starlight fits and organization of the synthesis results: Spectral fits for all spatial zones were performed using a base of SSP spectra from a combination of MILES and Granada models, spanning ages from 1 Myr to 14 Gyr and 4 metallicities between 0.2 and 1.5 solar. The results were packed in a coherent multi-layered FITS (or HDF5) file with the Python Califa Starlight Synthesis Organizer, pycasso.

  3. Post-processing tools: pycasso includes a reader module and a series of analysis tools to perform operations like zone-to-pixel image reconstruction, mapping the spectral and stellar population properties derived by starlight into multi-dimensional time, metallicity, and spatial coordinates, averaging in spatial and temporal dimensions, manipulation of the stellar population arrays, etc.

Some of the products of this whole analysis are standard in IFS work, like the kinematical field and emission line maps. The real novelty, of course, resides in the spatially resolved stellar population products recovered from the fossil record of galaxy evolution encoded in the datacubes. The 2D products range from maps of the stellar extinction and surface densities of stellar mass to more elaborate products like mean ages and metallicities, time averaged star formation rates and -parameters. These quantities are all normally used in the analysis of integrated spectra, but their application to IFS data raises some conceptual issues, related to the unavoidable fact that not all stars in a given spaxel were born there. On the other hand, one can use spatially resolved data to construct diagnostics unaplicable to integrated data, such as indices comparing the present and past SFR in different regions, or comparing the “present here” with the “present elsewhere”. These IFS-based variations of the traditional parameter were shown to be useful to highlight different aspects of the spatially resolved SFH.

1D profiles in temporal or radial coordinates were used as a means to help interpreting the results. In our example galaxy, these compressed views of the manifold reveal clear age and metallicity negative gradients, as well as spatially dependent mass growth speeds compatible with an inside-out formation scenario. Finally, radius-age diagrams and “snapshot” sequences were introduced as further means of visualization.

Clearly, the application of spectral synthesis methods to IFS datacubes offers new ways of studying galaxies and their evolution. Future communications will use the tools laid out in this paper to explore the astrophysical implications of the results for galaxies in the CALIFA survey.

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. RCF thanks the hospitality of the IAA and the support of CAPES and CNPq. ALA 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 (SFS, PSB and JFB), 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.


  1. The detrended std. dev. is the rms variation with respect to a linear fit, useful to remove the effect of the spectral slope in the window.
  2. Note that this is an apparently superfluous definition, since we already have the noise spectrum. This ”empirical “ is nevertheless useful, as explained in Section 3.3.
  3. denotes the noise in spaxel , defined as the detrended rms in the same Å range used to define the signal for spatial binning purposes.
  4. Notice that eq. 6 is nearly, but not exactly identical to what is called “adev” in starlight’s output, the difference being the use of instead of in the denominator. This avoids faulty (but not flagged) entries with very low flux affecting the statistics. This rarely occurs in CALIFA, and when it does the pixels invariably have large errors, and so are irrelevant for the fit, but can still affect the value of .
  5. The stellar mass in this equation differs from that in eq. 7, Fig. 4, and the definitions of and . The latter are corrected by the mass returned to the medium by stars during their evolution (thus reflecting the mass currently in stars), while in the computation of SFR’s this correction should not be applied. We distinguish the two types of stellar masses with the prime superscript ( versus ).
  6. A related index is the so called specific SFR, which differs from just by the factor in equation 12.
  7. For clarity, the cumulative mass is computed with the total mass which has participated in star formation, instead of the one corrected for returned mass. In other words, we cumulate instead of .


  1. Asari, N.V., Cid Fernandes, R., Stasińska, G., Torres-Papaqui, J.P., Mateus, A., Sodré, L., Schoenell, W., Gomes, J. M., 2007, MNRAS, 381, 263
  2. Bell, E.F., Wolf, C., Meisenheimer, K., Rix, Hans-Walter, Borch, A., Dye, S., Kleinheinrich, M., Wisotzki, L., McIntosh, D.H. , 2004, ApJ, 608, 752
  3. Bertin, E., Arnouts, S., 1996, A&AS, 117, 393
  4. Bica, E. 1988, A&A, 195, 76
  5. Bica, E., Alloin, D., & Schmitt, H. R. 1994, A&A, 283, 805
  6. Blanc, G. A., Gebhardt, K., Heiderman, A., et al. 2010, New Horizons in Astronomy: Frank N. Bash Symposium 2009, 432, 180
  7. Bruzual G., Charlot S., 2003, MNRAS, 344, 1000
  8. Cappellari, M., & Copin, Y. , 2003, MNRAS, 342, 345
  9. Cardelli, J.A., Clayton, G.C., Mathis, J.S., 1989, ApJ, 345, 245
  10. Charbonnel, C., Meynet, G., Maeder, A., Schaeller, G., Schaerer, D., 1993, A&AS, 101, 415
  11. Cid Fernandes, R., Mateus, A., Sodré L., Stasińska, G., Gomes, J.M., 2005, MNRAS, 358, 363
  12. Cid Fernandes, R. 2007, IAU Symposium, 241, 461
  13. Falcón-Barroso, J., Sánchez-Blázquez, P., Vazdekis, A., Ricciardelli, E., Cardiel, N., Cenarro, A.J., Gorgas, J., Peletier, R. F., 2011, A&A, 532, 95
  14. Folkes, S., Ronen, S., Price, I., Lahav, O., et al., 1999, MNRAS, 308, 459
  15. García-Benito, R., & Pérez-Montero, E. 2012, MNRAS, 423, 406
  16. González Delgado , R. M.; Cerviño, M., Martins, L. P., Leitherer, C.; Hauschildt, P. H., 2005, MNRAS, 357, 945
  17. Gutiérrez, C. M., Azzaro, M., Prada, F., 2002, ApJS, 141, 61
  18. Ilbert, O., Capak, P., Salvato, M., et al., 2009, ApJ, 690, 1236
  19. Kehrig, C., Monreal-Ibero, A., Papaderos, P., et al. 2012, A&A, 540, A11
  20. Kelz, A., Verheijen, M. A. W., Roth, M. M., et al. 2006, PASP, 118, 129
  21. Koleva M., Prugniel P., De Rijcke S., 2009, Astron. Nachr. /AN 999, No. 88, 1
  22. Lilly, S. J., Le Fevre, O., Hammer, F., Crampton, D., 1996, ApJ, 460, 1L
  23. MacArthur, L. A., González, J. J., & Courteau, S. 2009, MNRAS, 395, 28
  24. Madau, P., Ferguson, H.C., Dickinson, M.E., Giavalisco, M., Steidel, C.C., Fruchter, A., 1996, MNRAS, 283, 1388
  25. Moles, M., Benítez, N., Aguerri, J. A. L. et al., 2008, AJ, 136, 1325
  26. Moorthy, B.K., Holtzman, J.A., 2006, MNRAS, 371, 538
  27. Ocvirk, P., Pichon, C., Lançon, A., & Thiébaut, E. 2006, MNRAS, 365, 46
  28. Panter, B., Heavens, A. F., & Jimenez, R. 2003, MNRAS, 343, 1145
  29. Panter, B., Jimenez, R., Heavens, A. F., & Charlot, S. 2007, MNRAS, 378, 1550
  30. Pérez, E., Cid Fernandes, R., González Delgado, R. M., et al. 2013, ApJ, 764, L1
  31. Pérez-González, P.G., Rieke, G.H., Villar, V., et al., 2008, ApJ, 675, 234
  32. Pérez-Montero, E., García-Benito, R., Hägele, G. F., & Díaz, Á. I. 2010, MNRAS, 404, 2037
  33. Roth, M. M., Kelz, A., Fechner, T., et al., 2005, PASP, 117, 62
  34. Rudnick, G., Rix, Hans-Walter, Kennicutt, R.C., Jr., 2000, ApJ, 538, 569
  35. Sánchez, S.F., Kennicutt, R.C., Gil de Paz, A., et al., 2012, A&A, 538, 8
  36. Sánchez-Blázquez, P., Ocvirk, P., Gibson, B. K., Pérez, I., & Peletier, R. F. 2011, MNRAS, 415, 709
  37. Schaerer, D., Charbonnel, C., Meynet, G., Maeder, A., Schaller, G., 1993a, A&A, 102, 339
  38. Schaerer, D., Meynet, G., Maeder, A., Schaller, G., 1993b, A&AS, 98, 523
  39. Schaller, G., Schaerer, D., Meynet, G., Maeder, A., 1992, A&AS, 96, 269
  40. Tojeiro, R., Heavens, A. F., Jimenez, R., & Panter, B. 2007, MNRAS, 381, 1252
  41. Vazdekis, A., Sńchez-Blázquez, P.; Falcón-Barroso, J.; Cenarro, A. J.; Beasley, M. A.; Cardiel, N.; Gorgas, J.; Peletier, R. F., 2010, MNRAS, 404, 1639
  42. Verheijen, M. A. W., Bershady, M. A., Andersen, D. R., et al., 2004, Astron. Nachr., 325, 151
  43. York, D.G., Adelman, J., Anderson, J.E., et al., 2000, AJ, 120, 1579
  44. Walcher, J., Groves, B., Budavári, T., & Dale, D. 2011, Ap&SS, 331, 1
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description