text:
# A probabilistic formulation for Empirical Population Synthesis: Sampling methods and tests

## 1 Introduction

The stellar populations of galaxies carry a record of their star-forming and chemical histories, from the epoch of formation to the present. Their study thus provide a powerful tool to explore the physics of galaxy formation and evolution. Despite the detailed understanding of stellar evolution and of the properties of single stellar populations such as star clusters achieved over the past century, uncovering the history of galaxies through the study of their integrated light has proven to be a difficult problem (Pickles 1985; Worthey 1994 and references therein).

Two global approaches have been developed to tackle this issue. The first of these is Evolutionary Population Synthesis, which performs ab initio calculations of the spectral evolution of galaxies on the basis of stellar evolution theory, stellar spectral libraries, and prescriptions for the initial mass function, star formation rate and chemical evolution (Tinsley 1972; Larson & Tinsley 1974, 1978; Guiderdoni & Rocca-Volmerange 1987; Fioc & Rocca-Volmerange 1997; Charlot & Bruzual 1991; Bruzual & Charlot 1993; Leitherer et al. 1996, 1999 and references therein). The other approach is Empirical Population Synthesis (EPS hereafter), also known as ‘Stellar Population Synthesis with a Database’, which uses observed properties of stars or star clusters as a base on which to decompose the population mix in galaxy spectra (Spinrad & Taylor 1971; O’Connell 1976; Faber 1972; Pickles 1985; Bica 1988; Pelat 1997, 1998; Boisson et al. 2000). Of particular interest is the EPS method of synthesizing the spectra of galaxies using a base of star clusters. Perhaps the most appealing property of this method is the fact that the results do not hinge on assumptions about stellar evolution tracks nor initial mass functions, as these are ‘mother nature given’ in the empirically built base spectra. This technique has been pioneered by the work of Bica (1988, hereafter B88), where the stellar content of a large sample of galaxies was synthesized using the spectral base of star clusters described in Bica & Alloin (1986a, 1986b, 1987). Since then, several studies made use of this technique (Schmidt, Bica & Dottori 1989; Bica, Alloin & Schmidt 1990; Jablonka, Alloin & Bica 1990; Schmidt, Bica & Alloin 1990; Jablonka & Alloin 1995; De Mello et al. 1995; Bonatto et al. 1996, 1998, 1999, 2000; Schmitt, Storchi-Bergmann & Cid Fernandes 1999; Kong & Cheng 1999; Raimann et al. 2000).

In its original version (B88), this method made use of a large base of $`n_{}=35`$ elements, combined in different proportions ($`x_i`$; $`i=1\mathrm{}35`$) to synthesize a number of equivalent widths ($`W`$) of conspicuous absorption features. An uniform sweep of the parameter space along paths consistent with simple predefined chemical evolution scenarios was performed, and all combinations which produced $`W`$’s within 10% of the observed values were considered ‘solutions’. The final population vector $`𝐱`$ was taken to be the arithmetic mean of all such sampled solutions.

One of the main concerns raised by B88 method is that of the uniqueness of the solution, as first discussed by Schmidt et al. (1991, hereafter S91). With a 35-D population vector and at most 9 $`W`$’s to be synthesized, one has a highly degenerate algebraic system. As shown by S91, however, the original base was highly redundant, with several linearly dependent elements, and others which were practically so, as deduced by a principal component analysis. S91 then proposed the use of a reduced base of 12 elements. S91 also criticized B88’s uniform sampling of the $`𝐱`$ space (a ‘discrete combination procedure’ in their terms) and the use of a ‘mean solution’. Instead, they developed a constrained minimization procedure, which, as shown by their simulations with test galaxies built out of the base, is able to retrieve with good accuracy the correct population vector. Another difference of the S91 method compared to B88 is the possibility to search for solutions along the entire age $`\times `$ metallicity plane. B88 assumed that the chemical enrichment of the galaxy occurred during its formation ($`10`$ Gyr ago), with all stars younger than that having the same $`Z`$ as the oldest most metal rich component. Although this may be a reasonable approximation for the nuclear region of isolated galaxies, it does not allow for the effects of mergers, inflows or outflows, and is not applicable to high redshift sources, where the spectra usually integrate over the entire galaxy, including regions with different star formation and chemical histories (Jablonka et al. 1990; Jablonka & Alloin 1995).

The next improvement in this technique was the introduction of continuum colors in the synthesis. This was done by Schmitt, Bica & Pastoriza (1996), who, as B88, use a discrete combination procedure search for acceptable solutions, sweeping the parameter space along chemical enrichment paths. This version of EPS has been used in several subsequent studies (Bonatto et al. 1996, 1998, 1999, 2000; Kong & Cheng 1999). Another step was taken by Bonatto, Bica & Alloin (1995), who extended the star cluster base to the space ultraviolet (1200–3200 Å). They too use an arithmetic average to represent the solution of the synthesis.

In two recent papers, Pelat revisited the EPS problem of synthesizing $`W`$’s using an elegant formalism based on convex algebra concepts. In Pelat (1998) he shows how the algebraic solution of the EPS problem can be narrowed down to a sub-space of $`𝐱`$ delimited by a set of ‘extreme-solutions’, any convex combination of which yields an exact solution in the underdetermined cases (more parameters than constraints). His results shed new light into the issue of algebraic degeneracy, which has always haunted EPS due to the widespread belief that with more parameters than observables one is bound to fit the data in one way or another. He proposes a quick test of whether the $`W`$’s of a galaxy are synthesizable by computing the region in $`W`$-space spanned by all physical combinations of the base elements. We verified that in the case of the 12 elements base used by S91 this region is narrow (Leão 2001), such that small measurement errors are enough to place objects outside this ‘synthetic domain’, hence preventing the existence of mathematically exact solutions. In real applications, one therefore expects the degeneracy of EPS to be more of a statistical nature than algebraic. In these cases, as well as in overdetermined problems, Pelat (1997) demonstrates that the best model is to be found along the boundary of the synthetic domain by minimizing an appropriately defined $`\chi ^2`$-like distance. Uncertainties in the population vector can also be readily computed, as shown by Moultaka & Pelat (2000).

Boisson et al. (2000) recently applied Pelat’s method to a sample of 12 active galaxies from Serote Roos et al. (1998). They synthesize $`n_W=47`$ absorption lines with a base of $`n_{}=30`$ stars from the stellar spectral libraries of Serote Roos, Boisson & Joly (1996) and Silva & Cornell (1992). Their results nicely illustrate the power of this new EPS technique, which will certainly be extended to larger samples and other classes of galaxies in future studies. The B88 method and its variants, on the other hand, already have a large body of literature associated to it. Nevertheless, this popular and attractive EPS technique, while intuitively acceptable, has neither been adequately tested nor cast onto a sound mathematical formalism which supports its application. Furthermore, it gives no measure of the uncertainties or potential biases on the population vector, thus giving no simple way to prevent overinterpretations of the resulting population of the age $`\times `$ metallicity plane. These issues are the focus of the present paper.

Our main goals are to:

* Elaborate a mathematically consistent formulation of the EPS problem in the context of probability theory.
* Develop and test a sampling method to explore the parameter space in EPS problems. In particular, we aim to improve upon the Schmitt et al. (1996) technique of synthesizing both equivalent widths and colors, but within the context of a probabilistic formulation.
* Use the method to investigate the effects of measurement errors and of the specific set of observables used upon the results of the synthesis with the reduced B88 star cluster base described in S91. In particular, we wish to evaluate whether these factors introduce biases in the population mix inferred from this popular synthesis technique and at what level of detail can its results be trusted.

This paper is organized as follows: Item (i) of the above list is addressed in section 2. In section 3 we deal with point (ii). This is done by investigating the applicability of a Metropolis algorithm to sample the parameters probability distribution in EPS, and testing the method with simulations based on Bica’s spectral base. In section 4 we address item (iii) by means of extensive simulations which explore the ability to recover the detailed star-formation and chemical histories of galaxies with this base. Finally, section 5 summarizes our main results.

## 2 Formalism

### 2.1 Basic equations

EPS studies seek to find combinations of a spectral base which reproduce a given set of measured observables, often taken as the equivalent widths $`W_j`$ of $`n_W`$ conspicuous absorption features. The base consists of $`n_{}`$ elements representing well defined simple stellar populations such as star-clusters, with equivalent widths $`W_{ij}^{}`$ and corresponding continuum fluxes over the lines $`F_{ij}^{}`$ ($`j=1,\mathrm{}n_W`$, $`i=1,\mathrm{}n_{}`$) normalized at a reference wavelength. Denoting by $`x_i`$ the fractional contribution of the $`i`$-th base element to the total flux at the reference wavelength, one obtains a system of $`n_W`$ equations

$$W_j=W_j(𝐱)=\frac{_{i=1}^n_{}W_{ij}^{}F_{ij}^{}x_i}{_{i=1}^n_{}F_{ij}^{}x_i};j=1,\mathrm{},n_W$$
(1)
for the synthetic $`W`$’s (e.g., Joly 1974). This set of constraints can be augmented by synthesizing $`n_C`$ observed continuum fluxes $`C_k`$ ($`k=1,\mathrm{}n_C`$), provided allowance is made for reddening, here parametrized by the V-band extinction $`A_V`$. Assuming all $`n_{}`$ populations are equally reddened, one obtains

$$C_k=C_k(𝐱,A_V)=g_k(A_V)\underset{i=1}{\overset{n_{}}{}}C_{ik}^{}x_i;k=1,\mathrm{},n_C$$
(2)
where a distinction is made between $`C^{}`$ and $`F^{}`$, since the $`n_C`$ fluxes to be synthesized need not correspond to the wavelengths of the $`n_W`$ absorption lines. We hereafter refer to the $`C_k`$’s as continuum colors, as they are in fact ratios of continuum fluxes with respect to the continuum at the normalization wavelength. The function $`g_k(A_V)`$ reddens the normalized color at wavelength $`\lambda _k`$ by $`A_{\lambda _k}`$ according to a specified reddening law.

Finally, physical solutions must further satisfy the normalization and positivity constraints:

$$\underset{i=1}{\overset{n_{}}{}}x_i=1;x_i0\mathrm{for}i=1,\mathrm{},n_{},\mathrm{and}A_V0.$$
(3)
The normalization condition effectively reduces one degree of freedom. When modeling colors, however, one has to introduce $`A_V`$ as a further parameter, so that the number of parameters still is $`n_{}`$.

### 2.2 Probabilistic formulation

The data $`𝒟`$ to be modeled are thus composed of a set of $`n_{obs}=n_W+n_C`$ observables. The measurement errors in each observable, collectively denoted by $`\sigma `$, are assumed to be known. Given these, the problem of EPS is to estimate the population vector $`𝐱`$ and the extinction $`A_V`$ that ‘best’ represents the data according to a well defined probabilistic model. It is equally important to estimate the uncertainties in the parameters, as these prevent over-interpreting the resulting mixture of stellar populations at a level of detail not warranted by the data or by intrinsic limitations of the base.

The probability of a solution $`(𝐱,A_V)`$ given the data $`𝒟`$ and the errors $`\sigma `$, is given by Bayes theorem:

$$P(𝐱,A_V|𝒟,\sigma ,)=\frac{P(𝒟|𝐱,A_V,)P(𝐱,A_V|\sigma ,)}{P(𝒟|)}.$$
(4)
$``$ summarizes the set of assumptions on which the inference is to be made. They include: the mapping between parameters and observables is given by eqs. (1) and (2); $`𝐱`$ and $`A_V`$ must satisfy the constraints expressed in eq. (3); the stellar population in the target galaxy is well represented by the base elements; the observational errors are Gaussian; the reddening law is known.

The likelihood $`P(𝒟|𝐱,A_V,\sigma ,)`$ is a measure of how good (or bad) is the fitting to the data for model parameters $`𝐱`$ and $`A_V`$. Under the hypothesis of Gaussian errors, the probability of the data given the parameters is

$$P(𝒟|𝐱,A_V,\sigma ,)e^{},$$
(5)
with $``$ defined as half the value of $`\chi ^2`$:

$`(𝐱,A_V)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\chi ^2(𝐱,A_V)={\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{n_W}{}}}\left({\displaystyle \frac{W_j^{obs}W_j(𝐱)}{\sigma (W_j)}}\right)^2`$
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{k=1}{\overset{n_C}{}}}\left({\displaystyle \frac{C_k^{obs}C_k(𝐱,A_V)}{\sigma (C_k)}}\right)^2`$
The likelihood, as the total $`\chi ^2`$, separates into a $`W`$ and a color related term, of which only the latter depends on $`A_V`$.

The normalizing constant $`P(𝒟|)`$ in Bayes theorem (the ‘evidence of $``$’) is irrelevant to the level of inference discussed here, i.e., the estimation of $`𝐱`$ and $`A_V`$. It can be used, for instance, to compare different spectral bases.

$`P(𝐱,A_V|)`$ is the joint a priori probability distribution of $`𝐱`$ and $`A_V`$, and states what values the model parameters might plausibly take. For instance, physically acceptable population fractions should have priors that are zero in regions of the parameter space where the constraints of positivity and unit sum are not satisfied. A prior that is uniform on the parameters and includes these constraints is a non-informative prior, because it does not impose any constraints on the solutions besides those expressed by eq. (3). This is the prior that will be used in this work.

For a non-informative prior the posterior probability $`P(𝐱,A_V|𝒟,\sigma ,)`$ is simply proportional to the likelihood:

$$P(𝐱,A_V|𝒟,\sigma ,)e^{(𝐱,A_V)}$$
(7)
This expression contains the full solution of the EPS problem, as embedded in it is not only the most probable model parameters but also their full probability distributions. Furthermore, projected posterior distributions for any of the $`n_{}`$ parameters can be obtained by marginalizing eq. (7) with respect to all other parameters. For instance, the probability density of proportion $`x_i`$ of the $`i`$-th base element is

$`P(x_i|𝒟,\sigma ,)`$ $`=`$ $`{\displaystyle \mathrm{}P(𝐱,A_V|𝒟,\sigma ,)𝑑x_1\mathrm{}}`$
$`dx_{i1}dx_{i+1}\mathrm{}dx_n_{}dA_V`$
and equivalently for $`A_V`$. Similarly, one can construct joint posteriors for, say, the total proportion of all base components of same age or $`Z`$, or for the mean age and $`Z`$ of the stellar population. This would of course provide a coarser description of a galaxy’s stellar content than the individual posteriors $`P(x_i|𝒟,\sigma ,)`$, but that may be all that is possible under some circumstances, such as when only a reduced set of observables is available or when observational errors are large (section 4).

It is easy to see how the B88 and Schmitt et al. (1996) synthesis techniques fit into this general probabilistic formulation. By synthesizing the observables within $`10\%`$ ‘error boxes’, these authors implicitly assumed a box-car likelihood function $`P(𝒟|𝐱,A_V,\sigma ,)`$, whereas by performing arithmetic means over all accepted $`(𝐱,A_V)`$ combinations in a uniform grid they are implicitly sampling the corresponding posterior probability distributions. Also, the a priori constraints on the occupation of the age $`\times `$ $`Z`$ plane imposed by these authors (but relaxed in subsequent works) simply reflect their use of an informative prior. This gives some formal justification for their heuristically designed EPS method. It is therefore reasonable to expect that an implementation of EPS based on the formalism presented here should give results roughly compatible with those previously obtained. An uniform sweep of the ($`𝐱,A_V`$) space is however a very inefficient sampler for (2.2), so some improvement is needed there. This is discussed next.

## 3 Sampling Method and Tests

Despite its formal merits, at the computational level our probabilistic approach faces the same basic difficulty as previous EPS codes, namely, the high dimensionality of the problem. For spectral bases with astrophysically interesting resolution in age and metallicity the number of elements $`n_{}`$ quickly becomes large enough to render the exploration of the parameter space a non-trivial task. Before discussing methods to sample the ($`𝐱,A_V`$) space, we present the spectral base used in this work (§3.1) and describe how we deal with the uncertainties in the observables (§3.2).

### 3.1 The spectral base and observables

The base used for the synthesis here is that of S91. It contains $`n_{}=12`$ population groups, spanning five age bins—10 Gyr (which actually represent globular cluster-like populations), 1 Gyr, 100 Myr, 10 Myr and HII (corresponding to current star formation, and represented by a pure $`F_\lambda \lambda ^2`$ continuum based on the spectrum of 30 Dor)—and four metallicities—0.01, 0.1, 1 and 4 Z (Table 1). The observables in this base comprise the $`n_W=9`$ equivalent widths of the absorption lines CaII K $`\lambda `$3933, CN $`\lambda `$4200, G band $`\lambda `$4301, MgI $`\lambda `$5175, CaII $`\lambda `$8543, CaII $`\lambda `$8662, H$`\delta `$, H$`\gamma `$ and H$`\beta `$, as well as $`n_C=7`$ continuum fluxes at selected pivot wavelengths: 3290, 3660, 4020, 4510, 6630, 7520 and 8700 Å, all normalized to 5870 Å. According to Bica & Alloin (1986a, 1987) and Bica, Alloin & Schmitt (1994), the equivalent width windows and continuum points were defined based on very high signal to noise spectra of galaxies, with the express goal of using them to synthesize the stellar population of galaxies. These spectra were obtained by creating two average spectra, of blue and red galaxies, where the stellar population features can be clearly traced. The base values for these quantities have been previously published in Bica & Alloin (1986b, 1987) and Schmitt et al. (1996), and some were measured from data in Bica et al. (1994). These are recompiled in Table 2 , as they bear a direct impact on the analysis below. There are thus a maximum of $`9+7=16`$ observables to model with $`n_{}=12`$ parameters: $`n_{}1`$ population fractions plus $`A_V`$. However, in several of the tests presented below we shall make use of a reduced subset of observables, as observational data sets seldom cover the whole spectral range spanned by this base.

According to Bica & Alloin (1986a, 1987) and Bica, Alloin & Schmitt (1994), since the ultimate goal of the continuum points and $`W`$’s is to synthesize the stellar population of galaxies, the $`W`$ windows and continuum points were defined based on very high signal to noise spectra of galaxies. These spectra were obtained by creating two average spectra, of blue and red galaxies, where the stellar population featuress can be identified.

The reddening of the colors is modeled with the extinction law described in Cardelli, Clayton & Mathis (1989, with $`R_V=3.1`$).

### 3.2 Errors in the observables

The measurement errors in the observables play a key role in determining the structure of the likelihood function. Whereas such errors are available when analyzing observed galaxy spectra, they have to be postulated in the synthesis of test galaxies performed below. An alternative and potentially more attractive approach to deal with the errors is to marginalize eq. (2.2) over $`\sigma `$. This yields more conservative estimates for the parameters, insofar as they are independent of the actual observational errors. In this paper we follow the traditional approach of treating $`\sigma `$ as a relevant constraint in the synthesis process.

In order to insure a realistic correspondence between the values of the errors and the actual quality of the data we have adopted a prescription for the $`\sigma `$’s based on our experience in dealing with the measurement of $`W`$’s and $`C`$’s in galaxy spectra (Cid Fernandes, Storchi-Bergmann & Schmitt 1998). Two sources of uncertainty affect the measurements of $`W`$’s: (1) the noise within the line window, and (2) the uncertainty in the positioning of the pseudo-continuum fluxes $`C_k`$ at the pivot wavelengths, which affect $`W_j`$ because they define the continuum over the line. The effect of the first of these sources upon $`\sigma (W_j)`$ is straight-forwardly obtained through standard propagation of errors by specifying the noise level at $`\lambda _j`$, the spectral resolution $`\delta \lambda `$, and the size $`\mathrm{\Delta }_j`$ of the pre-defined window over which $`W_j`$ is measured. In Bica & Alloin (1986a) system the pseudo continuum is defined interactively by visual inspection of the spectrum around the pivot $`\lambda `$’s, which makes the estimation of $`\sigma (C_k)`$ non straight-forward. In Cid Fernandes et al. (1998) we found these errors could be roughly scaled from the signal-to-noise ratio measured in nearby ‘line-free’ regions, such that the $`S/N`$ in the pseudo continuum is $`(S/N)_{\mathrm{PC}}\rho _{\mathrm{PC}}(S/N)_\lambda `$, with $`\rho _{\mathrm{PC}}2`$–3. This leads to the following expression for the $`\sigma (W_j)`$’s:

$$\sigma ^2(W_j)=\mathrm{\Delta }_j\delta _\lambda (S/N)_\lambda ^2+(\mathrm{\Delta }_jW_j)^2\rho _{\mathrm{PC}}^2(S/N)_\lambda ^2$$
(9)
where the two terms correspond to the two sources of uncertainty discussed above. For simplicity, we further assume that the noise in the observed spectrum is constant for all $`\lambda `$’s, which results in $`(S/N)_\lambda =C_\lambda (S/N)_{5870}`$. The errors in the $`C_k`$’s thus become also constant:

$$\sigma (C_k)=\rho _{\mathrm{PC}}^1(S/N)_{5870}^1$$
(10)
Besides yielding realistic values for the $`\sigma `$’s, this recipe has the advantage of quantifying the quality of the data in terms of a single quantity: The signal-to-noise ratio at the normalization wavelength 5870 Å. Of secondary importance are the spectral resolution and $`\rho _{\mathrm{PC}}`$, which in the simulations below are kept fixed at 5 Å and 3 respectively.

In order to provide a quantitative notion of how a value of $`S/N`$ translates into $`\sigma (W_j)`$’s, we have computed these quantities for the observables corresponding to spectral groups S1 and S7 of B88, typical of the nuclei of early (red) and late (blue) type spiral galaxies respectively. For S1 our error recipe yields $`\sigma (W_{\mathrm{CaIIK}})=36\mathrm{\AA }/(S/N)_{5870}`$ and $`\sigma (W_{\mathrm{MgI}})=19\mathrm{\AA }/(S/N)_{5870}`$, so errors of $`1`$ Å or less are achieved with $`(S/N)_{5870}30`$. For S7 these lines have $`\sigma (W)15\mathrm{\AA }/(S/N)_{5870}`$.

### 3.3 Exploration of the parameter space: Sampling methods

The exploration of the ($`𝐱,A_V)`$ space can be approached from two alternative perspectives: (1) A minimization problem, employing methods to search for the global $`\chi ^2`$ minimum (maximum likelihood). Minimization techniques were explored in S91 and Pelat (1997) for the synthesis of $`W`$’s only. (2) A sampling problem, where one seeks to map out the posterior probability distribution (eq. 7).

The methods discussed in this paper are sampling methods, which do not explicitly search for a minimum. In the limit of small errors, however, one expects the posterior maps to peak at the most likely model. For large errors, on the other hand, the best model sampled may be well off the true global minimum, but in this case the probability distribution is so broad that the very meaning of ‘best model’ is questionable. In such cases it is arguably more important to estimate plausible ranges for the parameters than to carry out a refined search for the most likely values.

Our main motivation to explore the sampling approach is that most applications of EPS with Bica’s base to date focused on the estimation of mean parameters, whose meaning cannot be directly ascertained without estimating their uncertainty and/or biases. A critical re-evaluation of this method is thus crucial to establish to which level of detail the results of this popular technique can be trusted.

#### 3.3.1 Uniform sampling

In order to compute the individual posterior probabilities for each parameter, we first note that eq. (2.2) for the $`x_i`$ posterior is simply the mean probability $`P(𝐱,A_V|𝒟,\sigma ,)`$ over the space spanned by all possible values of $`x_{ji}`$ and $`A_V`$ for a fixed $`x_i`$. The simplest method to compute such probabilities is to divide the ($`𝐱,A_V`$) space into a uniform grid with $`\mathrm{\Delta }x`$ steps for the population fractions and $`\mathrm{\Delta }A_V`$ for the reddening parameter. One then approximates eq. (2.2) by the finite sum

$`P(x_i|𝒟,\sigma ,)`$ $``$ $`{\displaystyle \frac{_se^{(𝐱_s,A_{V,s})}\delta (x_ix_{i,s})}{_se^{(𝐱_s,A_{V,s})}}}`$ (11)
where $`s`$ denotes a point (a ‘state’) in the 13-D grid and the $`\delta (x_ix_{i,s})`$ term retains only those grid points where population $`i`$ contributes with fraction $`x_i`$ of the light.

This discrete sweeping of the parameter space is easy to implement, and was in fact the recipe followed by most works with Bica’s base to date (see section 1). A serious drawback of uniform sampling is its computational price. The number of grid cells is an extremely steep function of the resolution $`\mathrm{\Delta }x`$. A coarse $`\mathrm{\Delta }x=10\%`$ sampling yields $`3.5\times 10^5`$ $`𝐱`$-points, increasing to $`8\times 10^7`$ for $`\mathrm{\Delta }x=5\%`$ and more than $`4\times 10^{11}`$ for a $`2\%`$ resolution! To obtain the grid size one has to further multiply these numbers by the number of points in the $`A_V`$ grid, whose limits have to be pre-defined. Schmitt et al. (1999), for instance, performed a synthesis study with $`\mathrm{\Delta }x=5\%`$ and $`\mathrm{\Delta }A_V=0.06`$ for $`A_V`$ between 0 and 1.5, which yields $`2\times 10^9`$ grid points, perhaps the finest grid ever used in EPS with this 12 elements base. One obviously always aim for the best resolution possible, but a priori estimates of what resolution is necessary are not straight-forward. A ‘good’ resolution should produce probability distributions that are smooth over scales of $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }A_V`$. The width of these distributions is of course a function of the observational errors and of the set of observables being synthesized.

#### 3.3.2 Metropolis sampling

A further drawback of uniform sampling is that the algorithm spends most of the time in regions of negligible probability, which contribute little to the integral in equation (2.2), or its numerical version (eq. 11). A more efficient sweep of the parameter space would be to use a sampling scheme which traces the full probability distribution (eq. 7). One such ‘importance sampling’ scheme is the Metropolis algorithm (Metropolis et al. 1953), which samples preferentially regions where the posterior is large. This is the parameter-space exploration method adopted in this work.

Our implementation of the Metropolis sampler was as follows. Starting from an initial arbitrary point, at each iteration $`s`$ we pick one of the 13 variables at random and change it by an uniform deviate ranging from $`ϵ`$ to $`+ϵ`$, producing a new state $`s+1`$. If the variable is one of the 12 $`x_i`$’s, the whole set is renormalized to unit sum. Moves towards unphysical values ($`x_i<0`$ or $`>1`$, $`A_V<0`$) were truncated. Downhill moves (i.e., towards smaller $`\chi ^2`$) are always accepted, whilst changes to less likely states are accepted with probability $`e^{(_{s+1}_s)}`$, thus avoiding trapping onto local minima. This scheme, widely employed in statistical mechanics (Press et al. 1992; MacKay 2001 and references therein), produces a distribution which tends to the correct one as the number of samples $`N_s`$ increase. There is, however, no universal prescription to choose optimal values for $`ϵ`$ or $`N_s`$, so some experimentation is needed.

### 3.4 Performance tests

In order to test our Metropolis EPS code we performed a series of simulations for artificial galaxies generated out of the base. Two sets of test galaxies were used. The first is composed of the 26 galaxies used by S91 (see their Table 5) to test their minimization algorithm. As they did not synthesize colors, we have used a fixed value of $`A_V=2/3`$ to redden all their galaxies. The second set is composed of 100 test-galaxies whose population vectors and extinctions were generated randomly using a scheme which insures that in many cases a few components dominate the light. All simulations presented in this paper were performed sampling $`N_s=10^8`$ states, with the ‘visitation parameter’ $`ϵ`$ set to 0.005 for both the $`x_i`$’s and $`A_V`$. Experiments were also performed changing these quantities, and found to give nearly identical results. We note that less samples can be used if one is only interested in the mean ($`𝐱,A_V`$), as this naturally converges much faster than its full probability distribution. Simulations were performed for severals values of $`(S/N)_{5870}`$ and different sets of observables. Here we restrain the discussion to the illustration of the sampling method.

In Fig. 1 we plot the individual projected posterior probability distributions for each of the 12 $`x_i`$’s plus $`A_V`$ for two of S91’s test galaxies composed of widely different stellar populations. The posteriors were computed synthesizing all 9 lines and 7 colors in the base (‘Set A’, as defined below) and using $`(S/N)_{5870}=1000`$ to define the errors in the observables. This is clearly an unrealistic, highly idealized situation, but it serves as a test-bed for the method in the limit of perfect data. The simulations were started from the $`x_i=1/12`$ and $`A_V=0.5`$ point. The plot shows that the posterior probabilities are sharp functions peaking at the expected input value, marked by upper arrows, demonstrating that the Metropolis sampler converges to the most likely region. One notes, however, that the strong components in the old age bin (such as $`x_1`$, $`x_2`$ and $`x_3`$ in the ‘old galaxy’ example, drawn as solid lines) have broader individual posteriors. This effect is further discussed in the next section.

In the opposite limit of large errors, $`\sigma \mathrm{}`$ (the ‘infinite temperature’ regime, as it would be called in statistical mechanics) the likelihood loses its discriminating power and becomes $``$ constant over the whole parameter-space. Eq. (2.2) then tends to $`P(x_i)(1x_i)^{10}`$ for a $`n_{}=12`$ base. In this regime the data does not constrain the model at all, and this expression for $`P(x_i)`$ simply reflects the volume of the $`𝐱`$-space spanned by a given value of $`x_i`$. Not surprisingly, this distribution peaks at $`x_i=0`$, and has a mean value $`\overline{x_i}=1/12`$. We have verified that the Metropolis scheme recovers this distribution satisfactorily as $`S/N0`$, which serves as a further test of the adequacy of this sampling method. The skewing of several of the individual posteriors plotted in Figs. 3 and 8 towards small $`x_i`$ is partly due to this effect, which drags the components to an intrinsically larger region of $`𝐱`$-space.

The evolution of the Metropolis walk through the 13-D ($`𝐱,A_V`$) space is illustrated in Fig. 2. For clarity, only 6 components of $`𝐱`$ and just one state for each $`10^4`$ steps are shown. Note that projected $`P(x_i)`$ posteriors such as those in Fig. 1 are essentially histograms of the $`x_i`$ states sampled in walks as that in Fig. 2, but weighted by their $`\mathrm{exp}\chi (𝐱)^2/2`$ likelihood. The input values of the components shown are indicated at the right edge of the plot for comparison. This particular example corresponds to a simulation with $`S/N=60`$ and with all 16 observables synthesized. One sees relatively large excursions of the parameters from their true values. Furthermore, anti-correlations between some of the components ($`x_1`$ and $`x_2`$, $`x_8`$ and $`x_9`$) are clearly visible. This ‘mirror effect’ becomes much more pronounced as $`S/N`$ increases, whereas for low $`S/N`$ the broadening of the likelihood function washes out such correlations. This is a consequence of redundancies between several of the base components, as further discussed below. The top panels in Fig. 2 show the evolution of the $`\chi ^2`$ and how two of the synthesized observables, $`W`$(CaIIK) and $`C_{4020}`$, oscillate very much within their respective $`\pm 1`$ sigma error ranges. One therefore anticipates an excellent fit to the data, and indeed the $`\chi ^2`$ obtained for the mean ($`𝐱,A_V)`$ ‘solution’ in this example is just 0.8 (see also Fig. 7).

Overall, we conclude that the Metropolis algorithm is an efficient sampler of the parameter space for EPS problems. It provides a much more continuous mapping of $`𝐱`$ and $`A_V`$ than would be possible with a uniform grid. Furthermore, it is substantially more efficient computationally. A uniform grid with $`10^8`$ points would only sample each $`x_i`$ at $`7\%`$ steps and 10 values for $`A_V`$. By concentrating on the relevant regions of the parameter space, we achieve a much better resolution, as can be judged by the smoothness of the posterior distributions on scales of 1–2%.

One possible variation of this scheme is to combine the uniform grid and Metropolis routines, starting a Markov chain from each point on a uniform grid. Experiments with this alternative scheme were performed, and found to produce essentially the same results. This is to be expected, since we verified that the starting point has little effect upon the resulting posterior distributions. In fact, we performed a series of simulations using the input parameters as the starting point and obtained identical results. Though further improvements are possible, the simple Metropolis sampler already provides a substantial improvement over previous works with this base.

## 4 Tests of the base

Having developed a probabilistic formulation of EPS and a new sampling method, we now present a series of numerical experiments designed to evaluate the actual ability to recover the stellar population mix of galaxies using Bica’s base. The tests were performed with the same series of fictional galaxies described in section 3.4. Simulations were performed for $`(S/N)_{5870}=10`$, 30, 60, 100, 300 and 1000, and three different sets of observables:

* Set A: All 16 observables in the base.
* Set B: Only the 9 absorption line equivalent widths.
* Set C: The $`W`$’s of CaII K, CN, G-band and MgI+MgH plus the continuum fluxes at 3660, 4020, 4510 and 6630, all normalized to 5870 Å.

Set A is the ideal, as it uses all information in the base. Set B is used as a test of how much is actually gained by synthesizing colors as well as $`W`$’s, while Set C is composed of the observables used by Cid Fernandes et al. (1998) and Schmitt et al. (1999) to characterize the stellar content of active galaxies and their hosts through long-slit spectroscopy. These data will be used to perform a spatially resolved EPS study in a future communication.

### 4.1 Effects of the errors in the observables

We first investigate the effects of the measurement uncertainties in the observables upon the results of the synthesis. Qualitatively, one expects the degradation of the data quality to broaden the probability distributions. Furthermore, a shift and skewing of the $`P(x_i)`$’s of intrinsically large components towards small $`x_i`$’s is also expected for low $`S/N`$ due to intrinsically larger number of states in this region of $`𝐱`$ (the ‘dragging effect’ discussed in §3.4). This is exactly what is observed in Fig. 3, where we overplot the individual posteriors for $`S/N=1000`$, 100 and 10 for one particular test galaxy of S91, synthesized with Set A observables. As in Fig. 1, the true parameters are indicated by the top arrows. The broadening of posteriors in this case can be safely attributed to the errors alone, since, as demonstrated by Pelat (1998) most of the test galaxies in S91 (including the one in Fig. 3) have a unique algebraic solution in a $`W`$-only EPS problem, despite the fact that the number of degrees of freedom (11 if we do not model the colors) falls short of the number of observables (9 $`W`$’s).

Despite the large number of constraints in Set A, one sees that the individual posterior distributions are substantially broad even for a $`S/N=100`$ spectrum. The average standard deviation of the $`x_i`$’s for this combination of $`S/N`$ and observables is 3% over all simulations, but can reach more than 10% in individual components. Furthermore, the dominant components, 2, 6 and 9 in this example, are all badly affected by the errors. Fig. 3 shows that decreasing data quality progressively shifts $`P(x_2)`$ towards smaller values, whereas the opposite happens with $`x_1`$. This compensation is also seen between components in the 1 Gyr ($`x_5`$, $`x_6`$ and $`x_7`$), 100 Myr ($`x_8`$ and $`x_9`$), and 10 Myr ($`x_{10}`$ and $`x_{11}`$) age bins and is the same ‘mirror’ effect observed before in Fig. 2. A compelling visualization of these compensations is given in Fig. 4, where we plot the individual Metropolis states for the test galaxy employed in Fig. 3 in a $`x_i\times x_j`$ matrix. The $`S/N`$ of 300 in this plot was chosen for clarity purposes; the same structure is present for other values, with the scatter around the strongest correlations increasing steadily as $`S/N`$ decreases, to the point that at $`S/N=10`$ the anti-correlations between components 1 and 2, or 8 and 9 become clouds similar to $`x_1\times x_5`$ in Fig. 4. All anticorrelations identified in Fig. 3 ($`x_1\times x_2`$, $`x_6\times x_5`$, etc.) are clearly represented in this graphical version of the correlation matrix.

These effects are rooted in the internal structure of the base. By definition, the base elements must be linearly independent, and the S91 base complies with this algebraic condition. Yet, some of its elements are practically linearly dependent, in the sense that combinations of other elements can recover them to a high degree of accuracy. As the noise increases, the residuals between representing an element by its exact $`W`$’s and $`C`$’s and a combination including other elements with more extreme values for the observables become statistically insignificant, explaining the compensations detected above. To demonstrate this we synthesized the base elements themselves with our code. Results for Set A observables and three different values of $`S/N`$ are shown in Fig. 5. Each panel corresponds to one of these 12 extreme test galaxies (labeled B1…B12), with the mean synthetic values of each of the components plotted along the vertical axis. The empty, shaded and filled histograms correspond to $`S/N=100`$, 30 and 10 respectively. Some of the components, most notably numbers 2, 3 and particularly 6, are synthesized with large contributions of others even for $`S/N=100`$. For the $`x_6=100\%`$ model (top left panel of Fig. 5), for instance, we find $`\overline{x_6}=60\%`$ for this high $`S/N`$, with 35 of the remaining 40% redistributed among components 5 and 7. Yet, the 16 observables are very well reproduced by this combination, with a total $`\chi ^2`$ of just 3.8. This spreading of the light fractions becomes more pronounced for lower $`S/N`$’s, to the point that at $`S/N=10`$, components like 8 and 9, or 10 and 11 cannot be distinguished anymore, and end up dividing their contributions in roughly equal shares, with residuals spread over other components. The spread is much smaller in the synthesis of components 1, 4, 5 and 7, which represent extreme metallicities within age groups. For intermediate $`Z`$ systems (models B2, B3 and B6), however, these extreme components attract much of the percentage light contribution spilled over from neighboring populations of same age, so that their individual contributions in a real stellar population mix are as unreliable as the others.

This experiment demonstrates that for practical purposes the base is still linearly dependent, at least in a statistical sense, despite the effort of S91 in reducing the highly redundant original base of B88. This ‘noise-induced’, or ‘statistical linear dependence’, as we may call it, could indeed be inferred from the PCA analysis of S91, which showed that very little of the variance is associated with the last few eigenvectors. Within our probabilistic formulation, this ‘statistical linear dependence’ defines the structure of the likelihood function, which spreads first along directions in $`𝐱`$-space producing similar observables, carrying the mean $`𝐱`$ to more densely populated regions. This effect explains the redistribution of the probability observed in Figs. 3, 4 and all other simulations.

It is therefore extremely difficult to retrieve accurately all the input stellar population parameters even for excellent data. This is further illustrated in Fig. 6, where the input parameters are plotted against their mean values, as derived from the Metropolis runs. One hundred artificial galaxies are plotted in each panel, with open circles corresponding to Set A observables and $`S/N=1000`$, and crosses indicating the results for Set A and $`S/N=60`$. A systematic underestimation of strong components is seen for $`S/N=60`$, a good spectrum by any standard. This happens at the expense of an overestimation of the weaker components, which produces the large scatter seen in the bottom left corners of Fig. 6. The average uncertainties in the $`x_i`$’s are of order $`\sigma (x_i)3`$–4% for $`S/N=60`$. These are not enough to account for the large differences between input and output $`x_i`$’s depicted in Fig. 6, so the errors induce a true bias in $`\overline{𝐱}`$. Note also that, in accordance with the discussion above, this bias is smaller for models with strong contributions of components 1, 4, 5 and 12, because of their extreme locations in the space of observables.

Only for unrealistically large $`S/N`$’s, which far exceed the quality of the data used to build the base in the first place, one is able to break the ‘statistical linear dependence’ of the base by distinguishing fine details in the observables. This should not be interpreted as a failure of the method, as what the code actually synthesizes are the observables! These are very precisely reproduced by the mean ($`𝐱,A_V`$), as illustrated in Fig. 7 for the $`W`$’s of CaII K, CN, G-band and MgI and two continuum colors. The remaining observables, not show in the figure, are equally well reproduced. Furthermore, whilst $`S/N>300`$ is needed to recover accurately the model parameters from the mean solution, the agreement between the synthetic and measured observables is excellent for any $`S/N`$. Obviously, an even better agreement is obtained if instead of $`(\overline{𝐱},\overline{A_V})`$ the best model sampled during the Metropolis excursion is used to reconstruct the observables, as illustrated in the fourth column of plots in Fig. 7 (see §4.5).

In conclusion, this comparison of input and output observables shows that the spreading of the probability and the confusion between components seen in Figs. 3 and 5 is not an artifact of the method, but a consequence of the internal structure of the spectral base. As found in other EPS studies (O’Connell 1996 and references therein) synthesizing the observables is one thing, but trusting the detailed star-formation history and chemical evolution implied by the synthesis is an altogether different story. On the positive side, a notable fact about Fig. 5 is that, consulting Table 1, one realizes that the reshuffling of the strength of the components occurs preferentially among populations of same age, an effect also clearly seen in Fig. 4. Though some redistribution among base elements of different age and $`Z`$ in the older 1 and 10 Gyr bins due to the age-$`Z`$ degeneracy also occurs, this indicates that the age distribution may be well recovered by the synthesis (§4.3).

### 4.2 Effects of different sets of observables

The effects of which set of observables is used in the synthesis are illustrated in Fig. 8, where results for Sets A, B and C for the same test galaxy as in Fig. 3 and $`S/N=1000`$ are compared. The figure shows that it is very advantageous to synthesize colors along with equivalent widths, despite the fact that one increases the dimensionality of the problem with the inclusion of $`A_V`$. This can be seen by comparing the performances of Sets A and B. The information contained in the colors not only improves the estimation of $`𝐱`$ but also allows a very good determination of $`A_V`$, which, unlike the population vector, is always well recovered by the synthesis.

Fig. 8 also shows how important it is to provide as many observational constraints as possible to constrain the synthesis, as Set A recovers the parameters much better than Sets B and C. Decreasing the number of observables in the synthesis thus have the same overall effect as decreasing the $`S/N`$. The width of the posteriors for Set C are partially attributed to the algebraic degeneracy in this set (8 observables and 12 degrees of freedom), which implies the existence of exact solutions in a sub-space of ($`𝐱`$,$`A_V`$). That is not the case of Set A for the test galaxy in this example (Pelat 1998), but there, as for other sets, the ‘statistical linear dependence’ of the base is very efficient in spreading the likelihood even for such an idealized S/N ratio. The excellent agreement between the ‘observed’ and synthesized observables even for much worse S/N (Fig. 7) illustrates this point. Algebraic degeneracy is not as critical as the ‘statistical degeneracy’ induced by the combined effects of noise, limited data and the internal correlations in the base.

Studies which synthesize only $`W`$’s sometimes use the resulting population vector to compute a predicted continuum shape and derive $`A_V`$ a posteriori through the comparison with the observed continuum (e.g., B88; Boisson et al. 2000). From our Set B simulations, we find that this procedure yields values of $`A_V`$ 2–3 times less accurate than obtained with Set A. For $`S/N=60`$, for instance, Set A yields an rms dispersion in $`A_V`$ of 0.06 mag around the true value, whereas the dispersion for set B is 0.13 mag. Since this is a relatively small difference, these experiments validate the a posteriori computation of the extinction. Colors are more useful to constrain the population vector than to estimate $`A_V`$.

### 4.3 Age grouped results

The results of the tests above reveal a striking difficulty to accurately determine all 12 base components of Bica’s base in the presence of even very modest noise and/or when not all base observables are available for the synthesis. As indicated by the simulations, both the measurement errors and the use of reduced sets of observables act primarily in the sense of spreading a strong contribution in one component preferentially among base elements of same age. Grouping the population vector in age bins should thus produce more robust results. This expectation is confirmed in Fig. 9, where, analogous to what was done in Fig. 6, we plot the input values of the $`x_i`$’s against the output $`\overline{x_i}`$’s, but now for the five age binned groups, also for Set A simulations. The left panels correspond to the same values of $`S/N`$ used in Fig. 6, so that one can appreciate the enormous improvement achieved by grouping equal age elements. The right panels show that good results are also obtained for $`S/N`$ of 30 (circles), and that young populations ($`100`$ Myr) are reasonably well traced even for $`S/N`$ as low as 10 (crosses). Though some deviations survive, specially for the smaller $`S/N`$’s and the older age bins, the age-binned synthesis results are much more reliable than the results for the individual components, and, most importantly, do not require unrealistic $`S/N`$.

This conclusion also holds for Sets B and C, as illustrated in Fig. 10. Naturally, the restricted input data in these sets translates into a loss of information, and the agreement is not as good as for Set A. Note that Set C does a much better job in recovering populations of 100 Myr or less than Set B, a further example of how useful colors are in the synthesis. This happens because the color information (present in Set C but absent in B) impose stronger constraints on young populations than in older systems.

The reliability of age-grouped proportions contrasts with the badly constrained and biased results for the individual components of $`𝐱`$, and is one of the main results of our study. This result puts in perspective all previous interpretations of the synthesis with this base. This fundamental limitation of the base was in fact previously known, but was never studied in detail. Schmitt et al. (1999), for instance, preferred to summarize the results of their synthesis study in terms of age grouped proportions for populations younger than 1 Gyr, in consonance with the results above. Still, in that work we kept a distinction between the different $`Z`$’s in the 10 Gyr age-group. At this level of detail the population fractions are very uncertain. In fact, since Set C observables and $`S/N30`$–60 were used in that study, Fig. 10 indicates that it would be safer to group their 1 and 10 Gyr proportions. Yet, the proportions for the younger ($`100`$ Myr) populations found by Schmitt et al. (1999), and which constituted the focus of that paper, are trustworthy (Fig. 10).

#### 4.3.1 The age $`\times `$ metallicity degeneracy

A close inspection of the $`S/N=10`$ and 30 results in Fig. 9 reveals a compensation effect between the 1 and 10 Gyr age bins, albeit at a much smaller level than for the individual components (Fig. 6). This effect is much more pronounced in Sets B and C, as seen in Fig. 10.

The confusion between the 1 and 10 Gyr populations is related to the age-$`Z`$ degeneracy, which sets in precisely in this age range. This well known effect (O’Connell 1986, 1994; Worthey 1994) is also present in Bica’s base and affects the synthesis results. To visualize and quantify this effect we have computed the average log $`Z`$ and log age for both the input and synthetic $`𝐱`$ for our test galaxies. Given that $`𝐱`$ is a light fraction at 5870 Å, these averages ultimately represent flux-weighted ages and metallicities, which serve our current illustrative purposes, despite their ambiguous physical significance.

The ‘output - input’ $`\mathrm{log}\mathrm{age}`$ and $`\mathrm{log}Z`$ residuals are plotted against each other in Fig. 11. In these diagrams $`\mathrm{\Delta }\mathrm{log}\mathrm{age}>0`$ (overestimated age) implies redder colors and stronger metal lines, which tend to be counter-balanced by the bluer colors and weaker lines resulting from $`\mathrm{\Delta }\mathrm{log}Z<0`$ (underestimated $`Z`$). The age-$`Z`$ degeneracy is most obvious in the results for Sets B and C, for which the biases in age and $`Z`$ persist even for $`S/N=100`$. For Set A, on the other hand, the residuals are small: $`\genfrac{}{}{0pt}{}{_<}{^{}}`$ 0.2 dex for $`S/N`$ $`\genfrac{}{}{0pt}{}{_>}{^{}}`$ 60, which shows that the age-$`Z`$ degeneracy can be broken, at coarse the level of age and $`Z`$ resolution offered by the base, with enough information and good spectra. Note that this remark applies to the global averages defined above, not the detailed component by component ages and $`Z`$’s, which we concluded to be highly uncertain (Fig. 6).

### 4.4 Metallicity grouped proportions

The results for $`Z`$ binned proportions, plotted in Fig. 12, are not nearly as good as those for the age binned groups. With the exception of the $`Z=10^2Z_{}`$ bin, which is actually represented by just one element (Table 1), the scatter in the input $`\times `$ output $`Z`$-binned proportions is so large for all other $`Z`$’s that one is forced to conclude that no accurate description of the chemical history of galaxies can be afforded with this spectral base. Given this limited $`Z`$ ‘resolution power’, one might consider removing intermediate components such as 2, 3 and 6. However, as we concluded that only the age distribution can be assessed with the base, there is no obvious advantage in this further reduction. This situation can be improved by providing a priori constraints on the occupancy of the age-$`Z`$ plane, like those imposed by B88 based on chemical evolution arguments, since these effectively reduce the dimensionality of the base to typically 8 components, thus alleviating the confusion between components and producing better focused results. The validity of this prior input external to the synthesis process, nonetheless, has to be evaluated in an object by object basis. In this study we followed S91 in not imposing any such extra information, as this approach encompasses a larger class of evolutionary scenarios and thus contemplates more possible applications.

### 4.5 Best $`\times `$ Mean Parameters

Throughout these experiments, we have consistently used the mean parameters as an estimate of the result of the synthesis. This convention was deliberately chosen to follow the more widely employed version of EPS with this base and thus to allow a reassessment of previous results. Furthermore, the mean is a convenient and reasonable way to summarize the results insofar as it represents a ‘center of mass’ of acceptable solutions. The severe biases identified in the mean synthetic population vector, however, raise the question of whether alternative approaches should be pursued instead.

By definition, minimization oriented techniques would do a better job at recovering the input parameters. Indeed, we verified that the best model population vector sampled in the Markov chain is much closer to $`𝐱_{true}`$ than its expected value, an agreement which can only be improved given that the Metropolis sampler is not an optimal minimization procedure. Though the simulations confirmed our expectation that for small errors the mean and best solutions should both converge to the true parameters, it was surprising to realize how easily $`\overline{𝐱}`$ moves away from $`𝐱_{true}`$, and consequently also from $`𝐱_{best}`$. Very small uncertainties in the observables, at the levels of $`S/N>300`$, are enough to set powerful compensation effects into operation (e.g., Figs. 3 and 4), shifting the components of $`\overline{𝐱}`$ away from their true values due to the asymmetric redistribution of the probability in $`𝐱`$-space.

Minimization procedures can in principle overcome this bias, as illustrated by the success of the S91 and Pelat (1997) tests as well as by our (not optimized) results for $`𝐱_{best}`$. A further bias in mean solutions is that they always produce $`\overline{x_i}>0`$. Real galaxies, even if synthesizable by the base, are likely to have observables outside synthetic domain due to measurement errors. As shown by Pelat (1997, 1998), the best model in this case lies at the surface of the synthetic domain, where very often several of the $`𝐱`$-components are 0. In these aspects, minimization algorithms are more suitable than the more traditional sampling approach. However, the high susceptibility of the synthesis to the data quality indicates that minor perturbations in the input observables, within their errors, may seriously affect the estimation of the best model parameters. The stability of the best solution is thus likely to be more fragile than that of the mean solution, which takes into account the effects of the errors in the observables in the trade-off between an optimal and a statistically acceptable fit of the data. We have in fact detected this effect in a series of simulations, but a more careful study, with specialized minimization techniques is required to quantify it properly. It would be particularly interesting to carry out this test with the method developed by Pelat (1997) applied to the base used in this work. His Monte Carlo tests with a $`n_{}=10`$ stellar base and $`n_W=19`$ equivalent widths produced essentially unbiased population fractions even for $`S/N`$ as low as 10. Extending his formalism to synthesize galaxy colors would also be desirable.

In any case, because of the nature of the EPS problem and the structure of Bica’s base, relatively large discrepancies in the parameters do not necessarily translate into large differences in the synthesized quantities. In this respect, the difference between mean and best models is largely academic, as the mean solution already does an excellent job in fitting the observables (Fig. 7).

## 5 Summary

In this paper we have revisited the method of Empirical Population Synthesis with the aims of improving it both at the formal and computational levels, and exploring its efficacy as a tool to probe the population mixture of galaxies. Our results can be divided into three parts.

(1) A simple probabilistic formulation of the problem was presented, which puts this method, traditionally employed in a less formal way, onto a mathematical footing. It was shown how former applications of the method fit into the probabilistic formulation, thus providing some a posteriori justification for previous results.

(2) An importance sampling scheme, based on the Metropolis algorithm, was developed and tested. It provides a more efficient and smooth mapping of the probability distribution of the parameters than is possible with the commonly used uniform grid sampling.

(3) In the third part of the paper, we applied the formulation and sampling method to a series of test galaxies constructed out the $`n_{}=12`$ star-clusters base of S91. The tests explored the ability to reconstruct the model parameters, which carry a record of the star formation and chemical evolution in galaxies, in the presence of (i) observational errors and (ii) limited data. This study centered on the comparison of input parameters with the mean synthetic solution ($`\overline{𝐱},\overline{A_V}`$), as this is the most common form of EPS in the literature with this base, and a systematic evaluation of the consistency of this method has not been carried out before. The main results of this study can be summarized as follows:

* (a) The allowance for errors in the observables sparkle linear dependences within the base elements. This induces systematic biases in $`\overline{𝐱}`$, redistributing the probability in components with a large fractional contribution to the integrated spectrum among similar components. This happens preferentially for components of same age. Though this was a predictable result, it was surprising to realize how powerful this effect is. For any realistic $`S/N`$, the bias is such that in most cases one cannot trust the estimates of all 12 individual components of $`\overline{𝐱}`$ to any useful level of confidence.
* (b) Reducing the number of observational constraints, due to incomplete coverage of the 3300–9000 Å spectral interval spanned by the base, has the same overall effect as reducing the data quality, that is, to broaden, shift and skew the probability distribution of the $`x_i`$’s. We find that synthesizing colors as well as equivalent widths yields substantially better results than synthesizing $`W`$’s only. The need to account for the extinction as an extra parameter is largely compensated by the better focused $`𝐱`$ obtained by considering the color information, particularly for populations of 100 Myr and younger. Unlike for the individual population fractions, we do not detect systematic biases in the estimates of $`A_V`$.
* (c) Despite the difficulty to retrieve accurately the detailed population vector, the observables are very precisely reproduced by the mean solution. The mapping between the parameters and observables spaces is thus highly ‘degenerate’. The nature of this ‘non-uniqueness’ is not algebraic, as it happens also when the number of free parameters exceeds the numbers of observables and in undetermined cases where the solution is known to be unique. Instead, it is related to the possibility of mimicking to a statistically indistinguishable level the effects of certain base elements by linear combinations of other elements. This happens even for tiny observational errors, inducing a ‘spill over’ of the likelihood among base components, carrying the mean solution along.
* (d) The structure of the base reflects the fact that evolution is the main source of variance in stellar populations. Accordingly, we found that grouping the $`\overline{𝐱}`$ components in their five age bins produces much more reliable results, and we recommend this procedure when evaluating the detailed population vector obtained in previously published results with this base. Some confusion between the 1 and 10 Gyr populations also takes place due to the age-$`Z`$ degeneracy, but this can be minimized with an adequate spectral coverage (measuring all observables in the base) and $`S/N>30`$. Grouping the results for the four different $`Z`$’s in the base is also recommended, but caution must be exerted when interpreting the results since these are much less precise than the age-binned proportions.

The difficulties to retrieve accurately the true stellar population parameters stem from two sources: (i) The convention to represent the results of the synthesis by the mean solution, and (ii) limitations imposed noise-induced linear dependences in the base. Minimization techniques may be instrumental in overcoming the biases in the mean population vector, which ultimately limit the resolution of the 12-D base to its 5 age bins. Within our probabilistic formulation, other priors can be investigated. Maximum Entropy, for example, is a prior that has been applied successfully to inverse problems formally similar to population synthesis, like image restoration. However, more than improvements on the method, which are possible and desirable, further work on the base itself will certainly be needed. Indeed, the limited age and $`Z`$ ‘resolution power’ of the base employed here lies at the heart of all difficulties identified in this paper. Including other spectral indices, modeling the full spectrum and incorporating other relevant information such as mass to light ratios in the synthesis process are possibles ways to progress in this direction.

## ACKNOWLEDGEMENTS

We thank Catherine Boisson, Monique Joly, Jihane Moultaka and Didier Pelat for stimulating discussions on the virtues and drawbacks of EPS techniques. RCF thanks the hospitality of Johns Hopkins University, where this work was finalized. Support for this work was provided by the National Science Foundation through grant # GF-1001-99 from the Association of Universities for Research in Astronomy, Inc., under NSF cooperative agreement AST-9613615. JRSL acknowledges an MSc fellowship awarded by CAPES. Partial support from CNPq, FAPESP, PRONEX and FUNPESQUISA-UFSC is also acknowledged.