text:
# Evolution and Ageing

## 1 Introduction

Talking about biological ageing may be somewhat distressing, specially for those who are over forty, as are two of the authors. However, ageing is just one of the features of evolution, and to understand such an unavoidable mechanism it is necessary first to understand the paths of evolution. In fact, the concepts of evolution and ageing are connected by Darwin’s theory of selection of the fittest, published in 1859. Insofar that it concerns the general evolution mechanism, this famous theory states that: If genetically distinct individuals compete for limited resources, those more fitted to the environment will produce more offspring. Random mutations mix the genes, giving rise to new genetic combinations, and at every generation natural selection eliminates the less efficient ones, in order to continuously improve adaptation. However, since ageing is related to an age structure, the same theory concerning this subject is better enunciated as: A mutation endangering the life of an organism before the reproductive age is much more dangerous to the species than a mutation affecting it only late in life, when it has already produced enough offspring to warrant the perpetuation of its lineage. In both cases, the notion of fitness, one of the most debated quantities in population genetics, first introduced by Fisher  and Wright , seems to be closely related to the reproductive rate. Selection is also related to competition, which implies that the weaker are supposed to die. Thus, nature must offer a death mechanism, which turns to be the ageing process.

Perhaps evolution can be fully explained by Darwin’s theory. There is, however, a problem: where is the equation that would allow biologists, geneticists, mathematicians and even exotic physicists to understand how nature works and to publish more and more papers? Because such an equation doesn’t exist, different models have appeared in order to explain the origins of life and its evolution. According to Luca Peliti , these models can be divided into three groups. The first one concerns microevolution, that is, the evolution of individuals belonging to the same species or to closed ones. One example is the Eigen model for quasispecies , described in section 2. In such models the interaction among individuals is generally introduced through some mechanism of global competition.

The second group concerns coevolution, where two or more species interact strongly in such a way that the survival of one species depends on the survival of the other. The most common problem studied with these models is the prey-predator one, presented in a quantitative way by Lotka  and Volterra  many decades ago. Another example is the host-parasite interaction, which will appear later in this paper (section 4.4) as one possible explanation for the evolution of sex.

Finally, the third group corresponds to models for macroevolution or large-scale evolution, that deal with all species alive at the same time but with no particular interacting mechanism between them. The concept of fitness for such models cannot be simply related to the reproduction rate of the individuals, since each species has its own reproductive strategy. It is necessary to consider instead a continuously evolving fitness landscape, that changes whenever a species mutates or disappears. In order to survive, the remaining species are also continuously evolving, trying always to be close to the peaks of the landscape. The Bak-Sneppen model  is one of the simplest and most famous models of this third group. It considers the dynamics of large-scale evolution as a result of a self-organized critical process, responsible for the presence of scaling laws in macroevolutionary data.

If it is easy to group the different models for evolution, the same is not true when the subject is ageing. In 1990 the Russians gerontologist Zhores Medvedev  classified more than 300 theories constructed in order to explain the reduction of survival probability with advancing age. Fortunately, these theories can be grouped in two types: the theories of “Why” and the theories of “How”. The first type tries to understand ageing through a global perspective and to explain its disparities among the species. Theories of the second type search for the specific mechanisms or immediate causes of ageing, each one with a given degree of validity depending on the species.

The theories of “Why” are also divided in two groups. The first one attributes physiological unavoidable reasons for ageing, such as the harmful action of the oxygen radicals that are constantly produced in our bodies or the existence of a programmed cell death after a given limited number of cell divisions. The rate of metabolism of each species is probably the oldest example of this kind of theory. Since large mammals such as elephants move slowly and live much longer than small ones, which are very active and so have a high metabolic rate, for many years it was believed that this was the explanation for the disparities in the longevity of different species. However, it was discovered that there are branches of the same species that hibernate in the northern hemisphere and live as long as those branches that inhabit the southern hemisphere. There are also some kinds of birds that have an extremely long lifetime, despite the large amount of energy they need for flying.

Observations like the ones above have gradually increased the importance of the second group of the “Why” theories, that is The Evolutionary theories, proposed around 1950 by Peter B. Medawar, George C. Williams and other famous biologists (for a review on the history of ageing theories and many of its interesting features see the special issue of La Recherche, 322, July/August 1999). These are theories based in Darwin’s proposal of selection of the fittest . As mentioned before, an important ingredient for modeling using the evolutionary approach is to have an age structure with no reproduction during youth. One of the pioneers, and one of the simplest models using this strategy is the Partridge-Barton model , which considers only two age intervals, one from t=0 to t=1 for juveniles and another, from t=1 to t=2, for adults. Reproduction is possible only at ages 1 and 2 and is followed by death after age 2. However, when only deleterious inherited mutations are considered, this model leads to population meltdown, that is, the population dies out due to the accumulation of such mutations. The Penna model for biological ageing , described in section 4, is also based on the mutation accumulation hypothesis and is now by far the most used one to reproduce and understand different aspects of real population dynamics. It deals with many age intervals and has successfully explained the catastrophic senescence of salmon; why women live longer than men; why does menopause exists, and other biological phenomena. Besides giving results that are in agreement with the empirical Gompertz law of an exponential increase of mortality with age, it doesn’t lead to the population meltdown mentioned above.

This review is organized in the following way: the Eigen model is described in section 2, the Bak-Sneppen model in section 3 and the Penna model and its main results in section 4. In section 5 we present our conclusions.

## 2 The Eigen Model for Microevolution

The Eigen’s quasispecies model  is one of the archetypical representation of the Darwin’s postulates of evolution by natural selection. The model, originally developed to address the issue of explaining the origin of life on Earth, describes the dynamics of populations of replicating biological macromolecules under the influence of selection and mutation mechanisms. It predicts that if mutations occur frequently enough, the target of selection will no longer be a single individual, but rather an ensemble of genetically related individuals called a quasispecies. Moreover, when the mutation rate surpasses a critical value, known as the error threshold, the result is a complete loss of this polymorphic genetic structure. In this section we shall restrict ourselves to Eigen’s first theoretical model of molecular evolution based on deterministic chemical kinetic theory (see  for a complete review of the original formulation). Its conceptual elements and results may nevertheless be fairly well described, and perhaps even more appropriately, in the language of population genetics , in terms of the frequencies of haploid multi-locus individuals (see  for an excellent review on this subject).

### 2.1 The Chemical Ansatz

The original formulation of Eigen focus on a well-defined model system: the flux reactor. It comprises basically a reaction vessel, in which biological macromolecules are continually built up out of energy-rich $`T`$ monomers (triphosphates), that are required for macromolecular synthesis, and decay, after a certain time, back to their energy-deficient $`M`$ monomers (monophosphates). Furthermore, it is assumed that this system can exchange energy and matter with its surroundings, by regulating both the supply of energy-rich and energy-deficient monomers, as well as the total population of macromolecules. Each chemical component of this reaction system consists of a reproducing macromolecule modeled as a single string of $`L`$ digits $`I_i=(s_1^i,s_2^i,\mathrm{},s_L^i)`$, with the variables $`s_\alpha ^i`$ $`(i=1,\mathrm{},\kappa ^L;\alpha =1,\mathrm{},L)`$ allowed to take on $`\kappa `$ different values. Each of these values represents a different type of monomer used to build the molecule (among biological macromolecules $`\kappa =4`$ for nucleic acids G,A,C,T; $`\kappa =2`$ purines and pirimidines; $`\kappa =20`$ for proteins). It is reasonable, therefore, to ignore the genotype-phenotype distinction; since the relevant Darwinian entities are replicating macromolecules, genotype and phenotype are aspects of one and the same object, its molecular sequence (see  for a comprehensive overview). In this way, for one particular string (or macromolecule, or genome) of type $`i`$, the various events that could happen inside the flux reactor can be visualized as the following single chemical reaction steps,

$$(T)+I_i\stackrel{W_{ii}}{}2I_i$$
(1)
$$(T)+I_i\stackrel{W_{ji}}{}I_i+I_jji$$
(2)
$$I_i\stackrel{D_i}{}(M)$$
(3)
$$I_i\stackrel{\mathrm{\Phi }_0}{}0.$$
(4)
In these reactions, it is assumed that the available amount of $`M`$ and $`T`$ monomers is constant. The reaction (1) denotes the self-replication, or error-free replication, of molecule $`I_i`$, and the reaction (2) takes on the erroneous replication, or mutation, of $`I_i`$ that can lead to a new molecule $`I_j`$. The replication matrix $`W`$ takes into account the primary structure of the macromolecules. This feature is what distinguishes this model, a sequence space model, from models of population genetics . More specifically, its elements are given by

$$W_{ii}=A_iq^L$$
(5)
and

$$W_{ij}=\frac{A_j}{\left(\kappa 1\right)^{d(i,j)}}q^{Ld(i,j)}\left(1q\right)^{d(i,j)}ij,$$
(6)
where $`A_i`$ is the replication rate (or fitness) of molecules of type $`i`$, which tells us how fast new $`I_i`$ are synthesized in the next generation, independently of whether the copies are correct or not. The parameter $`d(i,j)`$ is the Hamming distance between strings $`i`$ and $`j`$, i.e., the number of monomers or positions at which these two sequences differ. Here, $`q[0,1]`$ is the fidelity parameter of each monomer, i.e., the probability of inserting the correct monomer for any position. Correspondingly $`\mu =(1q)`$ is the error rate (or mutation rate) per monomer, which is assumed to be the same for all monomers. Thus, self-replication and mutation are represented by two autocatalitic reactions, homogeneous and heterogeneous, respectively. The chemical decomposition of the molecular specie $`I_i`$ is represented by reaction (3), so that $`D_i`$ is a general decay rate parameter. The reaction (4) represents the efflux of the molecule $`I_i`$ by diffusion, where $`\mathrm{\Phi }_0`$ is a global diffusion flux that is assumed to be the same for all molecules.

The relevant variables of this dynamical system are the concentrations of each macromolecule $`x_i=\left[I_i\right]`$. Then, inside the reactor, the concentrations $`x_i`$ of molecules of type $`i=1,2,\mathrm{},\kappa ^L`$ evolve in time according to the following differential equations

$$\frac{dx_i}{dt}=\underset{j}{}W_{ij}x_j\left[D_i+\mathrm{\Phi }_0\right]x_i.$$
(7)
The quasispecies model can thus be described as a population dynamics model given by the full set of phenomenological differential equations.

The equations (7), however, do not yet lead to competitive selection in our model system. Moreover, in order to induce selection pressure in the system described by (7) it is necessary to impose some type of restriction to it. To appreciate the effect of the competition in the model system (that may or may not lead to selection), let us imagine the ideal situation where the molecular species reproduce themselves without error, so that $`\mu =0`$ and $`W_{ij}=0,`$ $`\left(ij\right)`$. In this case, equations (7) become

$$\frac{dx_i}{dt}=\left[W_{ii}D_i\mathrm{\Phi }_0\right]x_i.$$
(8)
It is trivial to see that the solutions for the molecular concentrations are exponentials if $`\mathrm{\Phi }_0`$ is constant. Thus, if there are $`k`$ distinct species of molecules inside the reactor, the concentration of all species with replication rate $`W_{kk}>`$ $`D_k+\mathrm{\Phi }_0`$ grows exponentially as time goes by, while all those with $`W_{kk}<D_k+\mathrm{\Phi }_0`$ die out. In this scenario, there is no competition and thus no selection. This behaviour of segregation in the model is show schematically in Figure 1(a) for the case of five types of binary molecules ($`\kappa =2`$, $`s_k=0,1`$), all with the same frequencies at $`t=0`$.

In this case of replication without error, we can now subject the molecular population inside the reactor to some form of global constraint. This can be made by keeping the total population of molecular species constant, $`_ix_i=N`$. For this purpose the global dilution flux $`\mathrm{\Phi }_0`$ will be adjusted in time so as to keep pace with the increase in total molecular concentration, i.e., it will be determined by the condition $`_idx_i/dt=0`$. Therefore, the global flux of molecules inside the reactor must satisfy the condition

$$\mathrm{\Phi }_0=\frac{\underset{i}{}(W_{ii}D_i)x_i}{N}.$$
(9)
Now, with this equation for $`\mathrm{\Phi }_0`$, the set of differential equations (7) and (8) is non-linear. In Figure 1(b) the temporal behavior of the same molecular populations of example of Figure 1(a) is schematically shown when we constrain the system to have a constant population. In fact, as $`\mathrm{\Phi }_0`$ increases due to the constant population condition, a higher number of molecules is segregated, and only the one with the highest productivity will survive. In the stationary state of this competition, the winner is called the master sequence  ($`I_m`$) and the stationary state of this system is termed selection equilibrium. The molecular selection process is, therefore, an environmental effect, through a constraint imposed on molecular population, and always leads to an unambiguous selection decision.

### 2.2 Quasispecies and error threshold

Let us discuss briefly the more general case, described by equation (7). If we allow for erroneous replication ($`\mu >0`$), considerations similar to the ones above lead to the following condition on $`\mathrm{\Phi }_0`$

$$\mathrm{\Phi }_0=\frac{\underset{i}{}\underset{j}{}W_{ij}x_j\underset{i}{}D_ix_i}{N}.$$
(10)
The selection equilibrium solution of (7) can be found in terms of the eigenvalues and eigenvectors of the replication matrix $`W`$. If we define the vector $`𝐲=(y_1,y_2,\mathrm{}y_L)`$, where the components represent the relative concentrations (or frequencies) of each molecular species in the population, $`y_i=x_i/_jx_j`$, we can write

$$W^{}𝐲=\lambda 𝐲,$$
(11)
where the diagonal elements of $`W`$ were modified to $`W_{ii}^{}=W_{ii}D_i`$ (see the appendices of  for details on how to transform the non-linear system described by (7) and (10) in a linear one). Then, we define a quasispecies precisely as the dominant eigenvector $`𝐲_{\mathrm{max}}`$ associated to the largest eigenvalue $`\lambda _{\mathrm{max}}`$ of the replication matrix $`W^{}`$. This eigenvector describes the exact population structure of the population: each mutant $`I_i`$ is present in the quasispecies with frequency $`y_i`$ ($`_iy_i=1`$). Note that the largest eigenvalue is exactly the average replication rate of the quasispecies, $`\lambda _{\mathrm{max}}=_iA_iy_i`$. In this scenario, the frequency of a given mutant within the quasispecies does not depend on its replicative value alone, but also on the probability with which it is produced due to erroneous replication of others molecular species. Moreover, depending on the mutation rate $`\mu =(1q)`$, the master sequence and some (or all) mutants coexist as a quasispecies. Then, one of the main outcomes of this model is that, in equilibrium, selection does not, in general, lead to a homogeneous population formed by a single kind of molecular individual: a particular sequence is no longer the outcome of selection. A set of genetically distinct sequences, forming a mutant distribution centered around the master sequence, is produced instead.

The mutation rate in this model is the parameter that controls the width of this distribution, that is, how much the quasispecies spreads over the space of sequences. Moreover, one immediate consequence of the existence of a maximum eigenvalue is the appearance of a threshold relation. Then, for the binary version of this model, as the error rate $`\mu `$ increases, two distinct regimes are observed in the population composition: the quasispecies regime, characterized by the master string and its close neighbours, and the uniform regime, where the $`2^L`$ possible sequences appear in the same proportion. The transition between these regimes takes place at the error threshold $`\mu _t`$, whose value depends on the parameters $`L`$ and $`A_i`$ . This becomes a genuine phase transition from an adaptive to a disordered neutral phase in the limit $`L\mathrm{}`$ .

Hence, to proceed further we must specify the replication rate of each sequence, or, in a more concise form, the fitness or replication landscape (this term was originally coined by S. Wright ), a line drawn by all $`A_i`$ points related to the sequence space. In particular, the quasispecies concept and the error threshold phenomenon are illustrated more neatly by the so-called single-sharp-peak landscape. This simple, and probably the most studied of the landscapes, can be constructed by ascribing the replication rate $`a>1`$ to the master sequence $`(1,1,\mathrm{},1)`$. The remaining sequences, differing by as little as a single mutated monomer, are at a disadvantage $`a1`$ as expressed by their replication rate $`a^{}=1<a`$. In Figure 2 we present the steady-state molecular frequencies in the case where $`L=30`$ and $`a=10`$, as a function of their Hamming distances $`d=0,\mathrm{},L`$ from the master sequence. That is, the $`2^L`$ possible sequences were grouped into $`L+1`$ different classes of sequences, characterized solely by the number of mutant monomers ($`s_\alpha =0`$) they have, regardless of the particular positions they occupy inside the sequences. Then, if we take as a reference the case $`\mu =0`$, where only the master sequence survives in the selection equilibrium, the parts (a) to (c) of the figure show three situations in the quasispecies (or adaptive) regime, $`\mu <\mu _t`$: in (a) the master sequence is the more frequent, in (b) the master sequence is no longer the more frequent, and in (c), near the error threshold, the number of master copies becomes strongly reduced, but the quasispecies is evolutionarily most versatile because it produces a wide variety of mutants without destabilizing the master sequence. Part (d) illustrates the uniform (or stochastic) regime, $`\mu >\mu _t`$, where the $`2^L`$ possible genomes appear in the same proportion or the $`L+1`$ classes are distributed by the binomial distribution $`y_d=\frac{1}{2^L}\left(\genfrac{}{}{0pt}{}{L}{d}\right)`$.

Therefore, in this simple example of fitness landscape, without errors ($`\mu =0`$) or with too much errors ($`\mu >\mu _t`$), there is no evolution. In the first case there are no mutants; in the second case there is no adaptation. Moreover, it is possible to show that the existence of a threshold relation, as pointed out above, correlates the mutation rate with a maximum length of the sequences that can be reproducibly maintained by selection (see  for a simplified case, but with a very intuitive appeal). For the single-sharp-peak landscape:

$$L<\frac{lna}{1q}$$

That is, once the fidelity degree $`q`$ of the copying mechanism is fixed, the molecule size cannot exceed a given value $`L_{max}`$. On the other hand, it is impossible to increase the fidelity of the copying mechanism without changing $`L`$; this is the so-called Eigen’s paradox. It poses a serious difficulty in envisioning life as an emergent property of systems of competing self-replicating macromolecules.

Some final comments regarding the formulation described above are in order. It is valid only in the limit where the total number of molecules $`N`$ goes to infinity. Of course, all real populations are finite, and they will not behave in the deterministic way expected for an infinite population . As pointed out by Higgs , the steady-state of a finite population is a dynamic one in which the population can continue to evolve, and therefore it is not equivalent to an infinite population model. Finally, it should be added that the error threshold itself is not a general phenomenon, and can be absent even for landscapes as simple as the Fujiyama one  (for a lucid discussion on the limits and further developments of this model see Ref. and references therein). For all other landscapes, such as the non-stationary one described in the next section, more elaborate approaches are required to clarify the real meaning of the error threshold. As a side remark, and to establish another connection with what follows, the evolution of a finite population in a random fitness landscape has been recently shown to exhibit a punctuated scenario .

## 3 The Bak-Sneppen Model for Macroevolution

In his excellent book “How Nature Works”  Per Bak argues that: Complex behaviour in Nature reflects the tendency of large systems with many components to evolve into a critical state where minor disturbances may lead to events, called avalanches, of all sizes. This delicate state evolves without interference of any outside agent and so is a critical self-organized state, that appears as a consequence of the dynamical interactions among individual elements of the system. Self-organized criticality is now a widespread concept and many different systems are known to evolve according to this dynamics. Particularly, the sandpile model  is the best one to explain what the famous avalanches are and how to measure them. Due to lack of space, we are forced to go directly to the Bak-Sneppen model and to the avalanches that appear in the large-scale evolutionary process that can explain, for example, the mass explosion of the Cambrian period and the mass extinction of the Cretaceous, when the dinosaurs disappeared.

In this model there are $`I`$ species, each one occupying one site of an unidimensional lattice (a ring). Each species has a random fitness $`0f_i1`$. The simulation evolves according to the following rule:

* search for the smallest $`f_i`$ corresponding to species $`i`$;
* change $`f_i`$, $`f_{i1}`$ and $`f_{i+1}`$ for 3 other values randomly chosen (mutation or extinction);
* return.

At the start of the simulation the fitness on average grows, although there are fluctuations up and down. However, after a transient period, the system reaches a stationary critical state where the fitness does not grow any further on average: all species have fitness above some threshold very close to $`2/3`$ (state of “stasis”). Consider a point in time when all species are over the threshold; at the next step the least fit species (right at the threshold) will be selected, eventually starting an avalanche or “punctuation” of mutation events. (In fact, whenever the fitness of a given species changes, it is possible to think that the species has undergone a mutation or that it has become extinct.) After a while, the avalanche stops, when again all species have fitness above the threshold. The avalanche size corresponds to the number of steps needed to recover the state of stasis (or equivalently, to the number of active species between two consecutive states of stasis). When this process is repeated for large systems (large number of species), one obtains the number $`N(S)`$ of avalanches of size $`S`$ given by the power law:

$$N(S)S^\tau .$$

This distribution means that there is no characteristic size for the avalanches, as would happen if instead of a power law it had an exponential behaviour. The larger the system, the larger the possible maximum size of an avalanche is. Small avalanches, in which a few number of species become active, are much more frequent than large ones; however, the probability that a system-sized avalanche occurs, activating all species, is not zero. Such a dynamical behaviour can explain the extinction of the dinosaurs without using any external agent, such as meteorites colliding with Earth (but of course does not exclude such a possibility).

The important point here is the fact that a given species, though highly fit, can be chosen to mutate because one of its neighbours has a low level of fitness. In this sense the species live in a continuously evolving fitness landscape, differently from the quasispecies of the Eigen model in which the fitness landscape is fixed: once a replication rate is attributed to a given macromolecule, it never changes.

Finally, in the next section we introduce the Penna model for biological ageing. It is also a model for microevolution, where individual interactions are introduced through a competition for food and space as in the Eigen model, but one that presents an age-structure, without reproduction at early ages.

## 4 The bit-string Penna model

The most successful computational model for age-structured populations is by far the Penna model . One of the reasons for its success relies on a particularly well-suited computational representation of a genome by means of a sequence of bits, the bit-string. When grouped into computer words these strings can be efficiently operated on by very fast logical and bit-wise CPU instructions. Conceptually, the model is extremely simple, but the results it shows are far from trivial. It has also proved to be flexible enough to be of value in a number of different problems in population dynamics, and the recent literature is eloquent proof of this statement.

The biological support for the Penna model comes from the mutation accumulation theory of senescence . In essence, this theory relates the evolution of senescence to the action of age-specific deleterious mutant alleles and the maintenance of some of these genes by the combined effects of mutation and selection pressures. In response to this conflict, deleterious mutant alleles with late ages of action would have higher equilibrium frequencies of affected individuals than genes with similar effects on survival early in reproductive life . A sufficiently large number of loci capable of mutating to deleterious alleles with age-specific effects would therefore generate a net decline in survival with advancing adult age.

### 4.1 Haploid asexual population

We will begin the description of the Penna model with its simplest implementation, the one that addresses haploid asexual populations. For a much more detailed explanation, together with a sample computer code, the reader is directed to Ref. . The basic structure of the model is the age-structured genome pool. Each genome from this pool is associated to an individual, and its pattern of alleles summarizes the genetic heritage acquired by this particular individual. The only genes represented are those with age-specific effects. Each gene appears in a particular position (locus) of the bit-string, and this position is associated with the age from which it becomes effective. A gene can appear in two alleles, for the $`2`$ possible values of the bit at one locus: a bit set to $`1`$ represents the deleterious allele, and it is set to $`0`$ for the non-deleterious variety.

This basic and simple structure allows the computation of the number $`A`$ of active deleterious mutations at any point of an individual life time: one has only to add up the bits of the bit-string from the first locus to the one corresponding to the actual age of that individual. The size $`S`$ of the bit-strings is the first parameter of the model, and determines the maximum theoretical age of the individuals. In this context, age is not to be understood as measured in human years, but rather in a species-specific time unit.

Selection pressure is modeled by the introduction of a threshold $`T`$ for the number of deleterious mutations that can be simultaneously active in a living individual’s genome. Usually, the probability of death because of genetic causes is assumed to be given by a step function $`\mathrm{\Theta }(TA)`$, with $`\mathrm{\Theta }(x)=1`$ if $`x0`$, which introduces a high degree of non linearity in the model. Smoother functional forms have also been used , but the fundamental results of the Penna model do not appear to be sensitive to this choice.

Death for non-genetic causes, representing the outcome of intra-species competition for the limited resources of the environment, is modeled by a density- and time-dependent quantity, the Verhulst factor. This is a mean-field death probability, given by $`V(t)=N(t)/Popmax`$, where $`N(t)`$ is the total population at the beginning of time step $`t`$ and $`Popmax`$ is another parameter of the model that quantifies the above mentioned environmental constraints on the size of the population. The introduction of the Verhulst factor, or some equivalent form of limiting factor for the total population, is a necessity in simulational models to avoid population overflow, although its biological motivation and implementation strategy have recently raised some interesting questions .

In an asexual population all the individuals are female. After reaching some minimum age $`Minage`$, and until she is $`Maxage`$ time periods old, every female in the population gives birth to $`B`$ offspring at every time step. At this moment, the genetic heritage of the mother is copied to her offspring, and mutations can occur. This process is simulated by generating, for each of the $`B`$ offspring, a genome cloned from that of the mother. On this bit-string, $`M`$ mutations are introduced in randomly chosen loci. It is usual to consider only deleterious mutations, since they are the overwhelming majority in nature. A clever coding trick is in order at this point: for the $`\mathrm{\Theta }`$ function implementation of the rule for genetic deaths, one might as well compute at this moment the programmed age of death, by adding up the value of the bits from locus $`0`$ until this sum reaches the threshold $`T`$; the last locus to be computed in this sum corresponds to the age at which this individual will die, if not sooner because of the Verhulst dagger. This trick avoids having to compute, at each time step, the number of active deleterious mutations for each individual, a time-consuming and redundant operation.

With the elements above, a typical simulation of the Penna model undergoes the following steps, described in a simple auto-explanatory meta-language:

* An initial population of $`Inipop`$ individuals is generated. The usual choices are either a mutation-free population, in which all the bits of all the bit-strings are $`0`$, or one in which each individual has a random number of mutations, between $`0`$ and $`S`$, in randomly chosen loci.
* FOR each time step:
* COMPUTE the Verhulst factor $`V`$;
* FOR each individual:

INCREMENT age by $`1`$;

IF age is smaller than the programmed age of death AND a randomly tossed number in the interval $`(0,1)`$ is larger than $`V`$ THEN reproduce, giving birth to $`B`$ mutated clones;

ELSE she dies: her genome is erased from the genome pool.

These dynamic rules are executed for $`Nsteps`$ time steps. In the last $`Asteps`$ averages are taken over the population, and constitute the outcome of the simulation. The quantities usually computed are:

* The age distribution of the population, i.e., a histogram of the number of individuals at each age, with some normalization.
* The survival probability as a function of age, defined as the ratio $`N(a+1,t+1)/N(a,t)`$ between the population with age $`a+1`$ at time step $`t+1`$ and the population that at the previous time step had age $`a`$.
* The mortality rate, defined as the logarithmic measure of the population decay, normalized so that its value at age 1 is zero and where only deaths due to genetic causes are considered.
* The genetic state of the population, which can be measured for instance by the fraction of defective genes in the population at each locus.

### 4.2 Early results

Among the many results obtained with the Penna model already published, we chose two of the most spectacular to comment on. The first one shows the agreement of the simulation results for the mortality rate with the empirically derived Gompertz law . This law was proposed in the late $`19`$th century to account for the observed mortality rate of the German population, and has since then been verified by a number of observations of both human and mayfly populations. It states that the increase of mortality with age is exponential. The Penna model was the first computational model for ageing that could reproduce this result. In Figure 3 we show data derived from a simulation - the particular parameters used can be found in the caption.

The second result deals with the catastrophic senescence observed in semelparous species. These are species that reproduce only once in a lifetime, such as the pacific salmon. Ageing in these species is called catastrophic because the females die soon after giving birth for the first and only time. To simulate a semelparous population with the Penna model one only needs to set $`Minage=Maxage`$, as opposed to having $`Minage<Maxage`$, which is the normal (iteroparous) case. These simulations could show that the catastrophic senescence effect is due to the lack of selection value of age-specific genes that become effective after the reproduction age . With the absence of selection, the mutation pressure turns on deleterious alleles for all these genes. In the computational representation provided by the Penna model, all loci associated with ages greater than the reproduction age become set to $`1`$ and the individual dies because of the accumulation of deleterious mutations in the first time step following its progeny. The consequences of this pattern of fixation of deleterious alleles can be clearly seen in the plot for the survival rate shown in Figure 4 which compares semelparous and iteroparous populations with equal parameters. The survival rate for the semelparous population has an abrupt decay to zero at age $`=Minage+1`$, in sharp contrast with the iteroparous one. Setting $`Minage=Maxage`$ also greatly simplifies the analytical formulation of the Penna model and allowed a formal derivation of the catastrophic senescence effect .

The analysis of Figure 4 -b suggests an interesting puzzle. The maximum age of an individual in this population is its last age of reproduction: there is no post-reproductive life!. The conflict between selection and mutations has a very simple outcome, and the individuals die as soon as they loose their function of perpetuating the species. This is not what is seen in nature, however. In human and some other mammal populations, females live after having ended their reproductive period (after menopause). To address this question the Penna model needs an extension to sexual diploid populations.

### 4.3 Diploid populations

These are species where the genetic information is carried by two homologous strains. The effectiveness of a deleterious allele at one locus now depends on a combination of information carried by the two strings. The concept of dominance appears in connection with this issue. If a dominant allele appears in one locus of any of the two strings, it is effective, irrespective of the allele that is present in the homologous locus of the other. For a non-dominant allele to be effective, on the other hand, it must be present in both homologous loci.

The appearance of diploid organisms in the life story of our planet marks also the onset of recombination and sex in the process of reproduction , although there are quite a few examples in nature of asexual reproducing diploid species where recombination is also present. This last strategy of reproduction, asexual but with recombination, is called meiotic parthenogenesis. Reproduction in diploid organisms involves the generation of a haploid cell, a process called meiosis, containing a subset of the genetic material of the parent that is to be transmitted to the offspring. To generate this haploid cell, alleles of the two homologous strains are reshuffled and recombined to form a single strain in a process called recombination. Diploidity and sex raise a whole new set of questions in population dynamics, some of which have been already studied in the framework of the Penna model and will be discussed in the sequel.

The extension of the Penna model to deal with diploid populations is rather straightforward . Now, the genome of each individual is composed of two bit-strings to be read in parallel: two alleles, one from each string, have to be taken into account to decide if the character of each age-dependent gene is deleterious or not. If this is an homozygote locus, i.e. if both alleles have the same character, then the gene has also this character. If the locus is heterozygote, the decision has to respect the dominance rule. To implement this rule, an extra $`S`$ bits long bit-string is generated at the beginning of a simulation indicating, for each locus, which is the dominant character. In $`D`$ of the $`S`$ loci, randomly chosen from a uniform distribution, the deleterious character is dominant. For these loci, a $`1`$ bit set in any of the two strains suffices to define the character of this gene as deleterious. For the $`SD`$ remaining loci, the deleterious allele is recessive, and the homologous bits of the two strings have to be both set to $`1`$ for any of these loci to have a deleterious effect.

Meiosis in the Penna model can be easily tailored to mimic nature. As a side example, apomyctic parthenogenesis, a diploid asexual mode of reproduction without recombination, would be represented by a simple random choice of one of the two genetic strains. For recombination, the usual is to select a random position out of the $`S`$ loci and cut the two strings at this position. Two new strings are generated by crossing the resulting four pieces: the left side coming from one of the strings is attached to the right side coming from the other. Of the two new strings, one is randomly chosen, and constitutes the genetic material to be inherited by the newborn. For meiotic parthenogenesis, this single strain is cloned: before mutations, the new genome is totally homozygote. For sexual species, a male individual is randomly selected and his genetic material also undergoes meiosis and recombination. The two strains, one coming from the female and another from the male, form the genome of the newborn. Its gender is now randomly chosen, with equal probabilities. On each of these two strings, $`M`$ (deleterious) mutations are added at randomly chosen loci. In Figure 5 we illustrate the procedures for reproduction of haploid asexual, diploid sexual and diploid meiotic parthenogenetic populations with simple examples.

Now we are in a position that allows us to address the puzzle with which we ended the last subsection, namely, why do women live as long as men instead of dying immediately after their reproductive period? And we show in Figure 6 the resulting survival rates when sexual reproduction is added to the model. Semelparous populations still suffer catastrophic senescence, even when the males do not loose their reproductive capacity at any age . But for iteroparous species, females now have a post-reproductive life and there is no difference between male and female survival rates. The presence of males in the population bring some benefits to the females, after all! But this observation does not ends our quest, for in the same Figure we also see that populations for which females do not undergo menopause, and are able to breed for their entire life, have a larger life span. Why would then nature “invent” menopause?

The answer to this question is now rather elaborate, and has been suggested already in the biological literature. It comes as a result of three new components, and the ability to acomodate and manipulate them is a demonstration of the flexibility of the model. These new components are the need for parental care of the newborn during a certain period of time, the reproductive risk that increases the death probability for a female at the moment of delivery, and the transformation of one of the parameters of the model, the maximum age of reproduction $`Maxage`$, in an individual and genetically acquired characteristic, now subject to mutations. The first new component is implemented in the model by requiring an infant to have a living mother in order to survive during its first $`Apc`$ periods of life. For the reproductive risk, a dependence on the number of active deleterious mutations $`A`$ was introduced. Thus, a female can die, with a probability essentially given by $`A/T`$, at the moment of delivery. If the newborn is a female, she inherits the same $`Maxage`$ of her mother with some probability $`P`$, or mutates to have it increased (decreased) by $`1`$ with probability $`(1P)/2`$. The Darwinian dynamics of the model is sufficient to produce a self-organization of the distribution of the menopause age of the female population, showing that inhibition of reproduction after a certain age is actually beneficial for the species as a whole . In Figure 7 the distribution of the age of menopause onset throughout the population is shown in a comparison between a population where parental care is not needed and one where an infant aged below a minimum would die if the mother was not any more present.

### 4.4 The evolution of sex

The maintenance of sexual reproduction among the great majority of species in nature, in spite of the inherent cost of having to produce males to assure reproduction, is still one of the great puzzles of biology. In fact, a simple reasoning shows that the need of two parents to generate even a single offspring should give asexual varieties a two-fold advantage over sexual ones . The advantage of sex relies on its ability to create greater genetic diversity, since the pool of alleles from which the newborn genome is extracted is different for each mating pair. Asexual reproduction, on the other hand, generates new genomes from a more limited pool, since there is only one parent involved. It is not clear though in what circumstances this greater diversity provided by sex would give it the upper hand against, say, meiotic parthenogenesis. One could argue that the species would benefit from the cloning of well-fitted genome, only possible in asexual reproduction, and that sex would make these genomes short-lived for exactly the same reasons it can create diversity. Sex could possibly be more efficient in getting rid of bad mutations, by bundling them together through mating and expelling them from the genetic pool through selection. These are as yet open questions, and the Penna model has been used to address them, even when the age structure of the population is not the issue. This is possible because of its particularly simple implementation of a selection mechanism, which is a general requirement for any model to be useful in this context.

For the evaluation of diversity in the Penna model one measures the number of different alleles, or bits, for each pair of genomes in the population. The resulting distribution of this so-called Hamming distance has a Gaussian character for any reproduction strategy. The comparison between sex and meiotic parthenogenesis show a similar width for the distributions, but a higher value for the distance where it peaks in the sexual case . A particularly interesting reflex of this greater diversity can be seen when a genetic catastrophe is provoked at a chosen time step by instantaneously setting an extra deleterious mutation at a chosen locus in all the individuals. An equivalent, and perhaps more easily translated into biological terms, procedure would be to decrease the threshold of deleterious mutations $`T`$ by one, representing a sudden depletion of vital environmental life protection resources. Thanks to their greater diversity, sexual populations are always resistant to this catastrophe, whereas simple asexual reproduction leads invariably to extinction, as shown in Figure 8. For meiotic parthenogenesis, the outcome is not so clearly cut, but extinction is a possibility .

A number of theories have been put forth to try to explain the evolution and maintenance of sexual reproduction. In the center of this debate is the so-called “Red Queen” hypothesis, that relies heavily on the ideas of diversity discussed above. In essence, it holds the action of genetically matching parasites as responsible for creating a rapidly changing environment. In this unstable ecology, only varieties that can mutate their genomic pool at least as fast as the adaptation of the parasites proceed can survive. The theory derives its name from this endless race, quoting from the Red Queen of Lewis Carol’s Alice in Wonderland: “It takes all the running you can do, to keep in the same place.” In fact, recent observations of competing varieties of a freshwater snail, Potamopyrgus antipodarum, have shown that there is a strong correlation between the prevalence of one reproduction regime and the concentration in its habitat of the trematode Microphallus, a parasite that renders the snail sterile by eating its gonads . Namely, the asexual variety is predominant where the parasite appears in small concentrations, whereas higher concentrations of the trematode forces the species to prefer a sexual regime.

This correlation could be shown to exist in simulations of a conveniently modified Penna model. The parasites are represented by a dynamically changing memory bank of genomes of some fixed number of entries. Each entry is modified if it comes into contact with the same genome twice in a row; in this case, it memorizes this pattern and stores it in the memory bank. At each time step, before the reproduction cycle, each female of the population is probed by a fixed number $`E`$ of randomly chosen entries of the parasite bank. If one of these entries is a perfect match for the female’s genome, she is rended sterile and can no longer reproduce. The number of parasite exposures $`E`$ is an indirect measure of the parasite concentration in the habitat. For the host population, the reproductive regime of the females is no longer a fixed character, but can mutate with some small probability. The simulations begin in the absence of the parasite infestation, and the initial population is set to have a sexual reproductive regime. As soon as the meiotic population appears, due to mutations in the reproductive regime, it overrides the sexual variety, in a demonstration of the two-fold disadvantage of sex above mentioned, and sex barely subsists due to infrequent back-mutations from the asexual variety. At some time step, the parasite infestation is turned on. The resulting predominant variety is going to depend solely on the intensity of this infestation, as measured by the exposure parameter $`E`$. For small values of $`E`$, the asexual variety has the upper hand. As $`E`$ is increased, a first-order transition is seen to a configuration dominated by the sexual population . Figure 9 shows the fraction of females in the population that reproduces sexually, as a function of the exposure parameter $`E`$. The sudden jump in this fraction signals the order of the transition.

## 5 Other applications of the Penna model

It is not difficult to find in the animal kingdom species that live and work in sexual pairs, but sometimes have an extra-pair relation, like the Scandinavian great reed warbler, chimpanzees, etc… (of course some disgusting men also belong to this category). As already shown by Martins and Penna , such a behaviour increases the genetic diversity and may lead to better fitted offspring depending on how females select the males for an extra-pair relation. However, as in the case of the parasites presented before, nature not always chooses the strategy leading to the highest reproduction rate. Another example is the California mouse , one of the rare monogamous species that have been found. These mouses live in a extremely cold place, and in order to survive the pups must be continuously heated by the body of one of the parents. When the male abandons the nest, the female very often kills the babies. Starting with a population with half of the males faithful and the other half non-faithful, Sousa and Moss de Oliveira  have shown that depending on the death probability of the abandoned offspring, the population may self-organize in a situation where all males are faithful despite of reproducing lesser. In their simulations, monogamy is paternally transmitted and exclusively related to parental care. In fact it is already known that there some genes responsible for maternal care that are paternally transmitted .

Another interesting result obtained with the Penna model concerns the higher mortality of males when compared to the females one. The mortality curve as a function of age for females is lower since birth until advancing ages (around 90 years), when both males and females mortalities become equal. Stauffer et al  introduced somatic mutations (that are not transmitted to the offspring) into the model, atributting to males a higher somatic mutation rate than to females. With this strategy they were able to reproduce the observed behaviour of the mortalities, with the females one lower than that of males until advancing ages, when the genetic mutations dominate and the two curves collapse. Penna and Wolf  obtained the same result atributting to the females a higher value of the limit number of genetic diseases $`T`$ than to the males. However, the best strategy was proposed by Cebrat , and confirmed by Schneider  et al., that modified the model in order to distinguish the double X chromosomes of the females from the single one of males. Considering the mutations in the single X male chromosome as dominant mutations, even more realistic results were obtained.

Finally, we must say that it has also been possible to predict some unknown effects using this model. For instance, it has been shown that in small populations, if a given percentage of males is periodically substituted by the same percentage of males coming from a large population, the extinction of small populations due to inbreeding may be avoided . Also a nice strategy for fishing some species in which fertility varies with size in order to warrant a larger stock without loosing money, has been recently proposed by Racco and Penna .

## 6 Conclusions

We have shown that different evolutionary models found in the literature are connected to the biological ageing ones through the Darwin’s theory of selection of the fittest. The models we have presented here are the most simple ones, but their results are far from trivial. The Eigen model for microevolution makes clear the difference between segregation and selection, as well as the connection between mutations and diversity. The Bak-Sneppen model exemplifies the importance of a continuously evolving fitness landscape to simulate large-scale or macroevolution.

Finally we have presented the Penna model for biological ageing and some of its most important results. Ageing is an unavoidable process (except for some rare individuals like D. Stauffer, who is eternally young) and has been extensively studied by many different scientists, since a very long time. Although the evolutionary theories for senescence have appeared around 1950, Monte Carlo Simulations on this subject started only after the publication of the Partridge-Barton analytical mathematical model in 1993. The Penna model is now the most widespread Monte Carlo technique to simulate and study the different aspects of population dynamics, including ageing. In this review we have focused attention on results concerning the differences between reproductive regimes and the advantages of sexual reproduction (although our arguments have not proved good enough to convince the rare individual just mentioned above that males are all alike but still useful).

## Acknowledgments

S.M.O. and J.S.S.M. thank P.M.C. de Oliveira, D. Stauffer, T.J.P. Penna, A.T. Bernardes, A.O. Sousa and S. Cebrat for many discussions and collaborations in previous papers. S.M.O.’s work is partially supported by the Brazilian Agencies CNPq, Faperj and Capes. D.A.’s work is supported by FAPESP. J.S.S.M.’s work is supported by DOE grant DE-FG03-95ER14499.


format: text

text:
# Ground state properties of the 2D disordered Hubbard model

## I Introduction

The electronic transport properties of disordered systems have been the subject of much investigation in physics. From the pioneering work of Anderson (1958) it is known that in three dimensions (3D) the eigenstates of a non-interacting electron gas in a random potential become localised at the Fermi energy above a critical value of the disorder strength $`W_c`$. In this regime, the eigenstates decay exponentially in space and hence cannot carry a current; thus the system is an insulator. For disorder strength, $`W`$, lower than $`W_c`$, the eigenstates are extended and diffusive transport takes place in the system in accordance with Ohm’s law. However, for two-dimensional (2D) systems, it was shown by the scaling theory of Abrahams et al that all states are localised for any disorder strength. Thus, it appears that there is no metal-insulator transition (MIT) for non-interacting electrons in 2D. The properties of non-interacting electrons in random potentials have since been studied systematically and the main physical effects have been understood. In this context, the experimental observation by Kravchenko et al of a transition from insulating to metallic behaviour, as seen in the resistivity as a function of temperature of the 2D electron gas, came as a great surprise to the community. The existence of this transition from insulating to metallic behaviour as a function of density has been confirmed by other groups. The experiments were carried out on very high mobility, low electron density ($`n_s`$) samples which correspond to a regime where the electron-electron interactions ($`E_{ee}`$) are much stronger than the Fermi energy ($`E_F`$), such that the dimensionless parameter $`r_s(E_{ee}/E_F)`$ lies in the range 5-50. This indicates the importance of electron-electron interactions in these systems. Further, it was shown experimentally that the application of an in-plane magnetic field ($`B_p`$) drives the system insulating. Since such a field can only couple to the spin, this experiment indicates the important role played by the spin degrees of freedom. While the early experiments have stimulated a spate of new experimental results, there has been no satisfactory theoretical explanation of the phenomenon of the 2D MIT, to date.

The effects of interactions in disordered systems have been studied from the metallic side in great detail, where the interactions are relatively weak. The study of interactions in the localised phase were mainly carried out within the mean-field approximation, which led to a number of important results, for example, the Efros-Shklovskii gap in the density of states near the Fermi level, but could not take into account quantum interference effects, important in this many-body system. The investigation of a simple model of two interacting particles (TIP) in the localised phase showed that short-range attractive/repulsive interactions can lead to destruction of localisation and propagation of pairs of particles on a length scale much larger than the one-particle localisation length. Thus, the effects of interaction on the localised phase are non-trivial and deserve a detailed study. This, however, is not an easy task. Indeed, even the analytical expressions for the matrix elements of the interaction in the localised phase are not known, hence, numerical studies of the problem become important.

Recent numerical approaches to the question have included the studies of persistent currents by exact diagonalisation of small 2D clusters and Hartree-Fock based calculations without and with residual interaction. These approaches led to some interesting indications but did not allow the study of sufficiently large systems and/or sufficiently many particles. Other approaches based on level spacing statistics of many-body states made possible the study of larger systems and showed the existence of ergodic (delocalised) states for low energy excitations, but not at the ground state. All these studies were carried out for spinless fermions. Recently, the properties of fermions with spin on a disordered 2D lattice were investigated using a finite temperature quantum Monte Carlo (QMC) method. The temperature dependence of the resistivity obtained numerically indicated a transition from insulating to metallic behaviour for sufficiently strong interactions and weak disorder strength. However, these calculations were carried out at finite temperatures and technical problems (”fermion sign problem”) did not permit the analysis of the ground state. This is not completely relevant to the experiments which were carried out at temperatures much below the Fermi energy and thus require a better understanding of the properties of the ground state, as in the general scenario of quantum phase transitions.

To investigate the properties of the ground state of the disordered interacting fermionic system, we choose the Hubbard model with site diagonal disorder. This could be considered as an important first step on the way to investigations of more complicated models with Coulomb interactions, which might be more appropriate for experiments at low densities. We study the ground state of our model on the square lattice by the projector QMC (PQMC) method. We use different characteristics to investigate the extent of the ground state wavefunction for a broad range of model parameters: disorder strength ($`W`$), interaction strength ($`U`$) and filling factor ($`\nu `$). The studies were carried out in the $`S_z=0`$ sector, with equal numbers of particles with up and down spins. This is thus the first numerical study of the ground state of a disordered, interacting system of fermions with spin.

This paper is organised as follows. In the next section, we describe the model and the method used. In the third section we present our results for the averaged Green function, the charge densities and the inverse participation ratios, all of which we use to characterise the ground state. We present a summary of our conclusions in the last section.

## II Model and method

The two-dimensional disordered Hubbard model on a square lattice is given by

$$\begin{array}{c}H=H_A+H_I\hfill \\ =\left(t\underset{ij,\sigma }{}\widehat{a}_{i,\sigma }^{}\widehat{a}_{j,\sigma }+\underset{i,\sigma }{}ϵ_i\widehat{a}_{i,\sigma }^{}\widehat{a}_{i,\sigma }\right)+U\underset{i}{}\widehat{n}_i\widehat{n}_i\hfill \end{array}$$
(1)
where the $`\widehat{a}_{i,\sigma }^{}`$ ($`\widehat{a}_{i,\sigma }`$) are the creation (annihilation) operators for a fermion of spin $`\sigma `$ at site $`i`$ with periodic boundary conditions, $`\widehat{n}_{i\sigma }`$ is the number operator for spin $`\sigma `$ at site $`i`$, $`t`$ is the hopping parameter, the Hubbard parameter, $`U`$, measures the strength of the screened interaction and $`ϵ_i`$, the energy of site $`i`$ is a random number drawn from a uniform distribution $`[W/2,W/2]`$, which parametrizes the disorder. The first two terms represent the Anderson Hamiltonian and the last term represents the interaction $`H_I`$. In the limit $`W=0`$, this Hamiltonian reduces to the usual Hubbard model. The filling factor $`\nu =N_p/(2\times N^2)`$, where $`N_p`$ is the number of fermions (particles) and $`N`$ the linear dimension of the system; thus the total number of sites is $`N^2`$.

We obtain the ground state properties of this model by the PQMC method. The PQMC method was initially developed for the Hubbard model and has been used to obtain reliable results for large lattices. The method can be generalised in a straightforward manner to include random site energies. We now present some details of the calculation for completeness and refer the reader to the literature for more detailed accounts.

The PQMC method consists in obtaining the true ground state $`|\psi _0`$ of the Hamiltonian(1) by projection from a trial wavefunction $`|\varphi `$ that is not orthogonal to the true ground state of the system,

$$|\psi _0=\underset{\mathrm{\Theta }\mathrm{}}{lim}\frac{e^{\mathrm{\Theta }\widehat{H}}|\varphi }{\sqrt{\varphi |e^{2\mathrm{\Theta }\widehat{H}}|\varphi }}.$$
(2)
The trial wavefunction is usually formed from the eigenstates of the non-interacting Hamiltonian (orbitals filled up to the Fermi level). In this case, we choose the eigenstates of the Hamiltonian $`H_A`$, thus including the random potential. To carry out Monte Carlo (MC) simulations of this quantum Hamiltonian, it is first necessary to map it onto an effective classical Hamiltonian. Thus, the projection operator $`\mathrm{exp}(\mathrm{\Theta }\widehat{H})`$ is first Trotter decomposed as $`\left(\mathrm{exp}(\mathrm{\Delta }\tau \widehat{H}_A)\mathrm{exp}(\mathrm{\Delta }\tau \widehat{H}_I)\right)^L`$ with $`\mathrm{\Theta }=\mathrm{\Delta }\tau \times L`$. This introduces a systematic error of order $`(\mathrm{\Delta }\tau )^2`$ due to non-commutation of $`\widehat{H}_A`$ and $`\widehat{H}_I`$. The interaction is then decoupled by a discrete Hubbard-Stratonovich (H-S) transformation, by the introduction of $`N^2\times L`$ Ising-like fields. Since the complete summation over these degrees of freedom is too time-consuming to be practical, the method reduces to a MC sampling of physical properties, which is the second source of error, the statistical error. It is important to note that during the MC process, each configuration of Ising spins is assigned a weight, which is interpreted as a probability. This quantity is positive definite only at half-filling for the uniform Hubbard model. The problems that arise from the non-positive-definite nature of this quantity are referred to as the ’fermion sign problem’ in the literature and are known to be particularly severe slightly away from half-filling in finite temperature methods and restrict the lowest temperature that can be attained in the simulation (in the clean limit).

We have studied system sizes of up to $`10\times 10`$ at quarter and one-eighth fillings (50 and 25 particles). We carried out extensive checks on the MC parameters to assure ourselves of convergence, as described in Ref., but in the presence of disorder. Thus, we chose $`\mathrm{\Theta }=3.0`$, with $`L=30`$, to have $`\mathrm{\Delta }\tau =0.1`$. This, with the symmetric Trotter decomposition reduces the systematic error to $`(\mathrm{\Delta }\tau )^30.001`$. We checked for statistical convergence of our data in several ways. By varying the number of sweeps after equilibration, we determined that 1000 MC sweeps are sufficient for equilibration and 2000 further sweeps for property estimates. We carried out measurements until the standard deviations on our values were of the order of the systematic error. We also tested our results against results obtained from exact diagonalisations for small system sizes and the results for the charge density, $`n_i=n_i+n_i=_\sigma \widehat{n}_{i\sigma }`$, presented in Fig. 1 show good agreement with the exact results. Convergence is of course the best for the ground state energy as compared to other physical quantities, and we have a relative accuracy of $`10^3`$ when compared to exact calculations. As for the effect of disorder on the sign problem, it was possible to study the $`10\times 10`$ lattice for $`U/t=6`$ for disorder strengths $`W/t`$ of up to 7-10. Our measure of the severity of the sign problem is to consider the quantity f = 1- (number of negative determinants)/ (total number of determinants). In all the cases considered, we have f = 0.999, which indicates that the sign problem is under control.

From the simulations, it is possible to obtain ground state expectation values of the single particle Green function, $`G_{ij}=_\sigma \widehat{a}_{i,\sigma }^{}\widehat{a}_{j,\sigma }`$, where the average is a MC average. Further, we can obtain ground state expectation values of other one- and two-body operators, such as the charge density and the charge-charge correlation functions. Each disorder realisation constitutes a full PQMC calculation. The properties are averaged over 16 disorder realisations. Thus, we have obtained the evolution of the Green function with distance, the charge densities and the inverse participation ratios. Our results will be described in the following section.

## III Results and discussion

To characterise the properties of the ground state, we study the correlation function defined as

$$C(𝐫)=\frac{1}{N^2}\underset{i,j}{}|\widehat{a}_{i,}^{}\widehat{a}_{j,}+\widehat{a}_{i,}^{}\widehat{a}_{j,}|^2\delta _{ij,𝐫},$$
(3)
where $`𝐫=𝐢𝐣`$ is the vector in the plane between the sites labelled $`i`$ and $`j`$, and the averages are carried out over the ground state eigenfunction and the different disorder realisations for all possible initial positions of $`𝐫`$, i.e. all corresponding $`i`$ and $`j`$. $`C(𝐫)`$ is simply related to the Green function $`G_{ij}`$ already defined in Section II. With this definition, $`C(0)=N^2_i(n_i+n_i)^24\nu ^2`$ in the limit of weak disorder and $`C(0)4\nu `$ in the strongly localised limit. The dependence of $`C(r)`$ on distance $`r`$ is related to the localisation of the eigenstate, i.e. we expect exponential decay of this quantity for localised states and slow decay at long distances for extended states. We also study the direction averaged correlation function $`\overline{C}(r)`$ which now depends only on the distance $`r=|𝐫|`$.

In Fig. 2 we show the decay of $`\overline{C}(r)`$ with $`r`$, for the 2D Anderson model, system (1) at $`U/t=0`$, for various disorder strengths. The change from flat behaviour with $`r`$ at weak disorder, when the eigenfunctions are delocalised in the finite sized system ($`W/t=2`$), to asymptotic exponential decay for stronger disorder, $`W/t10`$, (when the localisation length is smaller than the system size), is evident. We note that the initial non-exponential decay in the localised case is due to the fact that the ground state eigenfunction is a superposition of one-particle eigenstates of different energies. In fact, the one-particle localisation length $`l_1`$ of a state depends on its energy and therefore the many-body state, which is the Slater determinant of the 1-particle states up to the Fermi level, initially decays more rapidly. This is due to the low-lying states of smaller localisation lengths and it is only in the asymptotic limit that the decay of the many-body ground-state is determined by the maximum $`l_1`$ at $`E_F`$. This physical structure of many-body states complicates the observation of asymptotic exponential decay corresponding to the one-particle localisation length $`l_1(E_F)`$ at the Fermi level. Despite these complications, the asymptotic slope is seen to depend strongly on $`W`$ (Fig. 2), which is consistent with the exponential growth of $`l_1`$ with decreasing $`W`$ in 2D. In view of this, the investigation of the correlation function $`\overline{C}(r)`$ in the presence of interactions should tell us the impact of interactions on the localisation properties of eigenstates.

The dependence of $`C(𝐫)`$ on $`𝐫`$, shown in Figs. 3a,b for the 2D disordered, non-interacting Hubbard model model ($`U/t=0`$), also clearly shows localisation of the ground state eigenfunction. The decay (as seen from the contour plot Fig. 3b) is approximately symmetric in $`𝐫`$ which is due to averaging over different disorder realizations. Hence, it should be useful to study this quantity also in the interacting case, to clarify the ground state properties.

The behaviour of $`C(𝐫)`$ as a function of $`𝐫`$, for relatively strong interaction strength, ($`U=6t`$) is shown in Figs. 4a,b. The comparison with the non-interacting case (Figs. 3a,b) clearly shows that even such a strong interaction produces only a slight change in the correlation function. We observe similar behaviour for other disorder and interaction strengths ( $`5W/t10`$, $`0<U/t6`$, $`\nu `$ = 1/4, 1/8, results not presented here). This, in our opinion, provides direct evidence that even in the presence of strong interactions, the ground state of the system remains localised. This conclusion is futher supported by the data for the direction averaged correlation function, $`\overline{C}(r)`$ presented in Fig. 5. In fact, this direction average further smoothens fluctuations due to disorder. Indeed, even the introduction of relatively strong interactions ($`U/t`$ = 6) affects this function very weakly.

As an alternative test for localisation, we use another, more indirect method, similar to the approach presented in Ref.. As discussed in , we vary the amplitude of on-site disorder $`ϵ_i`$ for sites $`i`$ along one vertical line of the square lattice, as $`ϵ_ia\times ϵ_i`$ with $`a`$ = 1.1 and 1.3, corresponding to 10% and 30% change in disorder. We then study the charge density difference $`\delta \rho _x`$ produced by this perturbation, as a function of distance $`x`$ from the original line. We average over all sites with the same $`x`$ and additionally average $`\mathrm{log}|\delta \rho _x|`$ over 16 disorder realisations. The comparison of data for 10% and 30% variation shows that we are in the linear response regime. The results are presented in Fig. 6. For $`U/t=0`$, the response function shows a sharp drop from the initial peak followed by a slower decay at longer distances. This behaviour is qualitatively similar to the decay of the correlation function (Fig. 2, Fig. 5), for the same physical reasons as analysed above. Introducing interactions doesn’t at all affect the main structure of the curve which drops very quickly from the centre by more than one order of magnitude. We interpret this as a sign of a localised ground state. At the same time, we note that there is a slight difference introduced by interactions at the tails of the response functions. However this corresponds to a density variation less than 0.1%, which is at the limit of the accuracy of our calculation. In the light of the ensemble of data, we conclude that the ground state in the presence of interactions remains localised.

In a sense the Hubbard repulsion leads to local rearrangements of charge and does not seem to influence the long-distance properties of the system. This is clearly illustrated in Figs. 7a,b where the interaction leads to a more homogeneous charge distribution (note the change in the vertical scale) but doesn’t drastically change the global profile.

This point of view is further borne out by the data for the inverse participation ratio (IPR), $`\xi `$, presented in Fig. 8. The IPR, defined as $`\xi =(_in_i)^2/(N_p_i(n_i)^2)`$, gives the average number of sites visited per particle. $`\xi `$ is bounded from above at fixed filling, with $`\xi \xi _{max}=(2\nu )^1`$, corresponding to the weak disorder limit and from below with $`\xi \xi _{min}=0.5`$ in the strongly localised limit. We note that $`C(0)\times \xi =2\nu `$. We have studied the IPR as a function of system size and interaction strength at 1/4 (Fig. 8) and 1/8 (Fig. 9) filling. At $`U/t=0`$, the IPR naturally increases with decreasing disorder as states become more extended. At strong disorder ($`W/t5`$) $`\xi `$ is not sensitive to the system size since the states are localised and at fixed filling the IPR is counted per particle. The Hubbard interaction smoothly increases the IPR but does not introduce a qualitative change. This corresponds to the local reorganisation of charge introduced by $`U`$ leading to a more homogeneous charge density distribution, as discussed above. The data presented for 1/8 filling in Fig. 9 show qualitatively similar behaviour. However, the size variation is more restricted in this case ( $`N10`$ in our studies).

At this point, it is interesting to compare our studies with a recent finite temperature QMC study of a similar model. These authors considered the Hubbard model on a square lattice, with off-diagonal disorder, in contrast to our study. This is due to the fact that the method used, a finite temperature QMC method, suffers from a severe sign problem in the presence of diagonal disorder. For off-diagonal disorder, the situation becomes better, but the problem persists, restricting the lowest accessible temperatures. Their studies of several physical characteristics, including the conductivity (obtained by approximate analytic continuation of the imaginary time Green function), indicate the presence of an interaction induced metal-insulator transition in their model. However, this method is not adapted to analysis of the ground state properties. Our results are not in direct contradiction with this study, since it is fully possible that the ground state remains localised, while the low-lying excited states become delocalised. Such a situation has been observed in numerical studies of spinless fermions with Coulomb interactions on a 2D lattice with disorder. Our result directly demonstrates the localised nature of the ground state even in the presence of strong interactions. This is important in the framework of general studies of zero-temperature quantum phase transitions.

## IV Conclusions

We have used the projector quantum Monte Carlo (PQMC) method to study the ground state of the 2D disordered Hubbard model. This method allows us to study systems of up to 50 spin-1/2 fermions on a $`10\times 10`$ lattice, for interaction strengths $`U/t`$ up to 6 and a broad interval of disorder strengths. The comparison of several properties in the absence and presence of the Hubbard interaction allows us to conclude that interactions lead to local rearrangements of charge but do not destroy the localised structure of the ground state, within the range of parameter values studied.

These results indicate that short-range interactions are probably insufficient to bring about a quantum phase transition in the ground state of this system. Thus, it becomes important to consider the effect of long-range Coulomb interactions for electrons on a disordered lattice.

We thank M-B. Lepetit for useful discussions and the IDRIS, Orsay, for allocation of resources on their supercomputers.


text:
# Old White Dwarfs as a Microlensing Population

## 1. Introduction

Other authors in this volume cover the microlensing motivations much better than I can, so I shall suffice to remind you that one possible explanation of the microlensing events towards the LMC invokes a population of objects in the range $`0.30.8`$M. Potentially these could be either normal hydrogen-burning stars or white dwarfs, the burnt-out remnants of stellar evolution. To distinguish these populations, we turn to direct searches at optical wavelengths, since the latter population is $`10^310^4`$ as bright as the former. The number counts of faint red stars suggest that hydrogen burning stars cannot account for the microlensing population (Bahcall et al 1994; Graff & Freese 1996). The question I wish to address is how well one can constrain the white dwarf hypothesis by similar means.

## 2. White Dwarf Cooling

To derive a constraint from direct optical searches, we need to know how to recognise white dwarfs. White dwarfs of moderate age ($`<10`$ Gyr) emit approximately as black bodies. However, once the white dwarf cools to effective temperatures below $`5000`$ K, the atmospheric hydrogen resides in molecular form. This has dramatic consequences for the appearance of the white dwarf (Hansen 1998, 1999a; Saumon & Jacobsen 1999), because the dominant opacity source in such an atmosphere is collisionally induced absorption by H2, which absorbs primarily in the near infra-red. For black bodies of effective temperatures $`30005000`$ K the peak of the spectrum lies in the same wavelength region, so that the increased absorption leads to dramatic deviations from the traditional assumption of black body colours. The general trend of the colour evolution is that the flux is forced to emerge preferentially blueward of $`0.8\mu `$m. With the evolution of the black body flux redward, the peak of the spectrum for old white dwarfs lies around $`0.6\mu `$m.

A correct quantitative analysis of the atmospheric conditions is not only important for determining the colours of the white dwarf, but is also of critical importance for the cooling evolution. For white dwarfs cooler than $`T_{\mathrm{eff}}6000`$ K, the internal structure has no radiative zones. The bulk of the white dwarf mass resides in a degenerate core which is approximately isothermal due to the efficient thermal conduction by the electrons near the surface of the Fermi sea. This isothermal core is joined to a thin convective envelope which extends all the way to the photosphere. Thus, the entire white dwarf structure is critically dependant on the surface boundary condition and hence the atmospheric conditions. It must be noted that, while sophisticated analyses of white dwarf atmospheres have been around for several years (e.g. Bergeron, Saumon & Wesemael 1995), the white dwarf cooling calculations have lagged significantly behind, employing grey atmosphere boundary conditions. Thus, while one may trust the masses and gravities inferred in recent atmospheric analyses, one cannot completely trust the ages inferred for the older white dwarfs in most recent papers. To my knowledge, the only cooling calculations using boundary conditions based on proper atmosphere models are those of Hansen (1999a) and Salaris et al (2000). For isochrones and luminosity functions using the Hansen (1999a) models, see Richer et al (2000).

The models described above concern white dwarfs with pure hydrogen atmospheres. Empirically, a significant fraction of known white dwarfs appear to have no hydrogen in their atmospheres (Liebert et al 1979; Bergeron, Ruiz & Leggett 1997). When pure helium atmospheres reach effective temperatures $`<5000`$ K, they have no molecular component to provide significant opacity and, as a result, the photosphere lies at much higher densities, where the beginnings of pressure ionization start to provide a small ionized component. This higher density and the concomitant difference in the boundary condition results in much faster cooling for old white dwarfs with helium atmospheres. Furthermore, these objects are expected to appear approximately as black bodies (contrary to the hydrogen atmosphere case) but will be hard to detect because they cool much faster.

The matter of a component of the white dwarf population sporting helium atmospheres is one that is often not mentioned in the recent studies attempting to tie direct searches to microlensing populations. It is important to note that this represents an essentially unobservable (due to their much more rapid cooling) component and, as such, will always introduce some uncertainty into the comparison.

## 3. Observational Progress

### 3.1. The story thus far

The last couple of years have been fun for those interested in the question of halo white dwarfs. The realisation that old white dwarfs would appear somewhat bluer than black bodies meant that previous constraints on the local number density based on Hubble Deep Field point source counts were somewhat overstated. Much excitement was generated by the detection of proper motions of several faint blue point sources in the HDF North (Ibata et al 1999). The colours and magnitudes of these objects are consistent with old white dwarfs $`12`$ Gyr old at distances $`1`$ kpc. Given the radical consequences of this observation, some kind of confirmation would be nice. A third epoch, to confirm the measurement, has been approved but was postponed due to the HST gyro problems in December 1999. Nevertheless, the veracity of the HDF proper motions should be addressed within the next year or so.

In the interim, several other interesting results have come to light. Most importantly, Hodgkin et al (2000) presented a near-infrared spectrum of a nearby cool white dwarf, demonstrating the reality of the flux-suppression due to molecular hydrogen absorption. This provides a welcome check that the theory behind the predicted colour change is at least qualitatively correct. The high proper motion of this object also suggests membership of the Galactic halo. Another object showing such a long wavelength depression is LHS 3250 (Harris et al 1999)

A further interesting indirect argument supporting the white dwarf hypothesis comes from Mendez & Minniti (2000), who assert that the number counts of faint blue point sources (the class amongst which Ibata et al discovered proper motions) in the Hubble Deep Field South is approximately twice that in the HDF North. This is consistent with the objects (whatever they are) being Galactic in origin, since the southern field looks in a direction closer to the Galactic centre and consequently the stellar density is expected to be higher.

A note of caution, however, is sounded by Flynn et al (2000), who find little evidence for a population equivalent to that of Ibata et al in the Luyten proper motion survey. They conclude that hydrogen atmosphere white dwarfs cannot comprise a significant fraction of the Galactic halo. On the other hand, a new search of wide field photographic plates by Ibata et al (2000) has uncovered several high proper motion objects, at least two of which are spectroscopically confirmed to be white dwarfs and with some modicum of IR flux suppression.

In my not-completely-unbiased interpretation of the above, the evidence in favour of the white dwarf hypothesis is encouraging, although not conclusive. The recent results of Hodgkin et al and Ibata et al (2000) have demonstrated at least the existence of a population of high proper motion (and thus probably halo) white dwarfs. The question is now a matter of number density. It is on this point that we must await further data.

### 3.2. Future prospects

While the above results indicate that some white dwarfs do show halo kinematics, the fact that there is little evidence for a large halo white dwarf population in existing large-scale proper motion surveys casts some doubt that this population may be responsible for the microlensing results. This last statement relies heavily on the data set of Luyten (1979), so it would be nice to check this with another large scale survey.

Towards this end, a program led by Peter Stetson has begun at CFHT to image almost 20 square degrees to $`V25`$. Repeated with an interval of 2-3 years, this program should provide a conclusive test of the white dwarf halo hypothesis. Using the new MACHO mass fraction estimate ($`20\%`$) and assuming only half of all white dwarfs have hydrogen atmospheres, this program is still expected to find about 8 white dwarfs per square degree.

Another fascinating prospect is the hope of finding such old white dwarfs in the Sloan Digital Sky Survey (SDSS). Although this comprehensive large scale survey will not provide proper motions as a matter of course, some proper motion information is potentially available from comparison with Palomar Sky Survey plates. Furthermore, there will be a strip of sky in the southern hemisphere that will be multiply imaged. Even without proper motions, very old white dwarfs are potentially detectable by colours alone in the SDSS five-band photometric system. Photometric selection alone is unlikely to provide a complete census as only the oldest white dwarfs will stand out from the locus of other stars (the problem is that even the reddest bandpass doesn’t extend significantly beyond $`1\mu `$m). Nevertheless, very old white dwarfs are interesting in their own right.

Figure 2 shows a representative stellar locus and how one might select old white dwarfs based on colours. The colour selection suggests that only white dwarfs cooler than 3500 K can be selected in this manner. Using an i-band detection limit of 18, we find such objects can be detected to distances $`30`$ pc. If the entire MACHO fraction (taken to be 20%) were in 0.5$`M_{}`$ white dwarfs, 50% of whom had hydrogen atmospheres, and all had 3500 K temperatures, we would expect approximately 43 white dwarfs in the Sloan survey, or 1 every 480 square degrees. This is probably an overestimate. Hotter white dwarfs cannot be distinguished from the stellar locus (it is worth noting that the Hodgkin et al and Ibata et al detections all correspond to $`T_{\mathrm{eff}}35004000`$ K, so in the marginal region) and cooler objects are fainter and thus the effective volume is smaller. Nevertheless, this simple estimate offers encouragement that at least some white dwarfs are potentially detectable in this manner. More detailed calculations are underway to provide a more robust estimate.

Another large scale proper motion survey that should offer interesting constraints is that being undertaken by the EROS project (Goldman 1998). With a survey area of 350 square degrees and an I-limit of 20.5, we can use the same approximate scenario as above to estimate $`23`$ white dwarfs in their sample. This is a somewhat more realistic ballpark number in this case as the proper motion selection will allow detection of hotter objects.

The bottom line of the above is that, although uncertainty about mass distributions and chemical compositions make prognostication difficult, the fact that naive estimates lead to predictions of 10-100 detections in ongoing surveys suggests that we should see something if the hypothesis of a significant white dwarf halo is correct.

## 4. Beige Dwarfs

Although it was not my intention to address this subject at the conference, it was raised a couple of times, so I will conclude with a few remarks about “Beige” dwarfs (Hansen 1999b). Given that brown dwarfs are strongly ruled out by the microlensing timescales, red dwarfs by direct observations and white dwarfs (as well as neutron stars) are an uncomfortable fit due to issues of chemical pollution, it appears that there are no surviving baryonic candidates for halo MACHOs. However, there is one remaining possibility. Lenzuni, Chernoff & Salpeter (1992) demonstrated that one could circumvent the traditional hydrogen burning limit by starting with a brown dwarf and accreting material slowly enough that it could cool. In this fashion one could construct a degenerate hydrogen/helium object with mass $`>0.1M_{}`$. If one could build such objects to masses $`>0.3M_{}`$, they would make ideal MACHO candidates as they would be faint and would have no nuclear burning history and thus no chemical pollution problem. The term “Beige” dwarf comes from the superposition of brown dwarf and white dwarf characteristics. The primary problem for this scenario is that there is little evidence that this mode of “star” formation was ever important.

The existence of Beige dwarfs is observationally testable. They have similar radii to white dwarfs but can only exist at effective temperatures $`T_{\mathrm{eff}}<2000`$ K, so the detection of hotter objects would suggest white dwarfs. Indeed, the spectroscopically confirmed objects of Hambly et al and Ibata et al (2000) indicate temperatures $`35004000`$ K. Therefore these are unlikely to be Beige dwarfs.

#### Acknowledgments.

I would like to thank the conference organisers for a thoroughly enjoyable (and efficiently run) meeting. Support for this work was provided by NASA through Hubble Fellowship grant #HF-01120.01-99A, from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.

## References

Bahcall, J. N., Flynn, C., Gould, A. & Kirhakos, S. 1994, ApJ, 435, L51

Bergeron, P., Ruiz, M.-T. & Leggett, S. K. 1997, ApJS, 108, 339

Bergeron, P., Saumon, D. & Wesemael, F. 1995, ApJ, 443, 764

Flynn, C., Sommer-Larsen, J., Fuchs, B., Graff, D. S. & Salim, S. 2000, astro-ph/9912264

Goldman, B., in The Third Stromlo Symposium: The Galactic Halo, eds. Gibson, B. K., Axelrod, T. S. & Putman, M. E., ASP Conference Series Vol 165, p 413; San Francisco

Graff, D.S. & Freese, K. 1996, ApJ, 435, L516

Hansen, B. M. S. 1998, Nature, 394, 860

Hansen, B. M. S., 1999a, ApJ, 520, 680

Hansen, B. M. S., 1999b, ApJ, 517, L39

Harris, H. C., Dahn, C. C., Vrba, F. J., Henden, A. A., Liebert, J., Schmidt, G. D. & Reid, I. N. 1999, ApJ, 524, 1000

Hodgkin, S. T., Oppenheimer, B. R., Hambly, N. C., Jameson, R. F., Smartt, S. J. & Steele, I. A. 2000, Nature, 403, 57

Ibata, R. A., Richer, H. B., Gilliland, R. L. & Scott, D. 1999, ApJ, 524, L95

Ibata, R. A., Irwin, M., Bienayme, O., Scholz, R. & Guibert, J. 2000, ApJ, 532, L41

Lenzuni, P., Chernoff, D. F. & Salpeter, E. E., 1992, ApJ, 393, 232

Liebert, J., Dahn, C. C., Gresham, M. & Strittmatter, P. A. 1979, ApJ, 233, 226

Luyten, W. J. 1979, LHS Catalogue: A catalogue of stars with proper motions exceeding 0.5” annually, University of Minnesota, Minneapolis

Mendez, R. A. & Minniti, D. 2000, ApJ, 529, 911

Richer, H. B., Hansen, B., Limongi, M., Chieffi, A., Straniero, O. & Fahlman, G. G. 2000, ApJ, 529, 318

Salaris, M., Garcia-Berro, E., Hernanz, M., Isern, J. & Saumon, D. 2000, submitted to ApJ

Saumon, D. & Jacobsen, S. B. 1999, ApJ 511, L107

York, D. et al 2000, AJ submitted