text:
# Toward reconstruction of dynamics of the Universe from distant supernovae type Ia

## I Introduction

Recent observations of type Ia supernovae Riess et al. (1998, 2002); Perlmutter et al. (1999) supported by WMAP measurements of anisotropy of the angular temperature fluctuations Benoit et al. (2003) indicate that our Universe is spatially flat and accelerating. On the other hand, the power spectrum of galaxy clustering Parcival et al. (2002) indicates that about $`30\%`$ of critical density of the universe should be in the form of non-relativistic matter (cold dark matter and baryons). The remaining, almost two thirds of the critical energy, may be in the form of a component having negative pressure (dark energy). Although the nature of dark energy is unknown, the positive cosmological constant term seems to be a serious candidate for the description of dark energy. In this case the cosmological constant $`\mathrm{\Lambda }`$ and energy density $`\epsilon _\mathrm{\Lambda }=\mathrm{\Lambda }/(8\pi G)`$ remain constant with time and the corresponding mass density $`\rho _\mathrm{\Lambda }\epsilon _\mathrm{\Lambda }/c^2=6.44\times 10^{30}(\mathrm{\Omega }_\mathrm{\Lambda }/0.7)(h/0.7)g`$ $`cm^3`$, where $`h`$ is the Hubble constant $`H_0`$ expressed in units of $`100`$ $`\mathrm{km}`$ $`\mathrm{s}^1\mathrm{Mpc}^1`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7\pm 0.1`$. Although the cold dark matter (CDM) model with the cosmological constant and dust provides an excellent explanation of the SNIa data, the present value of $`\mathrm{\Lambda }`$ is $`10^{123}`$ times smaller than value predicted by the particle physics model. Many alternative condidates for dark energy have been advanced and some of them are in good agreement with the current observational constraints Ratra and Peebles (1998); Peebles and Ratra (2003); Caldwell et al. (1998); Kamenshchik et al. (2001); Armendariz-Picon et al. (2000); Bucher and Spergel (1999); Parker and Raval (1999). Moreover, it is a natural suggestion that $`\mathrm{\Lambda }`$-term has a dynamical nature like in the inflationary scenario. Therefore, it is reasonable to consider the next simplest form of dark energy alternative to the cosmological constant $`(w=1)`$ for which the equation of state depends upon time in such a way that $`p=w(a(z))\rho `$, where $`wp/\rho `$ is a coefficient of the equation of state parametrized by the scale factor or redshift. It has been demonstrated Szydlowski and Czaja (2003a, b) that dynamics of such a system can be represented by one-dimensional Hamiltonian flow

$$=\frac{\dot{a}^2}{2}+V(a),$$
(1)
where the overdot means differentiation with respect to the cosmological time $`t`$ and $`V(a)`$ is a potential function of the scale factor $`a`$ given by

$$V(a)=\frac{\rho _{\mathrm{eff}}(a)}{6}a^2,$$
(2)
where $`\rho _{\mathrm{eff}}`$ is the effective energy density which satisfies the conservation condition

$$\dot{\rho }_{\mathrm{eff}}=3\frac{\dot{a}}{a}(\rho _{\mathrm{eff}}+p_{\mathrm{eff}}).$$

For example, for the $`\mathrm{\Lambda }\mathrm{CDM}`$ model we have

$$\rho _{\mathrm{eff}}=\mathrm{\Lambda }+\rho _{m,0}a^3$$

$$p_{\mathrm{eff}}=\mathrm{\Lambda }+0.$$
(3)
Of course the trajectories of the system lie on the zero energy surface $`0`$.

Hamiltonian (1) can be rewritten in the following form convenient for our current reconstruction of the equation of state from the potential function $`V`$, namelly

$$(p,x)=\frac{p_x^2}{2}+V(x),$$
(4)
where $`p_x=L/\dot{x}`$, here overdot means differentiation with respect to some new reparametrized time $`t\tau :d\tau =|H_0|dt`$, $`xa/a_0=1+z`$.

For example, for mixture of non-interacting fluids potential $`V(x)`$ takes the form

$$V(x)=\frac{1}{2}\underset{i}{}\mathrm{\Omega }_{i,0}x^{23(w_i+1)},$$
(5)
where $`p_i=w_i\rho _i`$ for $`i`$-th fluid and $`i`$ $`w_i=\mathrm{const}.`$ (similar to the quiessence model of dark energy).

Due to the Hamiltonian structure of Friedmann-Robertson-Walker dynamics, with the general form of the equation of state $`p=w(a(z))\rho `$, the dynamics is uniquely determined by the potential function $`V(a)`$ (or $`V(x)`$) of the system. Only for simplicity of presentation we assume that the universe is spatially flat (in the opposite case trajectories of the system should be considered on the energy level $`=\frac{1}{2}\mathrm{\Omega }_{k,0}`$).

Let us note that from the potential function we can obtain the equation of state coefficient $`w=p/\rho `$,

$$w=\frac{1}{3}\left(1+\frac{d(\mathrm{ln}V)}{d(\mathrm{ln}a)}\right).$$
(6)
The term $`d(\mathrm{ln}V)/d(\mathrm{ln}a)`$ has a simple interpretation as an elasticity coefficient of the potential function with respect to the scale factor.

Thus from the potential function $`V`$ both $`\rho (a)`$ and $`p`$ can be unambiguously calculated

$$\rho (a)=6\frac{V(a)}{a^2},$$

$$p=2\frac{V(a)}{a^2}\left(1+\frac{d(\mathrm{ln}V)}{d(\mathrm{ln}a)}\right).$$
(7)
## II The expansion scenario from supernovae distances

As it is well known in a flat FRW cosmology the luminosity distance $`D_L`$ and the coordinate distance $`r`$ to an object at redshift $`z`$ are simply related as

$$a_0r=a_0_t^{t_0}\frac{dt^{}}{a(t^{})}=\frac{D_L(z)}{1+z},$$
(8)
($`c=8\pi G=1`$ here and elsewhere). From equation (8) the Hubble parameter is given by

$$H(z)=\dot{[\mathrm{ln}a]}=\left[\frac{d}{dz}\left(\frac{D_L(z)}{1+z}\right)\right].$$
(9)
It is crucial that formula (9) is purely kinematic and depends neither upon a microscopic model of matter, including the $`\mathrm{\Lambda }`$-term, nor on a dynamical theory of gravity. Due to existance of such a relation it would be possible to calculate the potential function which is:

$$V(a)=\frac{1}{2}H^2a^2=\frac{\left[\frac{d}{dz}\left(\frac{D_L(z)}{1+z}\right)\right]^2}{2(1+z)^2}.$$
(10)
This in turn allows us to reconstruct the potential $`V(z)`$ from SNIa data. Let us note that $`V(z)`$ depends on the first derivative with respect to $`z`$ whereas $`w(z)`$ is associated with the second derivative.

Let us also note that from of the potential function for a one-dimensional particle-universe moving in the configurational $`a`$ (or $`x`$)-space can be reconstructed from recent measurements of angular size of high-z compact radio sources compiled by Gurvits et al. Gurvits et al. (1999). The corresponding formula is

$$V(a)=\frac{\left[\frac{d}{dz}\left(D_A(z)(1+z)\right)\right]^2}{2(1+z)^2},$$
(11)
where the luminosity distance $`D_L`$ and the angular distance $`D_A`$ are related by the simple formula

$$D_L(z)=(1+z)^2D_A(z).$$
(12)
Since the potential function is related to the luminosity function by relation (10) one can determine both the class of trajectories in the phase plane $`(a,\dot{a})`$ and the Hamiltonian form as well as reconstruct the quintessence parameter $`w(z)`$ provided that the luminosity function $`D_L(z)`$ is known from observations.

Now we can reconstruct the form of the potential function (10) using a natural Ansatz introduced by Sahni et al. Sahni and Starobinsky (2000). In this approach dark energy density which coincides with $`\rho _{\mathrm{eff}}`$ is given as a truncated Taylor series with respect to $`u=(1+z)`$

$$\rho _{\mathrm{DE}}=\rho _{\mathrm{eff}}=A_0+A_1u+A_2u^2.$$

This leads to

$$H(z)=H_0\sqrt{\mathrm{\Omega }_{m,0}(1+z)^3+A_0+A_1(1+z)+A_2(1+z)^2}$$
(13)
and

$$\frac{D_L(z)}{(1+z)}=\frac{1}{H_0}_0^z\frac{dz}{\sqrt{\mathrm{\Omega }_{m,0}(1+z)^3+A_0+A_1(1+z)+A_2(1+z)^2}}.$$
(14)
The values of three parameters $`A_0`$,$`A_1`$,$`A_2`$ can be obtained by applying a standard fitting procedure to SNIa observational data based on the maximum likelihood method.

The potential function (10) written in terms of $`(1+z)`$ is

$$V(a(z))=\frac{\rho _{\mathrm{eff}}a^2}{6}=\frac{1}{2}H_0^2\left[\mathrm{\Omega }_{m,0}(1+z)+A_0(1+z)^2+A_1(1+z)^1+A_2\right]$$
(15)
or in dimensionless form

$$V(x(z))=\frac{1}{2}\left[\mathrm{\Omega }_{m,0}(1+z)+\mathrm{\Omega }_{\mathrm{\Lambda },0}(1+z)^2+\mathrm{\Omega }_{t,0}(1+z)^1+\mathrm{\Omega }_{k,0}\right],$$
(16)
where $`A_i/(3H_0^2)=\mathrm{\Omega }_{i,0}`$, $`(i=\mathrm{\Lambda },t,k)`$.

Our approach to the reconstruction of dynamics of the model is different from the standard approach in which $`w(z)`$ is determined directly from the luminosity distance formula. It should be stressed out that the latter approach has an inevitable limitation because the luminosity distance dependence on $`w(z)`$ is obtained through a multiple-integral relation that loses detailed information on $`w(z)`$ Maor et al. (2001). In our approach the reconstruction is simpler and more information on $`w(z)`$ survives (only a single integral is required). Our approach is also different from the concept of reconstruction of potential of scalar fields considered in the context of quintessence Caldwell et al. (1998).

The key steps of our method are the following:

1) we reconstruct the potential function $`V(a)`$ for the Hamiltonian dynamics of the quintessential universe from the luminosity distances of supernovas type Ia;
2) we draw the best fit curves and confidence levels regions obtained from the statistical analysis of SNIa data;
3) we set the theoretically predicted forms of the potential functions on the confidence levels diagram;
4) those theoretical potential which escape from the $`2\sigma `$ confidence level is treated as being unfitted to observations;
5) we choose this potential function which lie near the best fit curve.

Our reconstruction is an effective statistical technique which can be used to compare a large number of theoretical models with observations. Instead of estimating some revelant parameters for each model separately, we choose a model-independent fitting function and perform a maximum likelihood parameter estimation for it. The obtained confidence levels can be used to discriminate between the considered models. In this paper this technique is used to find the fitting function for the luminosity distance.

The additional argument which is important when considering the potential $`V(a)`$ is that it allows to find some modification in the Friedmann equations along the “Cardassian expansion scenario” Freese and Lewis (2002). This proposition is very intriguing because of additional terms, which automatically cause the acceleration of the universe Alcaniz and Lima (1999); Avelino et al. (2003); Zhu and Fujimoto (2002, 2003). These modifications come from the fundamental physics and these terms can be tested using astronomical observation of distant type Ia supernovae. For this aim the recent measurements of angular size of high-redshift compact radio sources can also be used Gurvits et al. (1999).

The important question is the reliable data available. We expect that supernovae data would improve greatly over next few years. The ongoing SNAP mission should gives us about 2000 type Ia supernovae cases each year. This satellite mission and the next planned ones will increase the accuracy of data compared to data from the 90s. In our analysis we use the availale data starting from the three Perlmutter samples (Sample A is the complete sample of 60 supernovae, but in the analysis it is also used sample B and C in which 4 and 6 outliers were excluded, respectively). The fit for the sample C is more robust and this sample was accepted as the base of our consideration. For technical details of the metod the reader is referred to our previous two papers Szydlowski and Czaja (2003a, b).

In Fig. 1 we show the reconstructed potential function obtained using the fitting values of $`A_i`$ as well as $``$. The red line represents the potential function for the best fit values of parameters (see TABLES 2, 3 and 4). In each case the coloured areas cover the confidence levels $`68.2\%`$ ($`1\sigma `$) and $`95.4\%`$ ($`2\sigma `$) for the potential function. The different forms of the potential function which are obtained from the theory are presented in the confidence levels. Here we consider one case, namely the Cardassian model. In this case the standard FRW equation is modified by the presence of an additional $`\rho ^n`$ term, where $`\rho `$ is the energy density of matter and radiation. For simplicity we assume that density parameter for radiation is zero (see TABLE. 1). The Cardassian scenario is proposed as an alternative to the cosmological constant in explaining the acceleration of the Universe. In this scenario the the Universe automaticaly accelerates without any dark energy component. The additional term in the Friedmann equation arises from exotic physics of the early universe (i.e., in the brane cosmology with Randall-Sundrum version $`n=2`$).

The solid lines in Fig. 1 represent the potential functions for the different Cardassian scenarios of an accelerated expansion of the universe.

For visualization the quality of performed fitting procedure the Hubble diagrams are presented on Fig. 2.

The general conclusion from our statistical analysis is that the potential function has maximum at some $`z=z_0`$ (or $`a=a_0`$) and that

$$\frac{^2V}{a^2}(z)<0.$$
(17)
In the next section it will be demonstrated that the reconstructed potential function contain all informations which we need for reconstruction of the dynamics of any cosmological model on the phase plane.

## III Application – towards the reconstruction of phase portrait and model parameters

The dynamics of the considered cosmological models is governed by the dynamical system

$`\dot{x}`$ $`{\displaystyle \frac{dx}{dt}}=y,`$
$`\dot{y}`$ $`{\displaystyle \frac{dy}{dt}}={\displaystyle \frac{V}{x}},`$ (18)
with the first integral for (18) $`=0\dot{x}^2/2+V(x)=0`$. The main aim of the dynamical system theory is the investigation of the space of all solutions (18) for all possible initial conditions, i.e. phase space $`𝕄`$. In the context of quintessential models with the equation of state $`p=w(a(z))\rho `$ there exists a systematic method of reducing Einstein’s field equations to the form of the dynamical system (18) Szydlowski and Czaja (2003a). One of the features of such a representation of dynamics is the possibility of resolving of some cosmological problems like the horizon and flatness problems in terms of the potential function $`V(x)`$.

The phase space $`𝕄`$ (or state space) is a natural visualization of the dynamics of any model. Every point $`P=(x,y)𝕄`$ corresponds to a possible state of the system. The r.h.s of the system (18) define a vector field $`𝐅(P)=[y,\frac{V}{x}]`$ belonging to the tangent space $`\mathrm{T}_\mathrm{p}𝕄`$. Integral curves of this vector field define one-parameter group of diffeomorphisms $`\varphi (P)`$ called the phase flow. In the phase space the phase curves (orbits of the group $`\varphi (P)`$) represent the evolution of the system whereas the critical points $`y=0`$, $`\frac{V}{x}=0`$ are singular solutions – equilibria from the physical point of view. The phase curves together with critical points constitute the phase portrait of the system.

Now we can define the equivalence relation between two phase portraits (or two vector fields) by the topological equivalence, namely two phase portraits are equivalent if there exists an orientation preserving homeomorphism transforming integral curves of both systems into each other. Following the Hartman-Grobman theorem, near hyperbolic critical points ($``$ $`i`$ $`\mathrm{Re}\lambda _i0`$, where $`\lambda _i`$ is the appropriate eigenvalue of linearization matrix $`𝒜`$ of the dynamical system) is equivalent to its linear part

$`\dot{x}`$ $`=y,`$
$`\dot{y}`$ $`=\left({\displaystyle \frac{^2V}{x^2}}\right)_{(x=x_0,0)}(xx_0).`$ (19)
In our case the linearization matrix takes the form

$$𝒜=\left[\begin{array}{cc}0,& 1\\ \frac{^2V}{x^2},& 0\end{array}\right]_{(x_0,0)}$$
(20)
Classification of critical points is given in terms of eigenvalues of the linearization matrix since the eigenvalues can be determined from the characteristic equation $`\lambda ^2(\mathrm{Tr}𝒜)\lambda +\mathrm{det}𝒜=0`$. In our case $`\mathrm{Tr}𝒜=0`$ and eigenvalues are either real if $`\frac{^2V}{x^2}|_{(x_0,0)}<0`$ or purely imaginary and mutually conjugated if $`\frac{^2V}{x^2}|_{(x_0,0)}>0`$. In the former case the critical points are saddles and in the latter case they are centres.

The advantage of representing dynamics in terms of Hamiltonian (1) is the possibility to discuss the stability of critical points which is based only on the convexity of the potential function. In our case the only possible critical points in a finite donain of phase space are centres or saddles.

The dynamical system is said to be structurally stable if all other dynamical systems (close to it in a metric sense) are equivalent to it. Two-dimensional dynamical systems on compact manifolds form an open and dense subsets in the space of all dynamical systems on the plane Peixoto (1962). Structurally stable critical points on the plane are saddles, nodes and limit cycles whereas centres are structurally unstable. There is a widespread opinion among scientists that each physically realistic models of the universe should possess some kind of structural stability – because the existence of many drastically different mathematical models, all in agreement with observations, would be fatal for the empirical method of science Thom (1977).

Basing on the reconstructed potential function one can conclude that:
1) since the diagram of the potential function $`V(a(z))`$ is convex up and has a maximum which corresponds to a single critical point – the quantity $`\frac{^2V}{a^2}|_{(a_0,0)}<0`$ at the critical point (saddle point) and the eigenvalues of the linearization matrix at this point are real with oposite signs;
2) the model is structurally stable, i.e., small perturbation of it do not change the structure of the trajectories in the phase plane;
3) since $`\ddot{a}=\frac{V}{a}`$ one can easily conclude from the geometry of the potential function that in the accelerating region ($`\ddot{a}>0`$) $`V(a)`$ is a decreasing function of its argument;
4) the reconstructed phase portrait for the system is equivalent to the portrait of the model with matter and the cosmological constant.

By an inspection of the phase portraits we can distinguish four characteristic regions in the phase space. The boundaries of region I are formed by a separatrix coming out from the saddle point and going to the singularity and another separatrix coming out of the singularity and approaching the saddle. This region is covered by trajectories of closed and recolapsing models with initial and final singularity.

The trajectories moving in region IV are also confined by a separatrix and they correspond to the closed universes contracting from the unstable de Sitter node towards the stable de Sitter node.

The trajectories situated in region III correspond to the models expanding towards stable de Sitter node at infinity. Similarly, the trajectories in the symmetric region II represent the universes contracting from the unstable node towards the singularity.

The main idea of the qualitative theory of differential equations is the following: instead of finding and analyzing an indiwidual solution of the model one investigates a space of all possible solutions. Any property (for example acceleration) is believed to be realistic if it can be attributed to a large subsets of models within the space of all solutions or if it possesses certain stability property which is shared also by all slightly perturbed models.

## IV Conclusion

The possible existence of the unknown form of matter called dark energy has usually been invoked as the simplest way to explain the recent observational data of SNIa. However, the effects arising from the new fundamental physics can also mimic the gravitational effects of dark energy through a modification of the Friedmann equation.

We exploited the advantages of our method to discriminate among different dark energy models. With the independently determined density parameter of the Universe ($`\mathrm{\Omega }_{m,0}=0.3`$) we found that the current observational results require the cosmological constant $`n0`$ in the Cardassian models. On Fig. 1 we can see that both in the case of sample A (Fig. 1c) and sample C (Fig. 1a) $`n`$ should be close to zero. Similarly if we assume that the density parameter for barionic matter is $`\mathrm{\Omega }_{m,0}=0.05`$ then $`n0.36`$ in the case of sample C (Fig. 1b) and $`n0.39`$ for sample A (Fig. 1d). Moreover, we showed (for the sample C of Perlmutter SNIa data) that a simple Cardassian model as a candidate for dark energy is ruled out by our analysis if $`0.36n0.50`$ for $`\mathrm{\Omega }_{m,0}=0.3`$ and if $`0.48n0.13`$ for $`\mathrm{\Omega }_{m,0}=0.05`$ at the confidence level $`2\sigma `$.

Therefore, the standard $`\mathrm{\Lambda }\mathrm{CDM}`$ seems to be more appropriate for a large range of model parameters than the dark energy models based on the first version of Cardassian models.

Our main result is that the structure of the phase space of accelerating models can be effectively reconstructed from SNIa data.

###### Acknowledgements.

M. Szydłowski acknowledges the support of Jagiellonian University Rector’s Scholarship. Authors are very gratefull to dr A. Krawiec and dr W. Godłowski for discussion and comments.


format: text

text:
# A Novel Evolutionary Formulation of the Maximum Independent Set Problem

## 1 Introduction

Let $`G=(N,E)`$ be an undirected graph of node set $`N`$ and edge set $`E`$ such that $`n=|N|`$ and $`m=|E|`$. An independent set (or stable set) of $`G`$ is a subset $`S`$ of $`N`$ containing no neighbors (nodes that are connected by an edge in $`G`$). The set $`NS`$ is a node cover of $`G`$, that is, a set of nodes that includes at least one of the two end nodes of every edge. We call every subset of $`N`$ whose nodes are all neighbors of one another in $`G`$ a clique. The complement $`\overline{G}`$ of $`G`$ is an undirected graph of node set $`N`$ in which two nodes are neighbors if and only if they are not neighbors in $`G`$. Clearly, $`S`$ is an independent set of $`G`$ if and only if it is a clique of $`\overline{G}`$.

We are concerned in this paper with the problem of finding a maximum independent set in $`G`$, that is, an independent set of maximum size. Equivalently, this problem can be viewed as asking for a minimum node cover in $`G`$ or a maximum clique in $`\overline{G}`$. Finding independent sets (or any of the other equivalent structures) of extremal size has several important applications. We refer the reader to the volume that resulted from the Second DIMACS Implementation Challenge  for various examples of application areas, and to  for further examples from coding theory.

The problem of finding a maximum independent set in $`G`$ is NP-hard (it is NP-complete when formulated as a decision problem ), and remains NP-hard even if we settle for solving it approximately within $`n^{1/4ϵ}`$ of the optimum for any $`ϵ>0`$ . That is, if $`\alpha (G)`$ is the size of a maximum independent set of $`G`$ (the independence number of $`G`$), then finding an independent set of size at least $`n^{1/4ϵ}\alpha (G)`$ is NP-hard.

Our formulation of the maximum independent set problem is based on the notion of an acyclic orientation of $`G`$, i.e., an assignment of directions to the edges of $`G`$ that leads to no directed cycles. Let $`\mathrm{\Omega }(G)`$ denote the set of all the acyclic orientations of $`G`$. For $`\omega \mathrm{\Omega }(G)`$, let $`D_\omega `$ be the set of all chain decompositions of the nodes of $`G`$ according to $`\omega `$, that is, each member of $`D_\omega `$ is a partition of $`N`$ into sets that correspond to chains (directed paths) according to $`\omega `$. For $`dD_\omega `$, let $`|d|`$ denote the number of chains in $`d`$. Our point of departure is the following equality, which relates the independence number of $`G`$ to its acyclic orientations :

$$\alpha (G)=\underset{\omega \mathrm{\Omega }(G)}{\mathrm{max}}\underset{dD_\omega }{\mathrm{min}}|d|.$$
(1)
By (1), $`\alpha (G)`$ is the number of chains in the chain decomposition of $`N`$ into the fewest possible chains, according to the acyclic orientation of $`G`$ for which that number is greatest. This result is a refinement of Dilworth’s theorem  and is illustrated in Figure 1, where two acyclic orientations of the same graph are shown alongside the corresponding minimum chain decompositions. In the figure, the bottommost acyclic orientation is the one whose minimum chain decomposition is greatest, thence $`\alpha (G)=2`$ for the graph in question.

As we see it, the greatest significance of (1) is that it spells out how the set $`\mathrm{\Omega }(G)`$ can be regarded as a set of individuals, the fittest of which yields the independence number of $`G`$.111A relation dual to the one in (1) indicates how to express the chromatic number of $`G`$ (cf. ) in terms of its acyclic orientations . The use of that other relation in an evolutionary approach to find the graph’s chromatic number has been developed by one of us and collaborators . Viewing the maximum independent set problem from this perspective is based on taking $`\mathrm{min}_{dD_\omega }|d|`$ as the measure of fitness for individual $`\omega `$ and on searching $`\mathrm{\Omega }(G)`$ for an individual of maximum fitness.222$`\mathrm{\Omega }(G)`$ is a set of vast dimensions. For example, the number of distinct acyclic orientations of $`G`$ is $`2^m`$ if $`G`$ is a tree, $`n!`$ if it is a complete graph (all nodes connected to all others), and $`2^n2`$ if it is a ring. In general, and remarkably, the number of members of $`\mathrm{\Omega }(G)`$ is given by the absolute value of the chromatic polynomial of $`G`$ (cf. ) applied to the negative unit .

We develop this evolutionary approach in the remainder of the paper, starting in Section 2, where the details of the new formulation are introduced, including how to compute an individual’s fitness and the evolutionary operators of crossover and mutation for $`\mathrm{\Omega }(G)`$. We then continue in Section 3 with a description of the graphs to be used in our experiments, whose results are reported in Section 4. Conclusions are given in Section 5.

## 2 The formulation

The key to a better understanding of how to use (1) is that $`\mathrm{min}_{dD_\omega }|d|`$, the purported measure of fitness for individual $`\omega \mathrm{\Omega }(G)`$, is in fact the size of an independent set that can be unequivocally obtained from orientation $`\omega `$. In order to see this, we associate a directed graph, call it $`D(G,\omega )`$, with $`G`$ oriented by $`\omega `$. This directed graph has $`2n+2`$ nodes: two distinguished nodes, called $`s`$ and $`t`$, and two nodes, called $`i^{}`$ and $`i^{\prime \prime }`$, for each node $`iN`$. For every $`iN`$, $`D(G,\omega )`$ has an edge directed from $`s`$ to $`i^{}`$ and another from $`i^{\prime \prime }`$ to $`t`$. For every $`(i,j)E`$ that is directed by $`\omega `$ from $`i`$ to $`j`$, in $`D(G,\omega )`$ an edge exists from $`i^{}`$ to $`j^{\prime \prime }`$. Node $`s`$ is therefore a source (a node with no edges directed inward) and node $`t`$ a sink (a node with no edges directed outward). An illustration of this construction is given in Figure 2, where the directed graphs corresponding to the two acyclic orientations of Figure 1 are shown.

We now regard $`D(G,\omega )`$ as a flow network whose edges either have unit capacity (those leaving $`s`$ or arriving at $`t`$) or infinite capacity (all others), and consider the maximum flow from $`s`$ to $`t`$, whose value we assume to be $`F`$. This flow is necessarily integral, as illustrated in Figure 2, where solid lines have been used to draw edges carrying unit flow and dashed lines those carrying zero flow. It also establishes a chain decomposition of $`G`$ according to $`\omega `$ into $`nF`$ chains. Take, for example, the topmost network of Figure 2, where $`F=4`$. If we only follow edges carrying unit flow, then we can easily trace the single chain $`3^{}1^{\prime \prime }`$, $`1^{}2^{\prime \prime }`$, $`2^{}4^{\prime \prime }`$, $`4^{}5^{\prime \prime }`$. For the bottommost network, $`F=3`$ and we get the two chains $`3^{}1^{\prime \prime }`$, $`1^{}2^{\prime \prime }`$ and $`4^{}5^{\prime \prime }`$. These, readily, are the minimum chain decompositions according to the two acyclic orientations.

What is left to note is that only edges leaving $`s`$ or arriving at $`t`$ can be saturated by the maximum flow, so the corresponding minimum cut is necessarily given by a group of such edges. In Figure 2, the minimum cut in each network is indicated by boxes enclosing the nodes other than $`s`$ or $`t`$ that are involved in the cut. Such nodes do necessarily constitute a minimum node cover of the edges that do not involve $`s`$ or $`t`$ in $`D(G,\omega )`$, and consequently induce a node cover in $`G`$ as well ($`\{1,2,3,4\}`$ for the topmost network, $`\{2,3,5\}`$ for the other) with a corresponding independent set ($`\{5\}`$ and $`\{1,4\}`$, respectively).

In general, then, we have the following . If $`F`$ is the value of the maximum flow from $`s`$ to $`t`$ in $`D(G,\omega )`$, then $`nF`$ is the number of chains in the minimum chain decomposition of $`G`$ according to $`\omega `$ and $`F`$ is the size of a node cover in $`G`$. Consequently, $`nF`$ is also the size of an independent set in $`G`$. In order to determine the actual nodes that constitute this independent set, it suffices to look at the minimum cut in $`D(G,\omega )`$ and at the node cover it induces on the edges of $`D(G,\omega )`$ that do not touch $`s`$ or $`t`$. This node cover corresponds to a node cover in $`G`$ as well, whose complement with respect to $`N`$ is then the desired independent set.

We now turn to the three key elements of our formulation, namely how to assess an individual’s fitness and how crossover and mutation operate. We have designed these elements in such a way that they can be used directly in most standard templates of fitness-maximization genetic algorithms . What is novel in our formulation is the adoption of $`\mathrm{\Omega }(G)`$ as the search space out of which populations are formed. Not only do the acyclic orientations of $`G`$ relate cleanly, as we have discussed, to the independent sets of $`G`$, but also they allow for evolutionary operators that are simple and yet effective in several aspects of the evolutionary search. We will come to them shortly.

### Fitness evaluation

For $`\omega \mathrm{\Omega }(G)`$, let $`f(\omega )`$ denote the fitness to be maximized over $`\mathrm{\Omega }(G)`$ while searching for a maximum independent set of $`G`$. By (1), we have

$$f(\omega )=\underset{dD_\omega }{\mathrm{min}}|d|,$$
(2)
that is, the fitness of individual $`\omega `$ is the number of chains in the minimum chain decomposition according to that orientation.

It follows from our preceding discussion that $`f(\omega )`$ can be assessed along the following steps:

1. Construct the flow network $`D(G,\omega )`$.
2. Compute the value $`F`$ of the maximum flow from $`s`$ to $`t`$ in $`D(G,\omega )`$.
3. Let $`f(\omega )=nF`$.

We note, with regard to Step 2, that only the value $`F`$ of the maximum flow is needed for fitness evaluation, not the actual flow. In other words, what is needed is the maximum total flow incoming to $`t`$, not the particular assignment of flows to all edges. When a push-relabel method is used to compute the maximum flow, it is a simple matter to separate the computation into two phases : the first phase computes $`F`$ but may leave excess flow at some nodes; the second phase corrects this by returning flow in order to eliminate excesses. Conveniently, one of the most successful implementations currently available of a maximum-flow algorithm does precisely this , and as such allows the computation to stop at the end of the first phase, right after $`F`$ has been found. What is also convenient is that, when push-relabel methods are thus implemented, at the end of the first phase the minimum cut is also known, which is useful for determining the independent set that corresponds to the best individual found during the evolutionary search, in the manner we indicated earlier in this section. We return to this in Section 4.

### Crossover

The crossover of the two individuals $`\omega _1,\omega _2\mathrm{\Omega }(G)`$ to yield the two offspring $`\omega _1^{},\omega _2^{}\mathrm{\Omega }(G)`$ is best described in terms of a linear representations of the individuals involved. For individual $`\omega \mathrm{\Omega }(G)`$, the representation we adopt, denoted by $`L(\omega )`$, is the sequence $`L(\omega )=i_1,\mathrm{},i_n`$, where $`i_1,\mathrm{},i_n`$ are the nodes of $`G`$. In this sequence, and for $`1x,yn`$, node $`i_x`$ appears to the left of $`i_y`$ (i.e., $`x<y`$) if $`(i_x,i_y)`$ is an edge of $`G`$ and is oriented by $`\omega `$ from $`i_x`$ to $`i_y`$. For example, both $`3,1,4,2,5`$ and $`3,4,1,2,5`$ are valid linear representations of the bottommost acyclic orientation of Figure 1. Clearly, $`L(\omega )`$ represents $`\omega `$ unambiguously, though not uniquely.

Now let $`L(\omega _1)=i_1,\mathrm{},i_n`$ and $`L(\omega _2)=j_1,\mathrm{},j_n`$. Let also $`z`$ such that $`1z<n`$ be the crossover point. Then $`L(\omega _1^{})=i_1^{},\mathrm{},i_n^{}`$ and $`L(\omega _2^{})=j_1^{},\mathrm{},j_n^{}`$, where

* $`i_1^{},\mathrm{},i_z^{}=i_1,\mathrm{},i_z`$;
* $`i_{z+1}^{},\mathrm{},i_n^{}`$ is the subsequence of $`j_1,\mathrm{},j_n`$ comprising all nodes that are not in $`i_1,\mathrm{},i_z`$;
* $`j_1^{},\mathrm{},j_z^{}=j_1,\mathrm{},j_z`$;
* $`j_{z+1}^{},\mathrm{},j_n^{}`$ is the subsequence of $`i_1,\mathrm{},i_n`$ comprising all nodes that are not in $`j_1,\mathrm{},j_z`$.

Notice that the $`L(\omega _1^{})`$ and $`L(\omega _2^{})`$ thus determined are valid linear representations of acyclic orientations, since, by construction, both sequences contain all nodes from $`G`$. Furthermore, $`\omega _1^{}`$ and $`\omega _2^{}`$ inherit edge orientations from $`\omega _1`$ and $`\omega _2`$ as follows:

* Edges joining nodes in the set $`\{i_1^{},\mathrm{},i_z^{}\}`$ to any other nodes are oriented by $`\omega _1^{}`$ exactly as by $`\omega _1`$.
* Edges joining nodes in the set $`\{i_{z+1}^{},\mathrm{},i_n^{}\}`$ exclusively are oriented by $`\omega _1^{}`$ exactly as by $`\omega _2`$.
* Edges joining nodes in the set $`\{j_1^{},\mathrm{},j_z^{}\}`$ to any other nodes are oriented by $`\omega _2^{}`$ exactly as by $`\omega _2`$.
* Edges joining nodes in the set $`\{j_{z+1}^{},\mathrm{},j_n^{}\}`$ exclusively are oriented by $`\omega _2^{}`$ exactly as by $`\omega _1`$.

### Mutation

Like crossover, our mutation operator is defined in terms of the linear representations of individuals. We use single-locus mutation at the mutation point $`z`$ with $`1zn`$. For individual $`\omega \mathrm{\Omega }(G)`$, the mutation operator turns node $`i_z`$ into a source, thus yielding another acyclic orientation $`\omega ^{}`$. That $`\omega ^{}`$ is indeed acyclic has been argued elsewhere : in essence, a directed cycle through $`i_z`$ would be required for $`\omega ^{}`$ not to be acyclic, which is impossible, $`i_z`$ being a source according to $`\omega ^{}`$.

We also refer the reader elsewhere (, Section 4) for an argumentation as to why this type of mutation does more for the evolutionary search than simply to allow occasional random jumps loosely intended to escape local optima. Specifically, what is shown is that, for any two acyclic orientations $`\omega ,\omega ^{}\mathrm{\Omega }(G)`$, there necessarily exists a finite sequence of mutations that turns $`\omega `$ into $`\omega ^{}`$. So this operator can be regarded as providing the search space $`\mathrm{\Omega }(G)`$ with an underlying fundamental connectedness that allows, at least in principle, every acyclic orientation to be reached regardless of where the evolutionary search is started.

## 3 The experimental test set

In this section we give a brief description of the benchmark graphs used in the experiments reported in Section 4. They have all been extracted from the DIMACS challenge suite . That suite is structured from the perspective of finding maximum cliques, so the graph types listed next characterize $`\overline{G}`$.

c-fat$`n`$-$`c`$ . This graph comes from fault-diagnosis problems  and its set of $`n`$ nodes is partitioned into $`n/c\mathrm{log}n`$ sets of approximately equal sizes. Edges are deployed so that every node is connected to every other node in its own set and in the two sets that are neighbors to its own (according to an arrangement of the sets into a ring).

johnson$`W`$-$`w`$-$`d`$ . This graph arises in problems from coding theory. It has $`n=\left(\genfrac{}{}{0pt}{}{W}{w}\right)`$ nodes, each node labeled with a $`W`$-digit binary number having exactly $`w`$ $`1`$’s. Two nodes are joined by an edge if the Hamming distance between their labels (the number of digits at which they differ) is at least $`d`$.

keller$`d`$ . This graph is derived from another with $`4^d`$ nodes which arises in connection with proving Keller’s 1930 conjecture false for high dimensions. The conjecture is that a “tiling” of Euclidean $`d`$-dimensional space by unit cubes necessarily contains two cubes meeting in a full $`(d1)`$-dimensional face.

hamming$`W`$-$`d`$ . This graph comes from coding-theory problems also and has $`n=2^W`$ nodes, each node labeled with a $`W`$-digit binary number. Two nodes are connected if the Hamming distance between their labels is at least $`d`$.

san$`n`$_$`f`$\_x . This graph is artificially constructed on $`n`$ nodes to have a maximum clique whose size is determined beforehand. Its number of edges is $`fn(n1)/2`$; x is only used to differentiate among instances.

sanr$`n`$_$`p`$ . This is a random graph on $`n`$ nodes, generated by adding an edge with constant probability $`p`$ between any two distinct nodes . It is expected to have dimensions close to those of san$`n`$_$`f`$\_x for $`f=p`$.

brock$`n`$\_x . This is a random graph on $`n`$ nodes, designed to have a maximum clique much larger than would be expected from the nodes’ degrees. The number x is used for instance differentiation.

p\_hat$`n`$-x . This is a random graph on $`n`$ nodes whose density is based on two parameters. This contrasts with the usual random graphs with fixed edge probability , leading to node degrees that are more spread and to larger cliques also. The number x differentiates among instances.

MANN\_ax . This graph gives the clique formulation of the instance $`A𝚡`$ of the Steiner triple problem. This formulation is obtained by a conversion from the set-covering formulation of that problem .

Details on the graphs we used in our experiments are given in Tables 1 and 2. For each $`\overline{G}`$, the tables give the values of $`n`$ and $`m`$ (the number of edges in $`G`$), as well as $`\alpha (G)`$, when known from design characteristics.

## 4 Experimental results

Henceforth, we refer as WAO to the algorithm that results from the formulation of Section 2. This denomination is an acronym after “Widest Acyclic Orientation,” as by (1) what is sought during the evolutionary search is an acyclic orientation whose minimum chain decomposition has the most chains over $`\mathrm{\Omega }(G)`$.

WAO iterates for $`g`$ generations, each one characterized by a population of fixed size $`s`$. After generating the last population, it outputs the best individual found during the entire evolutionary search. For $`k>1`$, the $`k`$th population is obtained from the $`k1`$st population as follows. First an elitist step is performed, resulting in the transfer of the $`fs`$ fittest individuals from the current population to the new, with $`0f<1`$. Then WAO performs the following iteration until the new population is full: with probability $`p_c`$, two individuals are selected from the current population and the crossover operator is applied to them, the resulting two individuals being then added to the new population; with probability $`1p_c`$, one single individual is selected from the current population and then is subjected to the mutation operator before being added to the new population.

In order to decide on an appropriate selection method, we ran several initial experiments on reasonably-sized graphs. From these experiments emerged not only the selection method of our choice but also the suite of parameters we would adopt in all further experiments (we discuss these later). The selection method we used in our experiments picks individuals proportionally to their linearly normalized fitness in the current population. For $`1ks`$, this means that the $`k`$th fittest individual—that is, $`\omega `$ such that $`f(\omega )`$ is the $`k`$th greatest—is selected with probability proportional to

$$g(\omega )=L\left(\frac{L1}{s1}\right)(k1).$$
(3)
Ties between two individuals are broken by taking the individual that was added to the current population first as the fitter one. In (3), $`L`$ is the factor by which the linearly normalized fitness of the fittest individual in the current population (the $`k=1`$ case) is greater than that of the least fit individual (the $`k=s`$ case); that is, $`L=g(\omega _1)/g(\omega _s)`$, where $`\omega _1`$ and $`\omega _s`$ are those two individuals, respectively.

We present our results in comparison to those obtained by the heuristics of the DIMACS challenge . This is not to say that the best results known to date are necessarily the ones obtained by those heuristics, since several new methods have appeared in the meantime for the maximum independent set problem under one of its guises (e.g.,  and their references). However, all those more recent methods invariably go back to the DIMACS challenge heuristics as references for comparison, so those heuristics serve as an indirect basis for other comparisons as well. One exception to this comparison rule we have adopted is the genetic-algorithm approach of , which, like WAO, employs a nontrivial crossover operator. We present next a brief description of all the eleven heuristics to which we compare WAO directly.

B&C . This is a branch-and-cut method for which cutting planes are generated based on the more general technique of . It starts with an integer programming formulation of the maximum clique problem, and proceeds from the initial relaxation by generating new cutting planes and incorporating them into the current linear program.

Cliqmerge . This heuristic is based on a procedure that finds a maximum clique in the subgraph induced by the nodes of two cliques when they are merged together. The essence of the method is to find a maximum bipartite matching in the complement of this subgraph.

Squeeze . This is a branch-and-bound algorithm for the maximum independent set problem. Its lower bounds are obtained through a reduction to the problem of minimizing a general quadratic $`0`$-$`1`$ function.

CBH . This is an interior-point approach (cf. ) to the determination of maximum independent sets. Following a continuous formulation of the problem, a relaxation of it is solved and the resulting solution is rounded by a heuristic based on .

RB-clique . This method uses backtracking “coordinates” as the entities on which restricted backtracking is to be applied while seeking a maximum clique. The restrictions to which the backtracking coordinates are subject are given as input.

AtA . This strategy employs recurrent neural networks (cf. ) to find maximum cliques. The crux of the approach is an adaptive procedure for the determination of appropriate threshold parameters and initial state for the neural network.

SA&GH . This is a blend of heuristics to find maximum cliques. In most cases it employs simply simulated annealing , but for very dense graphs the greedy heuristic of  is used.

XSD . This is a family of heuristic methods to find maximum cliques. The methods are all related to neural-network models and include deterministic and stochastic descent approaches, with or without an intervening learning step between restarts.

B&B . This is a branch-and-bound approach to the maximum independent set problem. Upper bounds are obtained through a procedure derived from edge projection, a specialization of the clique projection of .

XT . This is a family of three variants of tabu search , two deterministic and one probabilistic.

OCH . This is a genetic algorithm for the maximum independent set problem. Each individual is an $`n`$-digit binary number, each digit indicating whether the corresponding node is in the independent set or not. The centerpiece of the method is the so-called optimized crossover, which generates one optimal offspring based on the same merging procedure of Cliqmerge and one other having a random character.

All the experiments we report on were conducted with $`g=10n`$, $`s=1.5n`$, $`f=0.05`$, $`p_c=0.2`$, and $`L\{15,30\}`$. As we indicated earlier, these reflect policies and values that emerged from early experiments on reasonably-sized graphs. They by no means represent optimal decisions of any sort, since the number of possible choices is, naturally, far too large.

In all experiments, we also made use of the maximum-flow code of , which implements the algorithm of , to compute the value of $`f(\omega )`$ as explained in Section 2. As we also indicated in that section, it suffices for the maximum-flow computation to stop right after completing its first phase, since the value of the maximum flow, and hence the size of the independent set that the individual contributes, is already known at this point. We also recall that running the maximum-flow code through its first phase only is sufficient even for the best individual found during the whole evolutionary search. In this case, what we need is to enumerate the members of the independent set contributed by that individual, not simply to know its size, but this can be obtained by examining the minimum cut that is also already known as the first phase ends.

Our results are shown in Figures 3 and 4, and in Tables 3 through 8. The two figures show, for ten graphs selected from Tables 1 and 2, the evolution of the best fitness (as given by (2)) ever found for an individual as the generations elapse during the best of twenty independent runs (the one that eventually yielded the largest independent set for that graph). In other words, they show the size of the largest independent set yet identified. These ten graphs were selected because, in terms of what is shown in Figures 3 and 4, they have led the evolutionary search to behave either in a way that we found to be somewhat typical or a way that yields interesting insight. We return to this shortly.

Tables 3 through 8 show the results obtained by WAO alongside the results of the competing algorithms we outlined earlier in this section. Tables 35 refer to the graphs in Table 1, while Tables 68 refer to those in Table 2. For each graph $`\overline{G}`$, the tables give the value of $`\alpha (G)`$, when known from design characteristics, and the sizes of the maximum independent sets obtained by the algorithms on $`G`$, when available. The result reported for WAO on each graph is the best result found over twenty independent runs. For the other methods, the results reported are the best results they yielded, as published in . The number appearing in parentheses next to “WAO” in Tables 5 and 8 is the value of $`L`$ that was used to obtain the results listed in the corresponding column.

The several runs of WAO were executed on a relatively wide assortment of machine architectures, so we refrain from providing detailed timing data. Also, comparing running times to those of the other methods—obtained roughly one decade ago—would be cumbersome however we tampered with the numbers seeking to compensate for the technological gap. Given these constraints, all we do is mention that each of our runs tended to complete somewhere between very few seconds and a week, depending on the graph at hand.

WAO is a competitive method, by all that can be inferred from Tables 38. We have in the tables used a bold typeface to indicate the instances on which WAO performed at least as well as the best performers. Several entries are thus marked, and for many that are not WAO is seen to have fallen short by a very narrow margin.

Returning to the plots in Figures 3 and 4 may highlight some of the patterns that help explain success or failure at meeting the best performers’ figures during our experiments with WAO. Of the ten graphs to which those figures refer, three correspond to cases in which WAO missed by a narrow margin (san$`400`$_$`0.9`$\_1, sanr$`400`$_$`0.7`$, and brock$`200`$\_4) and one to a case in which it missed widely (brock$`400`$\_4). What seems to distinguish one group from the other is that in the latter case the plot becomes flat early in the evolution, perhaps signaling an inherent hardness at escaping some particularly difficult local maximum. For the three graphs in the former group, however, and in fact for the six graphs in the figures outside either group, evolution seems to lead to fitness growth more or less steadily along a comparatively larger number of generations, even though for the san graphs it nearly stalls for a significant number of generations before it gains momentum again. In the case of san$`400`$_$`0.9`$\_1, particularly, it appears quite likely that a few more generations would have bridged the narrow gap between the $`98`$ that WAO achieved and the $`100`$ of its best contenders.

## 5 Concluding remarks

We have in this paper introduced WAO, a novel evolutionary heuristic for the maximum independent set problem. WAO is based on a view of the problem that relates independent sets to the acyclic orientations of the graph, and seeks to identify an acyclic orientation that is widest (has the decomposition into the fewest number of chains that requires the most chains) over the set of all the acyclic orientations of the graph. It incorporates no additional sophistication into the usual evolutionary-algorithm machinery, but rather into the design of the individuals’ representations and the evolutionary operators, all based on complex graph-theoretic notions.

We have found our new heuristic to perform competitively when compared to several others on the DIMACS benchmark graphs. Notwithstanding this, there certainly is room for further investigation and improvements. For example, there may exist a better set of parameters for the evolutionary search, just as there may exist a better functional dependence of $`g`$ and $`s`$ on $`n`$ (and perhaps also on $`m`$, unlike what we adopted in our experiments). Likewise, it is also conceivable that the formulation itself may be improved by the incorporation of optimizations into the crossover or mutation operator, or even by the introduction of new operators.

### Acknowledgments

The authors acknowledge partial support from CNPq, CAPES, the PRONEX initiative of Brazil’s MCT under contract 41.96.0857.00, and a FAPERJ BBP grant.


text:
# Trapped surfaces and symmetries

## Acknowledgements

We thank the anonymous referee for useful comments. JMMS acknowledges financial support under grants BFM2000-0018 of the Spanish CICyT and no. 9/UPV 00172.310-14456/2002 of the University of the Basque Country. MM wishes to thank the Junta de Castilla y León for financial support under grant SA002/03

## References


text:
# Simulations in the ϵ-Regime of Chiral Perturbation Theory Poster presented by K.-I. Nagai at Lattice 2003 Preprint DESY 03-126, HU-EP-03/56, SFB/CPP-03-31

## Abstract

We discuss the potential of Ginsparg-Wilson fermion simulations in the $`ϵ`$-regime of chiral perturbation theory, regarding the determination of the leading low energy constants of the effective chiral Lagrangian. It turns out to be very hard to measure observables in the topologically trivial sector. There a huge statistics would be required, due to the frequent occurrence of very small eigenvalues. Moreover, contact with chiral perturbation theory is only established if the physical volume of the system is sufficiently large.

QCD at low energy can be described by an effective Lagrangian $`_{\mathrm{eff}}`$ for the quasi-Goldstone bosons of chiral symmetry breaking. $`_{\mathrm{eff}}`$ may then be evaluated by chiral perturbation theory ($`\chi `$PT). The case of a finite volume and very light quarks — so that the pion Compton wave length clearly exceeds the box length — is called the $`ϵ`$-regime. There the zero-mode is treated by collective variables, while the excitations (which do fit into the box) can be evaluated perturbatively in the $`ϵ`$-expansion . In this framework, observables strongly depend on the topology of the gauge background , hence they should be measured in fixed topological sectors.

$`_{\mathrm{eff}}`$ involves coupling constants which appear as free parameters in $`\chi `$PT. In principle they can be determined by comparison with numerical data from lattice QCD in the $`ϵ`$-regime. The same values are also relevant for the physical situation in a large volume (the standard $`\chi `$PT), hence their knowledge is very important in view of experiments and further numerical studies, especially for the chiral extrapolation. Their numerical evaluation, however, may be easier in the $`ϵ`$-regime, due to the usage of a small volume.

Such simulations in the $`ϵ`$-regime are now conceivable by using Ginsparg-Wilson fermions, which avoid additive mass renormalization and define the topological charge through the Index Theorem . Here we discuss the question how one could extract results for the low energy constants which occur in the leading order of $`_{\mathrm{eff}}`$, the pion decay constant $`F_\pi `$ and the scalar condensate $`\mathrm{\Sigma }`$.

We performed quenched QCD simulations with the Wilson gauge action and the overlap Dirac operator

$`D_{\mathrm{ov}}=\left(1{\displaystyle \frac{m}{2\mu }}\right)D_{ov}^{(0)}+m,`$ (1)
$`D_{\mathrm{ov}}^{(0)}=\mu (1+A/\sqrt{A^{}A}),A=D_W\mu ,`$
where $`m`$ is the quark mass and $`D_W`$ the Wilson operator. At $`\beta =6`$ and $`\beta =5.85`$ we set $`\mu =1.4`$ resp. $`\mu =1.6`$. We approximated the inverse square root by Chebyshev polynomials to an accuracy of $`10^{12}`$. To identify the low lying eigenvalues (EVs) of $`D_{\mathrm{ov}}^{(0)}D_{\mathrm{ov}}^{(0)}`$ we applied alternatively the Ritz functional method and the Arnoldi algorithm. The latter allows for the determination of many EVs ($`O(100)`$). The very lowest EVs determine the index $`\nu `$: there is usually a gap by several orders of magnitude between the numerical values of the zero EVs and the rest of the spectrum. The chirality of the zero modes reveals the sign of $`\nu `$. Two typical examples are given in Table 1.

On a $`10^3\times 24`$ lattice at $`\beta =6`$ we also compared the index with the topological charge determined from (standard) cooling. In this respect, we investigated 51 configurations: for 41 of them the index of $`D_{\mathrm{ov}}`$ and the charge determined from cooling coincided, and in the remaining 10 cases it differed by 1.

Next we mention that Random Matrix Theory (RMT) has been applied to QCD and it yields predictions for the probability distributions of the low lying eigenvalues . For the first non-zero EV these predictions are confirmed if the physical volume is sufficiently large, $`V>(1.2\mathrm{fm})^4`$ . In particular, we know from these probability distributions that in the topologically trivial sector there is a significant density of very small (non-zero) EVs. As an example, we found a configuration on the $`10^3\times 24`$ lattice at $`\beta =6`$ where the lowest EVs of $`D_{\mathrm{ov}}^{(0)}D_{\mathrm{ov}}^{(0)}`$ read: $`7.32e5`$, $`7.33e5`$, $`1.04e4`$ etc. The corresponding chiralities $`0.127`$, $`0.127`$, $`0.826`$ show clearly that this configuration has index 0. This is also confirmed from cooling, but the form of the cooled configuration suggests that charge 0 arises from the cancellation of a pair of topological objects with charges $`\pm 1`$, see Fig. 1. The index never involved such a cancellation. However, in the case of small non-zero EVs there seems to be a trend towards a cancellation picture from cooling, as in the example discussed here.

We now address the meson correlators as physical observables, from which one can try to extract $`F_\pi `$ and $`\mathrm{\Sigma }`$. The corresponding formulae for quenched $`\chi `$PT have been worked out in Ref. . For a first comparison to numerical data we refer to Refs. . The axial current correlation function on an $`L^3\times T`$ lattice,

$`A_\mu (t)A_\mu (0)_\nu ={\displaystyle \frac{F_\pi ^2}{T}}+2m\mathrm{\Sigma }_{|\nu |}(z)Th_1(\tau )`$ (2)
$`2h_1(\tau )=\tau ^2\tau +1/6,\tau =t/T`$
$`{\displaystyle \frac{\mathrm{\Sigma }_\nu (z)}{\mathrm{\Sigma }}}=z\left[I_\nu (z)K_\nu (z)+I_{\nu +1}(z)K_{\nu 1}(z)\right]+{\displaystyle \frac{\nu }{z}}`$
is most suitable for our purpose, because it only involves $`F_\pi `$ and $`\mathrm{\Sigma }`$ as free parameters. The determination of $`F_\pi `$ is relatively easy in the minimum at $`\tau =1/2`$. $`\mathrm{\Sigma }`$ is related to the curvature and far more difficult to extract, since the sensitivity of this curvature to a change of $`\mathrm{\Sigma }`$ is hardly visible over a wide range 111This situation is somehow complementary to the RMT comparison of the EVs, where $`\mathrm{\Sigma }`$ is easier to determine .. However, we saw also here that the physical volume must exceed the lower limit given before in the context of RMT, otherwise one does not obtain the curvature required by any acceptable value of $`\mathrm{\Sigma }`$. Once we are in the right regime, $`\nu =0`$ may look like the simplest case from eq. (2), but from the numerical point of view this is a nightmare. The problems in this sector are related to the danger of very small EVs. Ref.  pointed out before that they may be a serious problem, and our observations confirm this in a striking way. As an example, we show histories for the measurement of $`A_0(t)A_0(0)_0`$ and $`A_0(t)A_0(0)_1`$ in Fig. 2. In the topologically neutral sector there are strong spikes at those configurations which have very small EVs. Their height is maximal at very small quark mass, as Fig. 3 illustrates.

Obviously the occurrence of such spikes makes it very hard to measure expectation values — they stabilize only with a huge statistics. To get an estimate for this effect, we consider the contribution of the smallest (non-zero) EV alone to the scalar condensate. This contribution reads

$$\mathrm{\Sigma }_{\mathrm{min}}^{(\nu )}=\frac{1}{V}_0^{\mathrm{}}𝑑zP_\nu (z)\frac{2m}{m^2+(z/\mathrm{\Sigma }V)^2},$$
(3)
where we insert the probabilities $`P_\nu (z)`$ provided by RMT. We now generated fake values for the lowest EV with probability $`P_\nu (z)`$ and computed $`\mathrm{\Sigma }_{\mathrm{min}}^{(\nu )}`$ from them. The result can only be trusted when the standard deviation becomes practically constant. Fig. 4 shows the evolution of this standard deviation as the statistics is enhanced. We see that the sector $`\nu =0`$ would require a tremendous statistics of $`O(10^4)`$ configurations, but the situation is clearly better for $`\nu 0`$.

We conclude that there are several conditions for running conclusive simulations with Ginsparg-Wilson fermions in the $`ϵ`$-regime. The volume should obey $`V>(1.2\mathrm{fm})^4`$, and measurements have to be performed in a sector of non-trivial topology. The quark mass should be small for conceptual reasons, but taking it too small causes technical problems in the measurements.


text:
# The Nuclear Region of Low Luminosity Flat Radio Spectrum SourcesBased on observations collected at the Multiple Mirror Telescope on Mt. Hopkins.

## 1 Introduction

BL Lacertae objects (BL Lacs), together with Flat Spectrum Radio Quasars (FSRQs) are called blazars. Blazars display some of the most extreme behaviour observed in Active Galactic Nuclei (AGN), such as high and variable optical polarization; they also show a strong and rapid variability at radio, optical and X-ray wavelengths. BL Lacs usually have weak or non-observed optical emission lines.

Blazars exhibit characteristics indicative of relativistic beaming (Padovani & Urry 1990), which occurs when the relativistic emission (a jet) is oriented close to the observer’s line of sight (Blandford & Rees 1978). The extreme properties of blazars are generally interpreted as a consequence of non-thermal emission from the relativistic jet. A transition population between beamed BL Lacs and unbeamed radio-galaxies has not been detected (Rector et al. 1999); yet, a transition population of “low luminosity BL Lacs” (LL BLLs) was predicted to exist in abundance in samples like the one studied by Browne & Marchã (1993).

In the past few years there has been a growing interest in LL BLLs and other flat radio spectrum sources. Because they are relatively close to us, their host galaxies can be studied more easily. The knowledge of their host galaxies is of the utmost importance, not only as a test of the unification schemes, but also as a way to better understand the relationship between nuclear activity and host galaxy. BL Lac host galaxies seem to be luminous ellipticals, often in loose groups; in fact, from a morphological point of view, BL Lac host galaxies seem to be absolutely normal ellipticals (Urry et al. 2000; Falomo & Ulrich 2000). HST observations by Verdoes Kleijn et al. (1999) suggest that, apart from the presence of dust and gas distributions, the nuclei of FR I radio-galaxies also resemble those of normal bright ellipticals. While there are some studies concerning the host galaxy morphology, little is, however, known about the nuclear stellar populations of Low Luminosity Flat Radio Spectrum (LL FRS) sources.

In an effort to better understand such objects and their host galaxies, we undertook the study of the optical properties of a sample of low luminosity, core-dominated radio sources expected to contain a high fraction of objects observed at small angles to the line-of-sight. This sample is known to contain several weak emission line galaxies, which could be objects related to BL Lac phenomena, but observed at larger angles from the line-of-sight; it is also possible that some of the galaxies are “hidden” BL Lacs, whose nuclear emission is completely swamped by the host galaxy starlight (Dennett-Thorpe & Marchã 2000).

With this study, we had two main purposes:

(i) studying the nuclear stellar populations of low luminosity, flat radio spectrum nuclei such as LL BLLs;

(ii) to search for the clear optical signature of the underlying activity.

To do this, we must rely on the computation of the nuclear stellar populations with subsequent subtraction of their contribution to the total nuclear spectra. Our stellar synthesis is not done by using templates of normal nearby galaxies to estimate the stellar populations. In fact, we take the spectra of our objects to compute the stellar populations; in this way, the information on the stars that compose our galaxies is taken directly from the data, without any a priori assumption on the nature of the stellar populations.

In Sect. 2 we describe the sample from which we have selected the objects discussed in this paper; we present our own sub-sample, as well as the criteria behind its selection. Section 3 briefly describes the observational data analyzed in the paper. In Sect. 4 we introduce the stellar population synthesis method; we also describe how it has been applied to the data, presenting the results obtained for each galaxy. In Sect. 5 we summarize our results and conclusions.

## 2 The sample

Being interested in studying the optical properties of LL FRS sources, we undertook the analysis of a set of such objects taken from the sample discussed in Marchã et al. (1996). This sample presents objects with a redshift up to 0.1. Since it only contains nearby objects, the host galaxies are therefore easy to study.

Marchã’s et al. original sample contains 57 optically bright ($`V`$ 17), nearby ($`z`$ 0.2), radio sources; these were selected on the basis of the spectral slope ($`\alpha _\mathrm{r}`$ 0.5 between 1.4 and 5 GHz) and are all located north of 20°, being off the galactic plane by at least 12°. The sample has a flux density limit of 200 mJy at 5 GHz.

The main purposes of such a sample were:

(i) to find LL BLLs with radio luminosities comparable to X-ray selected BL Lacs (most of which we now refer to as “high frequency peak” BL Lacs);

(ii) to investigate the differences between BL Lacs and other flat spectrum radio sources.

Spectroscopy of these objects yielded a high fraction of BL Lacs ($``$ 35%) and a relatively high number ($``$ 20%) of sources with strong emission lines. The remainder of the sample is composed of objects with hybrid properties (i.e. objects showing simultaneously BL Lac properties, such as high and variable polarization, and broad emission lines, usually associated with the spectra of quasars) and by galaxies which do not show any conspicuous signs of activity, the main contribution to the their optical spectra being of stellar origin.

In this paper we focus our attention on the latter. These galaxies, although radio selected together with confirmed BL Lacs and objects displaying a “Seyfert-type” emission spectrum, following exactly the same criteria, are optically very different from them, showing no strong emission lines (as Seyferts do), nor a flat continuum and no emission or absorption features (as usually observed in BL Lacs). This persuaded us that there should be more to these sources than meets the eye.

Among the galaxies in Marchã’s sample which displayed such characteristics, we have selected 19 objects according to the following criteria:

(i) the dilution of the absorption features due to the power-law continuum must be small enough so that information can be inferred from the lines. We took the 4000 Å break contrast $`C`$ = 0.3111The 4000 Å break contrast $`C`$ is defined as $`(S^+S^{})/S^+`$, where $`S^+`$ and $`S^{}`$ are the fluxes measured redward and blueward of the break, respectively. Several tests with our population synthesis method have shown that for lower values of $`C`$ the stellar populations cannot be determined. as our lowest possible limit.

(ii) the signal to noise ratio must be high enough so that the continuum can be traced and the absorption features used.

Table 1 gives the list of the 19 galaxies selected for this study; for the sake of consistency, we have kept the same nomenclature as in Marchã et al. (1996).

## 3 Observational data

The data consist of optical spectra ($``$ 3500–8700 Å and $``$ 3500–7400 Å) obtained at the Multiple Mirror Telescope on Mt. Hopkins. The spectra were acquired through a 1$`\stackrel{}{.}`$5 slit, oriented along the parallactic angle, and have a spectral resolution of 20 Å (FWHM); they have been flux-callibrated using the spectrophotometric standards observed on each night. All spectra were reduced according to the standard procedure using IRAF and ESO/MIDAS packages. Further observational and reduction aspects are discussed in detail in Marchã et al. (1996).

The spectra have been corrected for the interstellar reddening using the Howarth (1983) Galactic reddening law; the values of $`E(BV)`$ for each galaxy were calculated using the maps of dust IR emission from Schlegel et al. (1998). A correction to the broadening of the lines due to velocity dispersion is not necessary, since our spectral resolution is low.

22footnotetext: The break contrast values given by Dennett-Thorpe & Marchã (2000) are slightly different than those calculated by Marchã et al. (1996); this is because they use the $`S_\nu `$, rather than the $`S_\lambda `$ flux. The selection of our targets is not affected by this difference.

## 4 Stellar population synthesis

The main purpose of a stellar population synthesis is to gather information on the stellar content of the galaxies under study; our synthesis also provides an evaluation of the dust content, the internal reddening being a parameter taken into account for the calculations.

In addition to these pieces of information, the stellar synthesis also allow us to remove the starlight contribution from the data and consequently to unveil the optical emission-line features present in the spectra. Even though a few weak emission lines could be previously detected in some of the galaxies in our sample, their optical to near-UV spectra are clearly dominated by the stellar continuum. In fact, the fraction of starlight in such objects is so important that it almost completely swamps any strong emission feature associated with the active nucleus.

The method used to compute the nuclear stellar populations and its application to our sub-sample of objects are described below.

### 4.1 The method

We have used a synthetic population algorithm developed by Pelat (1997). This algorithm uses a new mathematical method which gives a unique solution (Global Principal Geometrical solution, or GPG), contrary to the other methods widely used for population synthesis. It uses the equivalent widths (EW) of all the absorption features found in the spectra. Essentially, it considers a galaxy as being made up of a set of stars with different spectral types, luminosities and metallicities. This particular composition will carry its own signature in terms of the EW of the absorption lines; it then defines the composition by matching this signature as closely as possible to the observed EW, using the least squares method. The internal reddening is also a free parameter given by the method in an indirect way.

The accuracy of the fit is estimated by a parameter, the distance, which is defined as the sum, for all absorption features, of the difference between the observed EW and the synthetic EW calculated through the combination of the different stars that compose the stellar populations obtained. This value depends on the number of absorption features used and the smaller this distance, the better the fit. In addition, residuals estimated over the continuum help us to verify the accuracy of the solution found. The theoretical aspects of the method are described in Pelat (1997). This method was already applied to the synthesis of nuclear stellar populations in AGN and proved to give very good results (Serote Roos 1996; Serote Roos et al. 1998; Boisson et al. 2000).

The stars used for the synthesis were compiled from the Silva & Cornell (1992) library. These authors have made available 72 stellar types with a large coverage in spectral types, luminosity classes and metallicities. The library wavelengths range from 3500 Å to 9000 Å, thus covering quite well the spectral range of our data. Their resolution being 11 Å, we have convolved the stars in the library with a Gaussian with $`\sigma `$ = 7.1, in order to match the resolution of our spectra (20 Å). The stellar library has been chosen in order to cover the space temperature/gravity as much as possible without being degenerate. In fact, we cannot include as many stellar types as we would like to, in order to prevent stellar library degeneracy, i.e. having two different stellar types with spectral energy distributions similar enough to be indistinguishable in a mathematical sense. Some super metal rich stars (SMR) have also been included. In total, the stellar library used contains 30 stars.

We have identified and measured the EW of all absorption features ($``$ 80) present in the spectra, as well as the values of the continuum for each feature. Identification of the absorption lines and bands, together with the wavelength ranges used to measure the EW of the features, can be found in Table 2.

The wavelength ranges have been defined taking into account the shape of the absorption features in both hot and cool stars of the stellar library. The continuum level has been determined globally over the whole wavelength range. The error due to the uncertainty on the continuum level is dominant over all other measurements and statistical uncertainties. For well defined strong features (e.g. Na ID) this error is always less or equal to 1 Å in absolute value. It can go up to a few Angström for wide bands and strong blends.

The EW of the same spectral features have been measured in the spectra of the observed galaxies. For these, we discarded the intervals corresponding to the emission lines, as well as those showing atmospheric absorption features. We would like to stress the fact that discarding some intervals does not affect the results of the method, as we still have used a very large number of parameters (between 44 and 63 EW, depending on the objects).

### 4.2 Synthesis results

#### 4.2.1 Stellar content

Several studies performed in normal, nearby and distant ellipticals, in the field and in clusters (e.g. Pickles 1985; Carter et al. 1986; Couture & Hardy 1988, 1990, 1993; Munn 1992; Davies et al. 1993; Brown 1996; Trager 1997; Charlot et al. 1998; Mobasher & James, 2000; Minniti 2001) show old stellar populations, with a main contribution of late stellar types, typically main sequence stars. However, for the more luminous ellipticals a non-negligeable component of O-B main sequence stars seems to be present in many objects. The populations in the centre of the galaxies are generally super metal rich dwarf and/or subdwarf stars, contrary to the ones found in the bulge. A population gradient is usually found, with the populations being younger and richer in metals at the centre of the galaxies. Concerning galaxies in clusters, the younger populations of red giants and supergiant stars are usually found in the ellipticals located at the outskirts of the cluster.

In Fig. 1 we plot the observed spectra (black lines) superposed on the synthetic ones (grey lines); all spectra are given in the rest frame.

Table 3 shows the stellar populations obtained for the 19 objects studied, together with the value of the distance for each solution, the number of absorption features used in the synthesis, and the internal reddening derived.

Due to the nature of the synthesis method we use, it is important to have in mind several considerations when analysing the results obtained for each object. In particular, we do not pretend to derive the stellar content of the objects in its very details, but to have an idea of its main components in terms of hot and cold stars, young and old populations; note that the young population is given not only by hot stars, but also by supergiants.

In a general way, we find old stellar populations typical of normal elliptical galaxies. For more than half of the objects in our sample, we do not find much dust. Concerning the super metal rich population, we have made some tests in order to estimate the importance of this component. Except for 0116+319 and 1241+735, this population is not important; we find good solutions to fit the spectra without it, as well. The hot stars are never present in the solutions, the populations being dominated by G, K and M spectral types. Thus, the younger populations appear in the form of supergiant stars (typical ages of the order of $`10^7`$ years). Despite supergiants being a very special class (which can include various types of stars, from the very massive and very young ones to the AGB or post-AGB stars, i.e. very old stars in their final phases), the fact that we are observing them in the optical range signifies that the contribution from supergiantes we are seeing comes from the young stars (Bruzual & Charlot, 1993); in fact, AGB and post-AGB stars contribute mainly in the UV range.

For several objects in the sample (0733+597, 0848+686, 1144+352, 1146+596 and 1755+626) we find only a small contribution of supergiants, implying that, in these objects, the population is quite old with ages of the order of $`10^9`$$`10^{10}`$ years. For 0651+410, 0902+468, 1245+676, 1658+302 and 2320+203, the largest contribution comes from main sequence and supergiant stars, the contribution from giant stars being quite small ($``$ 10%); this means that the old component comes from a non-evolved population of main sequence stars. For 0651+410 and 1658+302, a strong (41%) young component in the form of supergiants appears. A lack in main sequence stars (contribution of $``$ 10%) is found for two of the galaxies in the sample: 1241+735 and 2202+363; in these cases, the population is distributed more or less equally between giant and supergiant stars.

NGC 3379, a normal elliptical was studied using the same method of computing the stellar populations (Boisson et al., 2000). For this galaxy we have found a very homogeneous old population (with almost no population gradients in the central regions of the galaxy) with a main component of main sequence stars (K and M spectral types) and some contribution from giant (K, M) stars. This population seems to be supermetal rich, i.e. composed of stars with a metallicity higher than solar. This metal enrichment is a direct proof of strong star formation in the past. Little dust is present in the solutions.

Comparing these results with the ones found in the present work, for our sample of galaxies, we verify that our objects have equally an old stellar population, the main contributors being G, K and M stars of luminosity classes V and III. The main difference is that our populations do not seem to be metal rich (except in the case of two objects: 0116+319 and 1241+735), contrary to the one found for NGC 3379. This can be due to the fact that only the central regions of NGC 3379 are discussed in Boisson et al.’s study, whereas in this work we are dealing with an integreted spectrum of the whole galaxy. Since star formation can be really intense in the central regions, this metal enrichment would be more proeminent there, thus justifying the difference found.

#### 4.2.2 Colours and ages

In order to have another constraint on the population ages, we have computed the $`UBVRI`$ colours of our objects. Only the $`BV`$ and $`VR`$ indicies were calculated; we did not use the $`U`$ and $`I`$ colours because these filters do not fall entirely in our wavelengh range, and therefore we would be losing too much light to achieve confident results. Using evolutionary synthesis models for elliptical galaxies, from the Bruzual & Charlot code (Bruzual et al. 1997), we have estimated the ages of the populations. Two different kinds of models were used: a single stellar population model, i.e. a model computing the spectral evolution of a single initial instantaneous-burst of star formation, and a model that computes the evolution of the stellar populations with longer time scales of star formation, i.e. where the star formation rate (SFR) decreases exponentially with time (that is with the available gas content), during 1 Gyr. For both models the results are essentially the same. We find that the colour indices of most spectra are compatible with an intermediate/old stellar population with ages ranging from $`10^9`$ to $`10^{10}`$ years. This is in agreement with the results of our population synthesis, corroborating the method and our conclusions.

#### 4.2.3 Dust content

A degenerency exists between the galaxy’s internal dust and the blue continuum of the spectrum, i.e. we can fit the observed continnum either with a given amount of dust plus a given blue slope, or just with a certain amount of dust (smaller in this case). This happens because both contributions have the effect of blueing the spectrum. In this work we have fitted all spectra in such a way that the continuum is also perfectly fitted by a value of $`E(BV)`$, thus assuming that all contribution to the continuum comes from the dust and stellar content of the host galaxy. Instead, we could have found a blue excess in the continuum by adding more dust to the solutions and then argue about this quantity coming from the active nucleus. However, in this case it would be difficult to quantify either the quantity of dust present in the object, or the slope of the blue continuum. Also, and from the values of the Balmer break contrasts of the objects, we estimate the non-stellar continuum to be negligeable, well within the error bars of the flux calibration of the spectra. This is why we chose not to include a contribution from an AGN to the continuum.

## 5 Summary and Conclusions

In an effort to better understand LL FRS sources and their host galaxies, we undertook the spectroscopic study of 19 objects selected from Marchã’s et al. 200 mJy sample; this was known to contain a high fration of objects oriented close to the line-of-sight, such as LL BLLs. This study focused on the stellar content of the host galaxies.

We performed stellar population synthesis aimed at identifying the nuclear stellar populations present in the objects; our work made use of a very reliable mathematical method, yielding a Global Principal Geometrical solution.

Our main result concerns the stellar populations found in these objects, which are in agreement with those found in elliptical galaxies; our results show that the nuclear populations are composed of old stars of solar metallicity or lower; the dust content is weak, which is also typical of “normal” elliptical galaxies. It is most interesting to note that similar stellar populations were reported in Low Ionization Nuclear Emission-Line Regions. In particular, Serote Roos (1996) and Boisson et al. (2000) have shown, in a study of the stellar populations in a sample of AGN with different types of activity, that LINERs present the more evolved populations of the sample, with little recent star formation, i.e. the central regions of these objects are largely dominated by a red evolved population of metal rich stars.

The present work demonstrates the importance of stellar population synthesis in the study of low luminosity radio sources. We have shown that stellar synthesis allows us to obtain information on the stellar populations of the host galaxies, therefore providing material for further studies on the connection between host galaxy and the active nucleus. Stellar synthesis also constitutes a precious tool allowing us to unveil the optical emission lines present in the spectra of our objects. The analysis of such features provides information on the nature of the nuclear emission regions and excitation mechanisms at work in LL FSR sources; such a study has been performed and is the subject of a companion paper (Paper II).

Acknowledgements: We are greatly indebted to D. Pelat for the use of the population synthesis GPG program. We thank Thibault Lejeune for kindly computing the colours for our objects as well as the evolutionary synthesis models used in this work. We thank M. Marchã and A. Caccianica for fruitful discussions, and P. Véron for valuable comments. M. Serote Roos acknowledges financial support from FCT, Portugal, under grants no. BPD/9995/96 and BPD/5684/01, and from FCT/ESO Project ref. PESO/C/PRO/1254/98. A. C. Gonçalves acknowledges financial support from FCT, Portugal, under grant no. SFRH/BPD/9422/2002.