The stellar evolution of Luminous Red Galaxies, and its dependence on colour, redshift, luminosity and modelling

The motivation for this work is two fold. In Tojeiro & Percival (2010) we constrained the dynamical evolution of Sloan Digital Sky Survey (SDSS) LRGs by assuming the stellar evolution of these galaxies was known, and followed the passive evolution model (with a small metal-poor component) of Maraston et al. (2009) (M09). Given the evidence that ETGs have different histories according to their luminosity, as well as evidence for recent star formation in at least intermediate-mass ETGs, the assumption of passive evolution after a single star burst for all LRGs is clearly an oversimplification. The first of our goals is therefore to provide the community with empirically derived colour evolution models for the SDSS LRG sample that depend on luminosity, colour and redshift.

Secondly, our approach of using the fossil record to infer their star-formation and metallicity histories, as well as dust content, and from there compute their colour evolution as a function of redshift makes minimal prior assumptions about the evolution of the galaxies. The method adopted is to compute the star-formation and metallicity histories using VESPA (Tojeiro et al., 2007) - a full-spectrum fitting code that assumes no prior information about the star-formation or metallicity history of a galaxy other than it is non-negative. Non-parametric full-spectral fitting has advantages and disadvantages when compared to other methods. The main advantage comes from a well time-resolved and completely unconstrained star formation history, which does not suffer from any problems and biases traditionally associated with SSP-equivalent ages and analyses (e.g. Trager & Somerville 2009). On the other hand, full-spectral fitting can be more sensitive to inaccuracies in the modelling (particularly dust), and spectrophotometric calibrations. Parameter degeneracies are not a problem, in the sense that they can be estimated and incorporated into any analysis.

There are also clear advantages in using the fossil record to compute the colour evolution of a sample of galaxies. Primarily, one does away with the circular methodology of pre-selecting a sample across a range of redshifts that one believes is itself a uniform population of galaxies. This is traditionally done according to an evolutionary model and after applying the necessary K+e corrections and survey selection functions, and culling the observed sample to one that is self-consistent (i.e., any galaxy in this culled sample could have been observed at any redshift probed by the survey, according to a model). The obvious danger of such an approach is that the deductions about the evolution of the sample simply reflect the assumptions used to define it. In our approach, however, evolutionary colour paths that stray from the colour selection box are acceptable. We therefore de-couple the sample of galaxies observed at any given redshift from galaxies observed at earlier epochs, and we are insensitive to changing survey selections. Of course, if one wants to interpret this evolution in terms of the overall evolution of a population of galaxies then sample selection becomes crucially important. We leave that for a series of follow-up papers. In this work we introduce a new set of empirical models for the evolution of galaxies labelled as Luminous Red Galaxies, as a function of their observed colour, luminosity and redshift.

This paper is organised as follows. We start by describing the dataset used in this work in Section 2; in Section 3 we describe in detail the several steps of our methodology; in Section 4 we describe the three sets of SPS models used throughout the paper; we follow by presenting our results in Section 5, which we interpret in Section 6; and we finally summarise and conclude in Section 7.

We stack galaxies of similar colour, luminosity and at the same redshift. The resulting spectrum gives a weighted evolution of all the galaxies within the stack.

We work on a fixed grid in redshift, from 0.15 to 0.5 in steps of 0.05. At each redshift, we split the sample into four luminosity slices, such that each slice has approximately the same number of galaxies. Within each redshift bin and luminosity slice, we split the galaxies into stacks based on $r-i$ colour, with the number of stacked objects chosen such that each stack contains a fixed number of objects, set to be around 1000¹¹1For the sake of clarity, we will continue to refer to each region bound by colour, luminosity redshift as a cell, to each column of cells for the same redshift as a redshift bin, and to each grid of colour-redshift for a given luminosity range as a luminosity slice.. This number was chosen so that the SNR of the stacked spectra would permit differentiation between ages of interest in this paper (see Section 3.3 for details). The typical photometric errors increase steeply into the blue, so the spread in $g-r$ colour of LRGs would be expected to be larger than that in $r-i$ even if LRGs did not have any intrinsic scatter in their colours. There is also a linear relationship between redshift and $g-r$ up to $z\lesssim 0.35$, which we capture by binning in redshift. We therefore stack in $r-i$ colour rather than $g-r$. The result is a number of non-uniform grids (one per luminosity slice) in $r-i$ colour, depending on the density of objects in the colour-redshift plane. We choose to do this in favour of a uniform grid in order to keep a roughly constant SNR of the resulting stacked spectra. The disadvantage is that sparse regions in these diagrams need to be larger in order to reach the target number of galaxies, and variations in the intrinsic nature of the galaxies within each cell may be larger. This does not invalidate the resulting evolution model for this cell, which will simply be a roughly luminosity-weighted average evolution for all the galaxies within the cell, although it may mask larger departures from the mean within the cell. The above procedure occasionally results in cells that are smaller, in $r-i$, than the typical error on $r-i$ for individual galaxies. This means that there is likely scatter of galaxies across cells, and that the VESPA solutions from adjacent cells should only change smoothly (or as smoothly as the observational uncertainty in whatever quantity we bin). In other words the true resolution in $r-i$ is given by the photometric errors, and not by the size of the cells.

The goal of this paper is to compute a set of empirical evolutionary models for the colour of LRGs, and once again we assume that within a cell this is independent of luminosity. This assumption is strengthened by us taking more than one luminosity range, but there may be residual dependencies within a given chosen range. The purpose of the stacking procedure, therefore, is to compute an optimal estimate of the underlying spectrum, of which we assume each galaxy represents a random realisation. We discard luminosity information within each cell by normalising all spectra to a common wavelength, and we compute an inverse-variance weighted average, per wavelength point:

where $F$ is the stacked spectrum, $f_{\lambda,i}$ is the flux point at wavelength $\lambda$ of galaxy $i$ and $\sigma_{\lambda,i}$ the corresponding error. The sum is done over galaxies in any given cell.

The resulting error in each wavelength point of the stacked spectrum is simply

Note that by weighting the spectra within each stack based on signal-to-noise we are not modelling the total luminosity of the galaxies within the bin when we fit to the stacked spectrum. Instead our results should be interpreted as giving the history of a weighted function of the galaxies within that bin, although this weighting is actually close to a luminosity weighting, as spectral S/N is closely related to luminosity.

The typical number of galaxies per stack is 1000, and we have a total of 124 stacks.

We use VESPA to interpret the stacked spectrum of each cell in terms of a star formation and metallicity history. The algorithm is described in detail in Tojeiro et al. (2007), and it has been applied to SDSS’s LRG and Main Galaxy samples resulting in a queryable database, described in Tojeiro et al. (2009). We refer the reader to these two papers for details but, for completeness, we present now a very brief summary of the method.

In short, VESPA inverts the equation $F_{\lambda}=\hat{F}_{\lambda}$, where $F_{\lambda}$ is the observed rest-frame flux, and the model flux is

to solve for the star formation rate $\psi(t)$ (solar masses formed per unit of time) and the age $t$ and metallicity $Z$ of the components of the stellar populations that each give luminosity per unit wavelength $S_{\lambda}(t,Z)$, per unit mass. The dependency of the metallicity on age is unconstrained, turning this into a non-linear problem. To get around this, we interpolate between the tabulated values of $Z$ given by the SPS models (see Section 4) giving a piecewise linear behaviour:

where $S(t,Z_{a})$ and $S(t,Z_{b})$ are equivalent to $S_{\lambda}(t,Z)$ in Eq. (3), but at fixed metallicities $Z_{a}$ and $Z_{b}$, which bound the true $Z$.

Solving the problem then requires finding the correct metallicity range. One should not underestimate the complexity this implies - trying all possible combination of consecutive values of $Z_{a}$ and $Z_{b}$ in a grid of 16 age bins would lead to a total number of calculations of the order of $10^{9}$, which is unfeasible even in today’s fast personal workstations. We work around this problem using an iterative approach, which we describe in Tojeiro et al. (2007).

As for the treatment of dust, we follow the two-parameter dust model of Charlot & Fall (2000) in which young stars are embebbed in their birth cloud up to a time $t_{BC}$, when they break free into the inter-stellar medium (ISM):

where $\tau_{\lambda}^{ISM}$ is the optical depth of the ISM and $\tau_{\lambda}^{BC}$ is the optical depth of the birth cloud. In previous runs of VESPA, we took $t_{BC}=0.03$ Gyrs. However, in this run we have reduced resolution at the young end (see Section 3.3 for a discussion on the revised age grid) and we do not resolve star-formation under $74$ Myr. We therefore take this boundary as $t_{BC}$, and note that removing this dust component altogether has minimal effect in our results and conclusions.

There is a variety of choices for the form of $f_{dust}(\tau_{\lambda})$. To model the dust in the ISM, we use the mixed slab model of Charlot & Fall (2000) for low optical depths ($\tau_{V}\leq 1$), for which

where $E_{1}$ is the exponential integral and $\tau_{\lambda}$ is the optical depth of the slab. This model is known to be less accurate for high dust values, and for optical depths greater than one we take a uniform screening model with

We only use the uniform screening model to model the dust in the birth cloud and we use $\tau_{\lambda}=\tau_{V}(\lambda/5500\AA )^{-0.7}$ as our extinction curve for both environments.

Recovering the star formation history, and its dust content, is solved by making assumptions about the form of $f_{dust}$ and $S_{\lambda}(t,Z)$ and by minimizing

Even though the problem has an analytical solution, a dataset perturbed by noise or that is otherwise deteriorated leads to instabilities in the matrix inversion and the recovered solutions can be entirely dominated by noise (Ocvirk et al., 2006). VESPA has a self-regularization mechanism which estimates how many independent parameters one should recover from a given a dataset that has been perturbed. The result is a parametrization which varies from galaxy to galaxy, depending on its SNR and wavelength coverage.

It is the loss of resolution in lookback time, for low SNR spectra, that prompts the need for stacking in data-space. In other words, given the self-regularization mechanism and the inherent loss of - noisy! - information in cases of poor-quality data, it is not equivalent to stack in data or solution space. So, in general, the average of the star formation histories of the individual galaxies is not the same as the star formation history of the average spectrum. This becomes less and less true for increasingly better data of the individual objects, but in the case of the LRGs the data quality is sufficiently poor - especially at higher redshift - that stacking becomes necessary.

The physical quantities output by VESPA for each stacked spectrum are:

• •

  star formation fraction in time interval $\Delta t_{j}$, $x_{j}$;

• •

  mass-weighted metallicity for mass formed in time interval $\Delta t_{j}$, $Z_{j}$;

• •

  dust attenuation for stars younger than 0.074 Gyrs, $\tau_{BC}$;

• •

  dust attenuation for all stars, $\tau_{ISM}$.

Given our choice for stacking method and the lack of overall normalisation, we cannot convert the mass fractions, $x_{j}$, into an absolute mass recovered in each time interval. This, however, does not pose a problem as we aim to compute a set of colour evolution tracks, for which an overall normalisation is not needed. For any given galaxy within a cell, the overall normalisation is given by its redshift and apparent magnitudes.

In its original configuration, the VESPA age grid runs between 0.002 and 14 Gyr, in logarithmically spaced time bins. This, however, means that the width of the final bin is around 5 Gyr. This in turn is comparable to the stretch of lookback time that the LRG sample probes (roughly between z=0.5 and z=0.1) and we therefore have very little sensitivity to the exact formation epoch of LRGs. The original configuration was so designed because typical SDSS data for individual objects does not normally allow us to distinguish between populations that are separated in age by less than the bin widths.

The large SNR obtained from stacking, however, opens up the possibility to use a revised grid that has more sensitivity to older populations, at the expense of losing time resolution at the young end (the number of bins and overall structure must remain the same, for reasons constrained by the algorithm). We therefore construct an alternative age grid that we use to analyse the LRGs with VESPA. The boundaries of this new grid are: 0.002, 0.074, 0.177, 0.275, 0.425, 0.658, 1.02, 1.57, 2.44, 3.78, 5.84, 7.44, 8.24, 9.04, 10.28, 11.52 and 13.8 Gyr. Whereas these seem ambitious, statistically the data should be able to differentiate between models in adjacent bins: the differences between models of adjacent ages (normalised to a common value at a given wavelength) are several times that of the statistical noise in the stacked data. In practice, however, systematic errors can dominate - see the next Section for details.

Finally note that given the high SNR of the stacked spectra, VESPA always runs to its highest time-resolution.

VESPA was designed to deal with limitations coming from photon-noise. However, we are purposely putting ourselves in the limit where the photon-noise is negligible, and the quality of the fits is dominated by limitations in the modelling. We do not have a set of models that give statistically good fits in the limit of high-SNR and this has two potential effects: i) an error budget dominated by systematics in the modelling inaccuracies; and ii) the recovered parameters are biased.

To deal with the model inaccuracies, we make the assumption that if the models (SSPs, dust and parametrization) could represent galaxies perfectly, then the recovered solution would match the stacked spectrum exactly. Under this assumption, the residuals from a fit in each cell can be interpreted as an overall error including statistical photon-noise errors and systematic model errors.

1. 1.

   we run VESPA on the stacked spectrum and errors given by equations (2), and calculate the residuals $r_{\lambda}=|\hat{F}_{\lambda}-F_{\lambda}|$;

2. 2.

   we smooth $r_{\lambda}$ with a box-car filter to keep the large-scale variation but attenuate pixel-to-pixel fluctuations;

3. 3.

   we take the smoothed $r_{\lambda}$ as the new error, $\sigma^{total}_{\lambda}$ and re-fit the stacked spectrum using this new error.

Fig. 3 shows a typical example of the residuals and the new estimated error. Therefore we allow extra freedom in the regions where the models cannot fit the data, by assuming that the residual between data and best-fit model gives an indication of the 1-sigma confidence region for systematic problems with the models, combined with the statistical error. As covariances will vary between different models, we assume this systematic is uncorrelated for flux values at different wavelengths, to avoid placing further biases on acceptable solutions. This estimate of the combined error replaces the instrumental statistical error in all of our fits. Our estimated systematic error will almost certainly be too small, because the best-fit parameters were chosen to match the data, and minimise this difference. Fundamentally, this is a problem that cannot be solved without better models. However the methodology outlined in this section should at least provide a approximation for the error and will be better than not allowing for this problem.

We translate this error into an error bar on the recovered physical parameters by constructing $N$ random realisations of the best-fit spectrum, using the error calculated as per this section, and analysing these random realisations using VESPA. From this set of $N$ solutions we can estimate the variance and co-variance of all the parameters of interest. More details are given in Tojeiro et al. (2007, 2009).

While these residuals provide a relative comparison of our current ability to model galaxies by comparing solutions obtained with different sets of models (SSPs, dust, etc), we must be careful not to interpret any variation in the recovered parameters as an absolute error bar - different models can be wrong in the same way. We analyse these differences in Section 6.

Finally, we should also note that the difference between models in adjacent bins (see Section 3.3) can be of the order of - or smaller than - this estimate of the systematic error. Ignoring the wavelength dependence, the average systematic error would correspond to roughly a 1 Gyr error on an 8 Gyr population.

VESPA assumes that the star-formation rate (SFR) is constant within a time bin, and this gives results with discontinuities in SFR across bin boundaries. This is an obvious over-simplification of the problem, as any star formation history is likely to be continuous over an arbitrary set of boundaries, if sampled at high enough resolution. For the purposes of this paper, this over-simplification is problematic because a discontinuous SFR will result in a physically discontinuous colour evolution.

To solve this problem, once VESPA has produce model parameters, we replace the top-hat bins of constant SFR with Gaussian distributions. The widths and heights are chosen to be such that the mass formed within each pair of the old sharp boundaries does not change by more than 0.3%. The mass integrated over all ages is kept constant. These new Gaussian representations for each population (previously given by a single top-hat bin) are sampled at a much higher rate in lookback time, such that we can calculate a smooth colour evolution. Each new Gaussian bin has the same metallicity as the old corresponding histogram.

We do not expect this change in the SFR as a function of time to make a large difference in the interpretation of the spectra. This is because, by construction, the original boundaries are such that the spectral evolution between each boundary is not too significant. The offset between the spectra calculated top-hat and Gaussian bins is either smaller or of the order of the statistical error shown in Fig. 3. The flux obtained with the Gaussian bins is only used to compute the observed colours and their evolution with redshift, and the effect of this difference in the computed colours is negligible.

Our goal is to compute the rest-frame spectrum of each cell (representing an ensemble of galaxies) at any point $t_{e}>t_{g}$, where $t_{g}$ is the lookback time of the cell in question. We compute

where $f(t,Z)$ is the flux of a population of age $t$ and metallicity $Z$, and $m_{k}=g(t_{k})\Delta t_{k}$ is the mass formed at $t_{k}$, assuming a SFR given by $g(t)$. For the case where $t_{e}=t_{g}$ then we simply recover the observed spectrum, as shown in the section above.

From $F_{\lambda}(t_{e})$ we can easily compute any observed frame colours, by applying the relevant K+e corrections (fully known). As noted before, we cannot compute apparent magnitudes due to a lack of normalisation. However, for any individual galaxy the normalisation is set by the observed apparent magnitude and redshift, and $F_{\lambda}(t_{e})$ gives then not only the colour, but also the luminosity evolution of that galaxy, assuming that the SFR in the galaxies follows the weighted average within the cell.

Dust evolution does not leave a footprint in a galaxy’s spectrum today, so it does not affect our modelling. However, in order to compare the colour tracks presented in Section 5.3.2 to observed populations of galaxies at different redshifts, we do need to treat the evolution of dust. Note that for the two-parameter dust model of Charlot & Fall (2000), in addition to the inter-stellar dust component, young stars are allowed extra dust extinction to account for the effects of their birth cloud, which takes some time to dissipate, of the order $t_{BC}$. Both of these components are expected to evolve with redshift. We can envisage three possible solutions to this problem:

1. 1.

   we can assume that the values we recover for $\tau_{BC}$ and $\tau_{ISM}$ remain valid for the entire history of the galaxy. I.e., when $t_{e}$ is such that a population becomes younger than $t_{BC}$ then we apply the value of $\tau_{BC}$ obtained at $t_{g}$, and we always apply the same value of $\tau_{ISM}$ to all population; and

2. 2.

   we can estimate the evolution of $\tau_{BC}$ and $\tau_{ISM}$ from values observed for particular galaxies as a function of redshift, and apply these accordingly.

3. 3.

   we can do the comparison in a dust-extinction corrected colour-colour or colour-redshift plane, in which case no dust evolution should be applied to the colour tracks.

There are advantages and disadvantages to all methods. In the case of (i), we are over-simplifying the problem by ignoring a potential evolution with redshift, but we are keeping the information that relates to each individual cell. In the case of (ii) we are using a model for the evolution of dust, but we are assuming that the galaxies that we observe at low- and high-redshift are part of the same population. This latter assumption is exactly what we are looking to avoid by tapping into the fossil record. Option (iii) is probably the most robust, but requires modelling of galaxies at high redshifts. In this paper we plot the colour-tracks over observed colours, and for simplicity we apply the correction as given by (i) above. For other applications, colour evolution can be supplied with or without dust evolution.

It is worth summarising what the outputs of the analysis described in this section are. They can be broadly grouped into two types. On one hand we have measured quantities that relate to the position of the cell at a given redshift and colour, e.g. dust information ($\tau_{ISM},\tau_{BC}$), mass-weighted metallicities or mass-weighted ages. On another hand we have the colour-tracks as a function of redshift, which are inferred evolved quantities. All of these results have a dependence on the SSP synthesis codes. We describe the SSP models we use in the next section.

With the BC03 models (Bruzual & Charlot, 2003) we adopt a Chabrier initial mass function (Chabrier, 2003) and Padova 1994 evolutionary tracks (Alongi et al., 1993; Bressan et al., 1993; Fagotto et al., 1994a, b; Girardi et al., 1996). We use the BC03 models with the empirical STELIB library (Le Borgne et al., 2003) in the optical (3200Å  to 9500Å) and the theoretical BaSeL 3.1 stellar library (Lejeune et al., 1997, 1998; Westera et al., 2002) to either side of that range.

The M10 models (Maraston & Stromback , submitted) are the new, high-resolution version of the M05 models (Maraston, 2005). The energetics and stellar evolution calculations are unchanged from the models in Maraston (1998) and M05, so the post main-sequence continues to be treated with the fuel-consumption theorem (Renzini & Buzzoni 1986; Buzzoni 1989; Maraston 1998), leading to a more robust treatment of advanced stages of stellar evolution, such as the TP-AGB phase. The M10 models and publication focus not only on delivering a high-resolution version of the M05 models, but also on exploiting differences arising from stellar libraries and their calibration. Therefore, models are provided in a range of stellar libraries, each with different coverage in age, metallicity and wavelength. Here we choose the MILES library (Sánchez-Blázquez et al., 2006), as it provides the best coverage in age and metallicity. The wavelength coverage was extended to the UV following the method in Maraston et al. (2009). We adopt a Kroupa IMF (Kroupa, 2001).

The Flexible Stellar Population Synthesis (FSPS) code (Conroy et al., 2009; Conroy & Gunn, 2009) takes a novel approach to the computation of SSP spectra. By parametrizing uncertain stages in stellar evolution and allowing these parameters to change freely, the authors attempt to both quantify our ignorance about stellar evolution in terms of derived galaxy properties, but also to calibrate their models by finding the combination of parameters that best describe certain observations. Furthermore, models are provided with a choice of isochrones and stellar library. Here we chose the latest Padova evolutionary tracks calculations (Marigo & Girardi, 2007; Marigo et al., 2008), and the MILES stellar library. An UV extension was obtained by using the theoretical BaSeL spectral library. In spite of the freedom provided by these models, here we use them simply at their default settings and with the combination of parameters that describe the TP-AGB phase that were found to best fit star cluster data by Conroy & Gunn (2009). The opportunity remains to explore further the flexibility of these models - we leave that for a follow up paper. We use a Chabrier (2003) IMF.

The metallicities implicit in the SPS models and quoted throughout this paper refer to [Fe/H] abundances. However, it is now observationally well established that ETGs and bulges have a larger abundance of the so-called $\alpha$-elements (O, Ne, Mg, Si, S, Ca, Ti) than that predicted by SPS models, given the measured Fe abundances (e.g Worthey et al. 1992; Davies et al. 1993; Fisher et al. 1995; Longhetti et al. 2000; Maraston et al. 2003; Thomas et al. 2003; Thomas et al. 2005 and references therein).

$\alpha-$elements are mostly produced through the collision of $\alpha$-particle nuclei in type-II supernovae of very massive stars, whereas Fe-peak elements are associated with type-Ia supernovae, typically associated with older and lower mass stars (e.g. Woosley & Weaver 1995). Given the different delay times for the two types of stellar explosions, different [$\alpha$/Fe] abundances have therefore been interpreted as a signature for the length of time over which the star formation happened (e.g. Pagel & Tautvaisiene 1995; Thomas et al. 2005). The mismatch with the Fe abundances predicted by SPS models therefore occurs because empirical stellar libraries are primarily constructed by stars in the solar neighbourhood, which has a different chemical enrichment and star formation history than that which is typical of bulges and massive elliptical galaxies.

Whereas some effort has gone into calibrating certain spectral indices for different $\alpha$-element abundances (e.g Thomas et al. 2003), this has not yet been done for the full spectrum. Different chemical abundances are therefore a known limitation of current SPS models based on empirical libraries, if one wants to use a full-spectral fitting technique. Theoretical libraries however, can provide some insight on how different chemical abundances ratios affect galactic spectra (e.g. Coelho et al. 2007). For the moment, such discrepancies (surely to be systematic when looking at a population such as the LRGs) are automatically bundled with other unknown systematic errors in Section 3.4.

In Fig. 4 we show a fit for one of the stacks, and in Fig. 5 the corresponding residuals in units of the stack noise, as given by Eq. (2). This fit was obtained using the FSPS models (results for the same stack, with Mastro and BC03 models can be found in Appendix A). In this case we can see a visually good match on the blue end, and an increasing tension between data and models towards the red. The detail of the residuals plot reveals that there are significant departures from the data at all wavelengths, but more notoriously blueward of 3800Å and redwards of 5800Å .

Fig. 5 also shows regions of the spectrum that were excluded from fits that give the results presented in this section. Inclusion or exclusion of these regions does not affect the results in this section in any significant way.

In this section we show how ages, dust and metallicity vary with colour, luminosity and redshift in our LRG samples. The luminosity ranges refer to those shown in Fig. 1.

When we show or quote errors in this Section, they refer to variations in colour, for a given redshift bin and luminosity slice.These variations are typically larger than the systematic error obtained as per Section 3.4, which is likely to be an underestimation of the true systematic error. The effect of the statistical error on the recovered solutions is negligible.

A summary of the results in this section can be found in Fig. 6.

We can use the full star-formation history of each cell to estimate the colour evolution of a galaxy within that cell. This information can be coupled with the apparent magnitude and exact redshift of the galaxy to calculate its luminosity evolution. In this section we present a selection of these colour tracks for a sample of galaxies selected in colour, luminosity and redshift. As before, we present the results using only one set of SPS models in this section (FSPS), and provide the corresponding plots using the other two SPS models in Appendix A.3.

Fig. 6 summarises the in-situ measurements presented in Section 5.2.

We see a different behaviour in the metallicity evolution when using different models, but all models consistently show LRGs to be metal-rich, as expected following the fiducial model. This difference is not a surprise - metallicity is known to be more model dependent than age (e.g. Panter et al. 2008; Tojeiro et al. 2009), and even the simple fact that the models sample the metallicities at different intervals will have some effect.

The sharp increase in metallicity seen when using BC03 models is most likely due to the age-metallicity degeneracy. BC03 predicts a modest amount of recent star-formation at low-redshift, but this fraction increases to higher redshifts (5th panel). The results is a mass-weighted formation time that gets younger with redshift, and an increasing metallicity. It is worth noting that the fraction of mass in young/intermediate stars is not generally uniformly distributed over the last 3 Gyr of the galaxies. Instead, VESPA generally assigns the mass to one or two age bins, and the particular bins with the most mass change from stack to stack. In the case of BC03 and M10, these are firmly at the intermediate ages regime. This is probably indicative of sporadic star formation episodes rather than a continuous star formation type of history, for any given stack.

The FSPS and M10 models paint a scenario much closer to the fiducial model of passive stellar evolution, and a constant high-metallicity. However, the M10 models require a small amount of intermediate star formation at low redshift (a little over 10% by mass for the faintest objects, and slightly less for brightest ones), which in turn drives down the mass-weighted ages of the faintest objects at low redshift. Note that this is different from the reason that drives down the mass-weighted ages of the faintest objects in FSPS - in this case it is the dominant, oldest populations that are slightly younger. This type of distinction is possible due to having time-resolved star formation histories, and it becomes crucially important if one aims to use the age-redshift relationship of LRGs to put tight constraints on $H(z)$ (e.g. Jimenez et al. 2003; Simon et al. 2005; Crawford et al. 2010; Stern et al. 2010).

Aside from the dust results, we see little coherent dependence on luminosity, in either age, metallicity or amount of recent star formation. Work on ETGs has consistently shown that fainter objects tend to have younger stars on average (e.g. Thomas et al. 2005; Carson & Nichol 2010; Zhu et al. 2010 ), but the reason why we do not see clearly in this study is easily understood. The range in luminosity of our LRG sample is much smaller than the traditional ETG samples, which tend to go down to fainter magnitudes, leaving us with a small lever-arm in luminosity with which to constrain this type of relationship. Furthermore, as the LRG sample is neither magnitude nor volume limited, any trend with luminosity over redshift is difficult to interpret. We do see a tendency for fainter objects to be younger, in both the FSPS and M10 cases (although for different reasons, as explained in the previous paragraph), but this is limited to $z\lesssim 0.3$. This is explained both by the lack of faint galaxies in the sample at larger redshifts, and by the larger range in magnitudes in the sample at low redshifts.

Dust extinction results are relatively consistent across different models. This is likely to be due to the fact that the dust modelling is common in all three cases, and that the large-scale changes due to the presence of dust are independent enough of large-scale changes due to age or metallicity. We consistently see a very clear separation in luminosity - fainter objects have more dust extinction than bright objects - and an increase towards high-redshift. In the case of Mastro and FSPS, we also see evidence for an upturn at low redshifts, suggesting that perhaps there is some contamination at his end by dustier galaxies (with some recent star formation in both cases, although up to 10 times more, in fractional mass, in the case of the Mastro models). Additionally, dust extinction is the only of our measurements to show a clear trend with colour - redder cells tend to have larger extinction, across all luminosities.

Our particular interest in this paper is to see how the solutions discussed in the previous section translate into the colour evolution of LRGs of different colour, luminosity and redshift. Inevitably, the model dependence observed in the previous section will manifest itself in the colour tracks.

The most obvious is the effect of recent star formation. An episode of star formation results in both a short and abrupt change in the colour tracks, and a reddening of the colours to higher redshift. The magnitude of these two effects depends on the size (in mass), age and duration of the episode of star formation. Metallicity and dust induce less abrupt changes in the colour evolution.

The interpretation of the results from the different sets of models is rather different. In the case of FSPS, with negligible amounts of recent star formation and a constant metal-rich population, the colour tracks follow the locus of observed objects very closely. As the selection box was designed around a passively evolving model, this is no less the expected behaviour. Furthermore, we see little dependence on colour or luminosity, suggesting that the differences noted in the previous section (namely dust and age of oldest population) have a small effect on the colour evolution.

It is a different matter in the case of the M10 or BC03 models. With BC03, the substantial and ever-present (in redshift, colour and luminosity) level of recent star formation means that star formation events constantly push the colour tracks outside the contours. Depending on the duration of the burst (and we have limited resolution with our grid) then an object remains outside the selection for box for a given amount of time, until the excessively blue light of the young stars subsides enough to make the object red once again. In principle, there is nothing wrong with this scenario, and it follows that if LRGs go through small and sporadic events of star formation, then their true number density is greater than the number density measured inside the selection box at any given redshift. The fraction of objects not observed due to this is also in principle calculable, with an analysis of this type and good, trustworthy modelling.

In the case of the M10 models, this effect is much less important due to the smaller and more restricted amount of recent star formation. M10 models tend to be redder (especially in $r-i$), but mostly lie within the locus of observed galaxies, even when a small amounts of star formations is detected.

Changes in colour tracks can not be attributed solely to the differences in the physical solutions. These were obtained by fitting the stacked spectra in the optical ($\lambda\gtrsim 3200$Å, depending on the models), and we compute the colours to a redshift range that surpasses the optical band, so UV extensions were used to accomplish this. Differences in the spectral libraries, and in the treatment of blue stragglers or blue horizontal branch stars affect the inferred colour evolution in this regime, even when if they have limited influence in the optical.

A natural question to ask is: which set of models best describes real galaxy spectra? Answering this question, however, is far from trivial. We can use the distribution of the average residual per pixel (as calculated in Section 3.4 and shown for one stack in Fig. 3), as the simplest possible measure of goodness of fit. We show the distribution of these values for our 124 stacks and for all three SPS models in Fig. 16. There is a factor of roughly 1.5 between the mean residual of FSPS and BC03 or M10 models. This suggests FSPS models provide a better match overall for fitting the high SNR data. The physical solutions from this set of models are also the most in line with the LRG fiducial model, although the M10 solutions also broadly fit in with the mostly passive, metal-rich typical LRG picture. However the interpretation of this result is not straightforward and there are a number of caveats.

We are in a regime where neither set provides a statistically good fit to the data, and we do not know what the true answer should be. This means that when interpreting the best-fit solution and its residual when compared against the data, we need to bear in mind that the fitting procedure can compensate for deficiencies in the models by choosing a wrong solution. This behaviour can be explicitly seen in the ages of LRGs, where the best fit solution in data-space is known to be incorrect (see Section 5.2.1 for more details). In that case, we have strong physical reason to impose a prior on the age of LRGs based on the age of the Universe obtained from independent methods. An approach based on strong physical priors, or on requiring different wavelength ranges to give consistent answers, can help to identify specific regions where each model is more likely to fail. These regions can in turn be interpreted in terms of stellar physics and library deficiencies, and be a useful tool for SPS models. Having a wider range of galaxy populations will also provide a broad test of the models. We leave this analysis for another publication.