text:
# Anthropic Reasons for Non-zero Flatness and Lambda

## I Introduction

The collection of considerations now known as the Anthropic Principles emerged from attempts by Whitrow  to understand why it is unsurprising that we find space to have three dimensions, and by Dicke  to understand the inevitability of Dirac ’large number’ coincidences in cosmology. Dicke recognised that it was unnecessary to introduce the idea of a time-varying gravitational constant in order to understand why we could not fail to observe that the number of protons in the observable universe is of order the square of the ratio of electromagnetic to gravitational force strengths. Subsequently, Dicke inspired a detailed observational and theoretical investigation of gravity theories in which the Newtonian gravitational constant becomes a space-time variable. he was partly motivated by apparent discrepancies between the predictions of standard general relativity and observations of the perihelion precession of Mercury. These discrepancies were subsequently ascribed to errors in the measurements of the shape and diameter of the Sun created by solar surface activity .

There have been many investigations of the apparent coincidences that allow complexity to exist in the universe (see ). Typically, they examine the stability of life-supporting conditions to small (or large) perturbations to the values of constants of Nature or to quantities fixed by cosmological ’initial’ conditions at $`t=0`$ or $`t=\mathrm{}`$. These in turn divide into studies of two sorts: first, those in which the hypothetical changes introduced to the ’constants’ are self-consistently permitted by the cosmological or physical theory employed; and second, those in which they are not. An investigation of the first kind might be one in which the cosmological initial conditions were enlarged to allow anisotropies or the possibility of a significant deviation from flatness. An investigation of the second type might note that a change in the observed value of the electron to proton mass ratio to another fixed value would make it difficult to produce ordered molecular structures. Studies of universes in which traditional ’constants’ of Nature are changed are restricted by the lack of self-consistent theories which allow all these possible changes to be accommodated. Without them, it is impossible to determine the possible knock-on effects of varying one constant on others.

There are some exceptions. Varying gravitation ’constant’, $`G`$, (or dimensionless constants formed with it like $`Gm^2/hc`$ for any mass $`m`$), can be studied using scalar-tensor gravity theories . A varying fine structure ’constant’ can be studied using the theory of Bekenstein and Sandvik, Barrow and Magueijo (BSBM) , . Moreover, the formulation of physical theories whose true constants inhabit more than three space dimensions provides a framework for the rigorous study of the simultaneous variation of their three-dimensional counterparts , , . Recently there has also been much interest in theories where a variation in the fine structure constant is due to a change in the light propagation speed. In another paper we propose various methods for experimentally distinguishing between these different theories.

Observational evidence for a variation in a traditional constant can be found without the need for a self-consistent theory of its variation simply by demonstrating incompatibility with the predictions of the standard theory. The most observationally sensitive ’constant’ is the electromagnetic fine structure constant, $`\alpha e^2/\mathrm{}c`$, and recent observations motivate the formulation of varying-$`\alpha `$ theories. The new many-multiplet technique of Webb et al, , , exploits the extra sensitivity gained by studying relativistic transitions to different ground states using absorption lines in quasar (QSO) spectra at medium redshift. It maximises the information extracted from the data set and has provided the first evidence that the fine structure constant might change with cosmological time. The trend of these results is that the value of $`\alpha `$ was lower in the past, with $`\mathrm{\Delta }\alpha /\alpha =0.72\pm 0.18\times 10^5`$ for $`z0.53.5.`$ Other investigations have claimed preferred non-zero values of $`\mathrm{\Delta }\alpha <0`$ to best fit the cosmic microwave background (CMB) and Big Bang Nucleosynthesis (BBN) data at $`z10^3`$ and $`z10^{10}`$ respectively but appeal to much larger variations. We have shown that the simplest theory which joins varying $`\alpha `$ to general relativity via the propagation of a scalar field can explain these observations together with the lack of evidence for a similar level of variation locally, 2 billion years ago, or at very high redshifts, $`z10^3`$ . In this paper we will show how this theory also provides some novel anthropic perspectives on the evolution of our universe or others.

There have been several studies, following Carter,  and Tryon , of the need for life-supporting universes to expand close to the ’flat’ Einstein de Sitter trajectory for long periods of time. This ensures that the universe cannot collapse back to high density before galaxies, stars, and biochemical elements can form by gravitational instability, or expand too fast for stars and galaxies to form by gravitational instability (see also ,  and ). Likewise, it was pointed out by Barrow and Tipler,  that there are similar anthropic restrictions on the magnitude of any cosmological constant, $`\mathrm{\Lambda }`$. If it is too large in magnitude it will either precipitate premature collapse back to high density (if $`\mathrm{\Lambda }<0`$) or prevent the gravitational condensation of any stars and galaxies (if $`\mathrm{\Lambda }>0`$). Thus existing studies provide anthropic reasons why we can expect to live in an old universe that is neither too far from flatness nor dominated by a much stronger cosmological constant than observed ($`\left|\mathrm{\Lambda }\right|10\left|\mathrm{\Lambda }_{obs}\right|`$).

Inflationary universe models provide a possible theoretical explanation for proximity to flatness but no explanation for the smallness of the cosmological constant. Varying speed of light theories  offer possible explanations for proximity to flatness and smallness of a classical cosmological constant (but not necessarily for one induced by vacuum corrections in the early universe). Here, we shall show that if we enlarge our cosmological theory to accommodate variations in some traditional constants then it appears to be anthropically disadvantageous for a universe to lie too close to flatness or for the cosmological constant to be too close to zero. This conclusion arises because of the coupling between time-variations in constants like $`\alpha `$ and the curvature or $`\mathrm{\Lambda }`$, which control the expansion of the universe. The onset of a period of $`\mathrm{\Lambda }`$ or curvature domination has the property of dynamically stabilising the constants, thereby creating favourable conditions for the emergence of structures. This point has been missed in previous studies because they have never combined the issues of $`\mathrm{\Lambda }`$ and flatness and the issue of the values of constants. By coupling these two types of anthropic considerations we find that too little $`\mathrm{\Lambda }`$ or curvature can be as poisonous for life as too much.

## II Time variation of $`\alpha `$

First, consider a simple theory with varying $`\alpha e^{2\psi }/\mathrm{}c`$ where $`\psi `$ is a scalar field that can vary in space and time. A generalisation of the scalar theory proposed by Bekenstein  described in ref.  to include the gravitational effects of $`\psi `$ gives the field equations

$$G_{\mu \nu }=8\pi G\left(T_{\mu \nu }^{matter}+T_{\mu \nu }^\psi +T_{\mu \nu }^{em}e^{2\psi }\right),$$
(1)
and the $`\psi `$ field obeys the equation of motion

$$\mathrm{}\psi =\frac{2}{\omega }e^{2\psi }_{em}.$$
(2)
We have defined the coupling constant $`\omega =(\mathrm{}c)/l^2`$, where $`l`$ is the length scale down to which the theory is accurately coloumbic. It is clear that $`_{em}`$ vanishes for a sea of pure radiation since then $`_{em}=(E^2B^2)/2=0`$. We therefore expect the variation in $`\alpha `$ to be driven by electrostatic and magnetostatic energy-components rather than electromagnetic radiation. In order to make quantitative predictions we need to know how much of the non-relativistic matter contributes to the RHS of Eqn. (2). This is parametrised by $`\zeta _{em}/\rho `$, where $`\rho `$ is the energy density , and for baryonic matter $`_{em}=E^2/2`$. In previous papers  we showed how the cosmological value of $`\zeta `$ (denoted $`\zeta _m`$) is largely determined by the nature of dark matter. To accommodate for a lower $`\alpha `$ in the past, as preferred by the data, the dark matter constituents need to have high magnetostatic energy content (One possible contender would be superconducting cosmic strings which have $`\zeta _m1`$). In line with our recent work and the observational data we will in this paper confine ourselves to negative values of $`\zeta _m`$.

Assuming a homogeneous and isotropic Friedmann metric with expansion scale factor $`a(t)`$ and curvature parameter $`k`$ we obtain the field equations ($`c1`$)

$$\left(\frac{\dot{a}}{a}\right)^2=\frac{8\pi G}{3}\left(\rho _m\left(1+|\zeta _m|e^{2\psi }\right)+\rho _re^{2\psi }+\frac{\omega }{2}\dot{\psi }^2+\rho _\mathrm{\Lambda }\right)\frac{k}{a^2},$$
(3)
where the cosmological vacuum energy $`\rho _\mathrm{\Lambda }`$ is a constant that is proportional to the cosmological constant $`\mathrm{\Lambda }8\pi G\rho _\mathrm{\Lambda }`$. For the scalar field we have

$$\ddot{\psi }+3H\dot{\psi }=\frac{2}{\omega }e^{2\psi }\zeta _m\rho _m$$
(4)
where $`H\dot{a}/a.`$ The conservation equations give for the non-interacting radiation, and matter densities $`\rho _r`$ $`e^{2\psi }a^4`$and $`\rho _ma^3,`$ respectively. This theory enables the cosmological consequences of varying $`\alpha `$, to be analysed self-consistently rather than by changing the constant value of $`\alpha `$ in the standard theory, as in the original proposals made in response to the large numbers coincidences .

The cosmological behaviour of the solutions to these equations was studied by us ,  for the $`k=0`$ case and is shown in Figure (1). The evolution of $`\alpha `$ is summarised as follows:

1. During the radiation era $`\alpha `$ is constant and $`a(t)t^{1/2}`$. It increases in the dust era, where $`a(t)t^{2/3}`$, until the cosmological constant starts to accelerate the universe, $`a(t)\mathrm{exp}[\mathrm{\Lambda }t/3],`$ after which $`\alpha `$ asymptotes rapidly to a constant, see fig.(1)

2. If we set the cosmological constant equal to zero then, during the dust era, $`\alpha `$ will increase indefinitely. The increase however, is very slow with a late-time solution for $`\psi `$ proportional to $`\mathrm{log}(2N\mathrm{log}(t))`$, see fig.(2). N is defined as $`N2\zeta _m/\rho _ma^3`$, a positive constant since we have confined ourselves to $`\zeta _m<0`$.

3. If we set the cosmological constant equal to zero and introduce a negative spatial curvature ($`k<0`$) then $`\alpha `$ increases only during the dust-dominated phase, where $`a(t)t^{2/3}`$, but tends to a constant after the expansion becomes curvature dominated, with $`a(t)t`$. This case is illustrated in fig.(3).

From these results it is evident that non-zero curvature or cosmological constant brings to an end the increase in the value of $`\alpha `$ that occurs during the dust-dominated era***In some Friedmann universes with initial conditions unlike our own there can be power-law growth of $`\alpha `$ during the radiation era, . In such universes the same general effects of negative curvature and positive $`\mathrm{\Lambda }`$ are seen. They still halt any growth in $`\alpha (t).`$ Our initial conditions are chosen so as to give a present day value of $`\alpha 1/137`$. The initial value of alpha would have to be several orders of magnitude lower in order to obtain the power-law growth. Hence, if the spatial curvature and $`\mathrm{\Lambda }`$ are too small it is possible for the fine structure constant to grow too large for biologically important atoms and nuclei to exist in the universe. There will be a time in the future when $`\alpha `$ reaches too large a value for life to emerge or persist. The closer a universe is to flatness or the closer $`\mathrm{\Lambda }`$ is to zero so the longer the monotonic increase in $`\alpha `$ will continue, and the more likely it becomes that life will be extinguished. Conversely, a non-zero positive $`\mathrm{\Lambda }`$ or a non-zero negative curvature will stop the increase of $`\alpha `$ earlier and allow life to persist for longer. If life can survive into the curvature or $`\mathrm{\Lambda }`$-dominated phases of the universe’s history then it will not be threatened by the steady cosmological increase in $`\alpha `$ unless the universe collapses back to high density.

## III Anthropic Limits on $`\alpha `$

We have seen that varying-$`\alpha `$ cosmologies with zero curvature and $`\mathrm{\Lambda }`$ lead to a monotonic increase in $`\alpha `$ with time. Here we summarise the principal upper limits on $`\alpha `$ that are needed for atomic complexity and stars to exist. There are a variety of constraints on the maximum value of the fine structure compatible with the existence of nucleons, nuclei, atoms and stars under the assumption that the forms of the laws of Nature remain the same. The running of the fine structure constant with energy due to vacuum polarisation effects leads to an exponential sensitivity of the proton lifetime with respect to the low-energy value of $`\alpha `$ with $`t_{pr}\alpha ^2\mathrm{exp}(\alpha ^1)m_{pr}^110^{32}yrs.`$ In order that the lifetime be less than the main sequence lifetime of stars we have $`t_{pr}<(Gm_{pr}^2)^1m_{pr}^1`$ which implies that $`\alpha `$ is bounded above by $`\alpha <1/80`$ approximately .

The stability of nuclei is controlled by the balance between nuclear binding and electromagnetic surface forces . A nucleus ($`Z,A`$) will be stable if $`Z^2/A<49(\alpha _s/0.1)^2(1/137\alpha ).`$ In order for carbon $`(Z=6)`$ to be stable we require $`\alpha <16(\alpha _s/0.1)^2.`$ Detailed investigations of the nucleosynthesis processes in stars have shown that a change in the value of $`\alpha `$ by $`4\%`$ shifts the key resonance level energies in the carbon and oxygen nuclei which are needed for the production of a mixture of carbon and oxygen from beryllium plus helium-4 and carbon-12 plus helium-4 reactions in stars . These upper bounds on $`\alpha `$ are model independent and were considered in more detail in refs. ,  and . However, sharper limits can be found by using our knowledge of the stability of matter derived from analysis of the Schrödinger equation. Stability of matter with Coulomb forces has been proved for non-relativistic dynamics, including arbitrarily large magnetic fields, and for relativistic dynamics without magnetic fields. In both cases stability requires that the fine structure constant be not too large.

The value of $`\alpha `$ controls atomic stabilityNote that if the electron mass and velocity of light are varied along with the value of $`\alpha `$ then the eigenvalues of the non-relativistic Schrödinger equation can remain invariant and atomic structure is unchanged . Here, we break the scale invariance by varying only $`\alpha .`$. If $`\alpha `$ increases in value then the innermost Bohr orbital contracts and electrons will eventually fall into the nucleus when $`\alpha >Z^1m_{pr}/m_e.`$ As $`\alpha `$ increases, atoms all become relativistic and unstable to pair production. In order that the electromagnetic repulsion between protons does not exceed nuclear strong binding $`e^2/r_n<\alpha m_\pi `$ is needed and so we require $`\alpha <1/20.`$ It is also known that atomic instability of atoms with atomic number $`Z`$ occurs in the relativistic Schrödinger equation when the fine structure constant is increased in value to $`\alpha =\frac{2}{\pi Z}.`$ However, when the many-electron and many-nucleon problem is examined with the relativistic Schrödinger theory there is a bound on $`\alpha `$ for stability that is independent of $`Z`$ . If $`\alpha <1/94`$ then stability occurs all the way up to the critical value $`\alpha =\frac{2}{\pi Z}`$, whereas if $`\alpha `$ $`>128/15\pi `$ the ’atomic’ system is unstable for all values of $`Z`$. In the presence of arbitrarily large magnetic fields, which aid binding by creating a two-dimensional form for the potential, matter composed of electrons and nuclei is known to be unstable if $`\alpha `$ or $`Z`$ is too large: matter is stable if $`\alpha `$ $`<0.06`$ and $`\alpha `$ $`<0.026(6/Z)^{1/2}`$,,.

If stars are to exist, their centres must be hot enough for thermonuclear reactions to occur. This requires $`\alpha `$ to be bounded above by $`\alpha ^2<20m_e/m_{pr}.`$ Carter has also pointed out the existence of a very sensitive condition $`\alpha ^{12}\left(m_e/m_{pr}\right)^4Gm_{pr}^2,`$that must be met if stars are to undergo a convective phase, although this stringent condition no longer seems to be essential for planetary formation .

The results collected above show that there are a number of general upper limits on the value of $`\alpha `$ if atoms, molecules, and biochemistry are to exist. These bounds do not involve the gravitation constant explicitly. Other astrophysical upper bounds on $`\alpha `$ exist in order that stars be able to form but these involve the gravitational constant.

## IV Time variation of $`G`$

A similar trend can be found in relativistic cosmologies in scalar-tensor gravity theories. Consider the paradigmatic case of Brans-Dicke (BD) theory to fix ideas. The form of the general solutions to the Friedmann metric in BD theories are fully understood ,. The general solutions begin at high density dominated by the BD scalar field $`\varphi G^1`$ and approximated by the vacuum solution. At late times they approach particular exact power-law solutions for $`a(t)`$ and $`\varphi (t)`$ and the evolution is ’Machian’ in the sense that the cosmological evolution is driven by the matter content rather than by the kinetic energy of the free $`\varphi `$ field. There are three essential field equations for the evolution of $`\varphi `$ and $`a(t)`$ in a BD universe

$`3{\displaystyle \frac{\dot{a}^2}{a^2}}`$ $`=`$ $`{\displaystyle \frac{8\pi \rho }{\varphi }}3{\displaystyle \frac{\dot{a}^{}\dot{\varphi }}{a\varphi }}+{\displaystyle \frac{\omega _{BD}}{2}}{\displaystyle \frac{\dot{\varphi }^2}{\varphi ^2}}{\displaystyle \frac{k}{a^2}}`$
$`\ddot{\varphi }+3{\displaystyle \frac{\dot{a}}{a}}\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{8\pi }{3+2\omega }}(\rho 3p)`$
$`\dot{\rho }+3{\displaystyle \frac{\dot{a}}{a}}(\rho +p)`$ $`=`$ $`0`$
Here, $`\omega _{BD}`$ is the BD constant parameter and the theory reduces to general relativity in the limit $`\omega _{BD}\mathrm{}`$ and $`\varphi =G^1`$constant.

In the radiation era the scale factor approaches the standard general relativistic behaviour for large times:

$$a(t)t^{1/2};G=constant$$
(5)
After the dust density dominates the dynamics the expansion approaches a simple exact solution with

$$a(t)t^{(2n)/3};Gt^n,$$
(6)
which continues until the curvature term takes over the expansion. Here, $`n`$ is related to the constant Brans-Dicke $`\omega _{BD}`$ parameter by

$$n\frac{2}{4+3\omega _{BD}}$$
(7)
and the usual general relativistic Einstein de Sitter universe is obtained as $`\omega _{BD}\mathrm{}`$ and $`n0.`$ If the universe is open, ($`k=1`$), then the negative curvature will eventually dominate the gravitational effects of the dust and then the BD model approaches the general relativistic Milne model with constant $`G`$

$$a(t)t;G=constant$$
(8)
Again, we see the same pattern of behaviour seen for the evolution of $`\alpha `$ in the BSBM theory. The smaller the curvature term, so the longer the dust-dominated era lasts, and the greater the fall in the value of $`G`$, and the smaller its ultimate asymptotic value when the curvature intervenes to turn off the variation. In general, in such cosmologies, if there exists a critical value of $`G`$ below which living complexity cannot be sustained, then a universe that is too close to flatness will have a smaller interval of cosmic history during which it can support life.

So far, we have discussed only the independent variation of $`\alpha `$ and $`G`$. What happens if they both vary at the same time? Previous studies of varying constants have only examined the time-variation of a single ’constant’. We have produced a unified theory , which incorporates the BSBM varying $`\alpha `$ and BD varying $`G`$ theories discussed above. When both $`\alpha `$ and $`G`$ are allowed to vary simultaneously in this theory we find  that our general conclusions still hold, although the quantitative details are changed. During the dust era of a flat Friedmann universe with varying $`\alpha (t)`$ and $`G(t),`$their time-evolution approaches an attractor in which the product $`\alpha G`$ is a constant and

$$\alpha G^1t^n$$
(9)
where $`n`$ is given by eq. (7). Thus we see that the $`G`$ evolution is left unchanged by the effects of varying $`\alpha `$, but variation of $`G`$ changes the time evolution of $`\alpha (t)`$ from a logarithm to a power-law in time. As before, the longer the dust era lasts before it is ended by deviation from flatness or zero cosmological constant, the longer the time-increase of $`\alpha `$ continues, inevitably leading to values that make any atom-based complexity impossible.

## V  Discussion

We have shown that some theories which include the time variation of traditional constants like $`\alpha `$ and $`G`$ introduce significant new anthropic considerations. A theory which self-consistently introduces the space-time variation of a traditional constant scalar quantity is strongly constrained in form by the requirements of causality and second-order propagation equations . Typically, this requirement leads to equations for the driving scalar, $`\phi `$ that have the form $`\mathrm{}\phi `$ proportional to linear combinations of the energy-momentum components. Explicit examples are provided by the Bekenstein-Sandvik-Barrow-Magueijo and Brans-Dicke theories. This structure ensures that the evolution of the ’constant’ whose variations are derived from those of $`\phi `$ is strongly dependent upon the material or geometrical source governing the background expansion dynamics. In the case of varying $`\alpha `$ we have shown elsewhere ,  that this ties the epoch after which time-variations in $`\alpha `$ become very small to the time when the cosmological constant starts to accelerate the expansion of the universe. In these theories there is therefore the possibility of a habitable time zone of finite duration during which a constant like $`\alpha `$ or $`G`$ falls within a biologically acceptable range.

Surprisingly, there has been almost no consideration of habitability in cosmologies with time-varying constants since Haldane’s discussions  of the biological consequences of Milne’s bimetric theory of gravity with two timescales, one for atomic phenomena, another for gravitational phenomena . Since then attention has focussed upon the consequences of universes in which the constants are different but still constants. Those cosmologies with varying constants that have been studied have not considered the effects of curvature or $`\mathrm{\Lambda }`$ domination on the variation of constants and have generally considered power-law variation to hold for all times. The examples described here show that this restriction has prevented a full appreciation of the coupling between the expansion dynamics of the universe and the values of the constants that define the course of local physical processes within it. Our discussion of a theory with varying $`\alpha `$ shows for the first time a possible reason why the 3-curvature of universes and the value of any cosmological constant may need to be bounded below in order that the universe permit atomic life to exist for a significant period. Previous anthropic arguments have shown that the spatial curvature of the universe and the value of the cosmological constant must be bounded above in order for life-supporting environments (stars) to develop. We note that the lower bounds discussed here are more fundamental than these upper bounds because they derive from changes in $`\alpha `$ which have direct consequences for biochemistry whereas the upper bounds just constrain the formation of astrophysical environments by gravitational instability (for alternative scenarios see ref. ). Taken together, these arguments suggest that within an ensemble of all possible worlds where $`\alpha `$ and $`G`$ are time variables, there might only be a finite interval of non-zero values of the curvature and cosmological constant contributions to the dynamics that both allow galaxies and stars to form and their biochemical products to persist.

### a Acknowledgements

We would like to thank Bernard Carr, Andrei Linde, and Martin Rees for discussions. HBS acknowledges support from the Research Council of Norway. The numerical work shown herein was performed on COSMOS, the Origin 2000 supercomputer owned by the UK-CCC.


format: text

text:
# Theoretical modeling of prion disease incubation

## I INTRODUCTION

Understanding the factors which regulate the incubation times for infectious prion diseases is important for assessing the risk of illness after potential exposure as well as for developing treatments which can delay disease onset. There are several striking aspects to prion disease incubation, which are not well understood: (i) The incubation times can run into years and decades , and yet, at the laboratory scale have been found to be highly reproducible. In fact, the reproducibility of incubation times with dose has been used as an independent measure of infectivity titre . (ii) There seem to be distinct stages for the disease incubation after intracerebral inoculation: (a) Following rapid initial clearance, there is a ‘lag’ phase (also termed ‘zero’ phase ) during which there is little or no infectivity, and (b) an exponential growth or doubling phase, during which the infectivity increases exponentially with a well-defined doubling period . Understanding the lag-phase is clearly important as any treatment strategy is more likely to succeed before the exponential growth phase takes over. (iii) As the dose of infection is increased in the laboratory, the incubation times become sharply defined and saturate to a dose independent value, but as the dose is reduced a broad distribution and a logarithmic dose dependence results . Such a broad distribution has also been found in epidemiological studies of Bovine Spongioform Encephelopathy (BSE) in England . (iv) For infection across species, there is a ‘species barrier’ and the first passage takes considerably longer to incubate than subsequent passages . (v) While prion aggregation has been observed in vitro, the aggregates are neuro-toxic but not infectious . These are the issues which motivate our study.

The purpose of this paper is to test the extent to which a purely physico-chemical model can capture the main reproducible features of prion disease incubation. In particular, we emphasize the importance of the aggregate morphology in determining the statistics of incubation times. Our basic hypothesis is that the ‘lag phase’ is determined by growth of misfolded protein-aggregates from initial small seeds (acquired through infection) to a typical ‘fissioning dimension’, whereas subsequent aggregate-fissioning and spreading leads to exponential growth and the doubling phase. For a single seed, the lag phase develops a broad but well defined distribution, which we can calculate via a microscopic statistical model. Thus, when the infection is very dilute, there is a broad distribution of incubation times. At higher doses of infection, self-averaging due to independent growth from many seeds leads to sharply defined incubation times. The dose dependence and its saturation, as well as the ratio of lag time to doubling time depends on the morphology of the aggregates, i.e., whether one has linear fibrils or compact higher-dimensional aggregates. In this sense, details of incubation time distributions provide an indirect means to infer early growth morphologies. Alternative theoretical models which deal with the above issues have also been developed in the literature.

The extent to which such a model explains the experimental phenomenology would help address the following questions: (i) Are the incubation times dominated by the nucleation and growth of misfolded protein aggregates?. (ii) Are the two phases of prion disease incubation, the lag phase and the exponential growth phase, controlled by the same process i.e. aggregation of misfolded proteins?  (iii) Assuming that aggregate growth controls the incubation time scales, we are led to ask: What is the aggregate morphology during early growth and how does it influence the dose incubation curves (iv) Does current data inform us about the characteristic size of the aggregates? (v) What are the practical epidemiological impacts from our model incubation time distributions? Can we use the model to constrain estimates of the number of infections of vCJD for example?

Employing statistical simulations of prion aggregation (based upon cellular automata rules) we argue here that several features of prion disease can be explained by exploring the statistics of the two phases of the disease incubation. Within our model, we show that compact two-dimensional aggregates can provide the observed broad distribution of incubation times for dilute doses and at the same time account for the typical large difference between lag time and doubling time. We present analytic calculations which provide a functional form for the distribution that can be used in further epidemiological studies, and we use these results to infer the dose dependence of the incubation time. Furthermore, our analysis shows how the dose-incubation curve can be related to experimental measurements of the time-course of infectivity and in particular, testable predictions can be made using our model for the dose dependence of the lag phase. Finally, we apply our model to epidemiological data for BSE (mad cow disease), iatrogenic CJD infections associated with dura mater transplants, and a vCJD cluster from the United Kingdom, from which we conclude that in all three cases the incubation time is dominated by slow aggregation from a few small (tens of nanometer scale) starting seeds. This size scale is comparable to small animals, but the estimated attachment rates are slower by an order of magnitude or more. We argue that the incubation time for vCJD is likely to be at the low end of previous estimates, implying an infectious toll in the hundreds rather than hundreds of thousands. We speculate that the differing attachment rate among species is regulated by the PrPc concentration which in turn is metabolically controlled.

We organize our paper as follows: Section II discusses a model incubation time distribution, which illustrates how our basic picture relates to dose incubation curves in prion diseases. Sections III-V deal with microscopic models related to protein misfolding, aggregation and fissioning. In Sections VI, VII, and VIII, we come back to the dose incubation curves, connections to epidemiological data and the disease phenomenology and discuss them in context of our models and present our conclusions.

## II A Model Distribution and Dose Incubation Curves

We first illustrate the ‘bare bones’ of our proposed picture for prion disease incubation by using a model distribution, where calculations can be done analytically. As mentioned in the previous section, the lag phase corresponds to aggregation of misfolded prions from the initial seed upto a ‘fissioning size’. This is a stochastic process and correspondingly there will be a distribution of aggregation times. Subsequent to this aggregation, we assume (based upon experimental observations) that the number of seeds increases exponentially with a well defined doubling time ($`t_2`$). This exponential growth continues until the number of seeds reaches a critical value which signals the onset of clinical symptoms and the end of the incubation period.

In this section, we assume the distribution of aggregation times from a single seed, $`P(t)`$, to be given by

$`P(t)`$ $`=`$ $`0t<t_0`$ (1)
$`=`$ $`{\displaystyle \frac{1}{n_1t_2}}t_0<t<t_0+n_1t_2`$ (2)
$`=`$ $`0t>t_0+n_1t_2`$ (3)
Here $`t_0`$ and $`n_1`$ are the parameters for the distribution. It has a sharp onset at time $`t_0`$ and its width is taken to be some multiple ($`n_1`$) of the doubling time $`t_2`$. When there are many seeds present, each initial seed will start fissioning into two new seeds after an aggregation time according to the above distribution. We will assume that once a seed has fissioned once, it continues to fission or effectively double every $`t_2`$ time steps, which is independent of the above distribution. The microscopic basis for this assumption will be explained in section V. The incubation time is given by the mean time taken for a given initial number of seeds ($`D_i`$) to reach a critical number ($`D_f`$) following the above processes.

The dose dependence of the the incubation time is calculated through the following steps:

1) First we calculate the mean ‘first arrival’ time i.e. the mean time taken for the first aggregate fissioning event. Let the cumulative probability for the ‘first arrival’ time for $`D_i`$ seeds be given by $`F^{(D_i)}(t)`$. Since each seed grows independently, this can be related to the cumulative probability for the ‘first arrival’ time for a single seed via the relation

$$F^{(D_i)}(t)=1(1F^{(1)}(t))^{D_i}$$
(4)
The mean first arrival time $`t_1`$ is given by solving $`F^{(D_i)}(t_1)=\frac{1}{2}`$. For the simple probability distribution discussed above, the mean first arrival time is well-approximated by the expression

$$t_m(D_i)=t_0+\frac{n_1}{2D_i}t_d$$
(5)
2) We now proceed to calculate the time spent in the doubling phase, i.e., the time taken till the number of seeds reaches $`D_f`$. All the aggregates formed after fissioning are assumed to further fission in time $`t_2`$. Besides these, initial seeds continue to ‘arrive’, i.e., aggregate to the fissioning size and thus join the number seeds that are doubling. Let the number of such seeds ‘arriving’ in the time interval $`mt_2<t<(m+1)t_2`$ be given by $`\alpha (m)D`$. Then, $`\alpha (m)`$ is given by

$$\alpha (m)=F^{(1)}((m+1)t_2)F^{(1)}(mt_2)$$
(6)
3) For our model distribution, the above equation can now be used to calculate the mean incubation time, which is approximately given by

$$t_i=t_0+\frac{n_1}{2D_i}t_2+[\mathrm{l}og_2(D_f)\mathrm{l}og_2(1+\frac{D_i}{n_1})]t_2$$
(7)
The above equation gives the dose-dependence of the incubation time for the model distribution. It should be noted that this expression already explains several generic features seen in experimental dose-incubation curves (DICs) and in the microscopic models we present in later sections. These are: 1) logarithmic dose-dependence at high doses and 2) deviations from logarithmic behavior at low doses. From the above expression for the incubation time, it can be seen that the dominant contribution to the incubation time at high $`D_i`$ comes from the ‘doubling’ phase which gives a logarithmic dose-dependence. However at low doses ($`Dn1`$), the time spent in the doubling phase does not change appreciably with dose. Instead, the variation in the incubation time comes from the dose-dependence of the the ‘first arrival’ time i.e. the lag phase time. This is reflected as a deviation from the logarithmic behavior in the DIC which is seen in experimental DICs . Our model thus makes the testable prediction that deviations from logarithmic behavior in the DIC should correspond to increases in the lag phase.

We now proceed to develop microscopic models for the initial aggregation process. The next section discusses the cellular automata approach to this problem.

## III Cellular Automata Simulations

Theoretical modeling of incubation times  starting at the molecular-level is all but impossible with a twenty order of magnitude span between molecular motion time scales and those of disease onset. On the other hand, kinetic theory allows one to model long-time processes but ignores short distance spatial fluctuations, important in nucleation and growth. We have developed a lattice-based protein-level cellular-automata approach, which bridges these two methodologies . Previously, we used it to calculate aggregation-time distributions, which compared favorably with the incubation-times inferred from BSE data . We also showed that playing with the rules in such simple models can be a “cheap” way to suggest, constrain and guide treatment protocols.

Our models consist of dilute concentrations of proteins diffusing on the lattice and interconverting between their properly folded state (PrPc) and the misfolded state (PrPSc) . Our key assumption is that the conformational state of a protein depends on the amount of water around it. A monomer isolated from others (surrounded by the omnipresent water) stays in its properly folded state. However, when proteins are surrounded by other proteins, thus excluding water from parts of their neighborhood, they can change conformations and go into a misfolded state (involving $`\beta `$ sheet bonding). A key parameter of our model is the coordination, $`q_c`$, at which the misfolded conformation PrPSc becomes stable. Only misfolded monomers may remain stably in a cluster, possibly breaking away from a cluster when they fold back into the PrPc form.

Assuming aggregation happens on the cell-surface, we choose a 2d hexagonal lattice. There is strong evidence that the local coordination environment of the sphingolipids to which the prion proteins attach is, in fact, hexagonal. The lattice structure and the detailed protein motion are not crucial in our model. At each time step, proteins can move randomly by at most a unit lattice spacing. The magnitude of the time-step is set by the time for a single monomer to misfold. It is implicitly assumed that proteins co-adsorb with each other followed either by a conversion in shape or separation. It is this conversion process which sets the unit of time.

By playing with the cellular automata rules it is possible to get different aggregate morphologies and aggregation time distributions. First, we consider the case where proteins are isotropic objects. We have performed a large number of runs at values of $`q_c=`$ 1, 2, and 3, with different monomer concentrations (held fixed during the simulation). The aggregation time is defined as the time required to grow from initial seed of size $`𝒜_i`$ to a final size $`𝒜`$. The lower coordination rules effectively remove the nucleation barrier, leading to frequent nucleation of new clusters. Typical aggregate configurations (and stable seeds) are shown in Fig. 1.

From these studies, we see that lowering the critical coordination provides too rapid a growth for prion aggregates, and with no nucleation barrier at these concentrations. In contrast, for $`q_c=3`$, the aggregation is very slow and it is characterized by a broad incubation time distribution, with a clear separation in time-scales for seeded and unseeded (i.e. infectious and sporadic) cases.

We can also obtain one dimensional fibril growth by considering the proteins to be anisotropic. For example, on a square lattice, we can get fibrils by: (i) identifying a preferred bonding face to our simple point proteins, now made into squares. (ii) We choose a critical coordination of 2. (iii) We make edge bonding of proteins with adjacent preferred faces to be quite strong under coordination q=1 (i.e., the conversion probability is higher than 50%), and somewhat less strong for face to face meeting of proteins. We assume zero conversion probabilities for all other faces. By choosing three kinds of faces with appropriate rules, we can obtain equivalent results on the hexagonal lattice. These rules assure fibril growth which is dimer dominated (see Fig. 1).

## IV Stochastic Analysis of Aggregation

In the low concentration limit, the aggregation results from a sequential addition of proteins to the initial seed. However, addition of monomers is not always stable. Given the rules, various stages of the aggregate size and shape require a pair of proteins (a dimer) to arrive simultaneously, in order to attach in a stable manner. Thus, the entire process can be approximated by one of stochastic sequential addition of monomer and dimer units. As the concentration, $`c`$, goes to zero, the monomer addition rate is proportional to $`c`$, whereas the dimer addition rate is proportional to $`c^2`$, and thus the growth will involve a minimum number of dimers and these will provide the dominant contribution to the growth times.

The growth to a final size $`𝒜`$ from an initial size $`𝒜_i`$ involves sequential addition of $`n`$ units. The probability for the successive additions at intervals $`t_1`$, $`t_2`$, $`\mathrm{}`$, $`t_n`$ is

$$P(t_1,t_2,\mathrm{},t_n)=\underset{j=1}{\overset{n}{}}p_je^{p_jt_j},$$
(8)
where the rate for the $`j`$th unit, $`p_j`$, depends on the geometry of the aggregate and the kind of unit (monomer or dimer) to be added. Hence, the probability distribution associated with the total growth time is

$$P(t)=_0^{\mathrm{}}𝑑t_1\mathrm{}_0^{\mathrm{}}𝑑t_n\underset{j=1}{\overset{n}{}}p_je^{p_jt}\delta (t\underset{i=1}{\overset{n}{}}t_i).$$
(9)
This integral can be evaluated by standard methods for arbitrary $`p_i`$. We note the answers for two cases:

1) The attachment probabilities are identical for each unit i.e $`p_j=p`$ for all $`j`$. In this case the probability distribution is the Gamma distribution:

$$P(t)=\frac{p(pt)^{n1}}{(n1)!}e^{pt}$$
(10)
2) The rate, $`p_j=p+jp^{}`$, increases linearly with $`j`$. In this case we obtain the Beta distribution in $`e^t`$ :

$$P(t)=Ae^{(p+p^{})t}(1e^{p^{}t})^{n1}$$
(11)
In 1d fibril growth, dimers are attached one by one with the area available for attachment of dimers staying constant, and thus Eq. 10 applies. For 2d compact growth, slow dimer attachments have to be combined with rapid filling up of rows by monomers. In the low concentration limit, the rate is limited by dimer attachment probabilities which increase linearly with the number of dimers already attached, thus leading to Eq. 11.

At finite concentrations, the monomer attachment times can no longer be neglected, and a more accurate treatment of the time scales in the 2d case requires a convolution of probabilities for monomer attachment times with those for dimer attachment times. The geometrical counting of number of monomers and number of dimers needed to grow to the desired size is straightforward and can be used to develop accurate fits to the numerical data (See Figure 2 A).

An important aspect of our 2d model is the asymptotic compression of the distributions at low concentrations. The initial stage of the growth is extremely slow and the process speeds up significantly as the aggregate grows. Thus, the mean aggregation time, $`t_m(1)`$ to go from an initial seed $`𝒜_i`$ to a final size $`𝒜`$ can be much larger than the typical aggregate-doubling time $`t_2`$ to go from size $`𝒜/2`$ to $`𝒜`$. Fig. 3 shows the ratio $`t_2(1)/t_m(1)`$ for different concentrations and different final sizes $`𝒜`$. The crossover to monomer dominated behavior ($`t_2/t_m(1)\frac{1}{2}`$) is indicated at the highest concentration whereas at low concentrations this ratio can be much smaller.

## V Aggregate-Fissioning

Fissioning of aggregates leads to exponential growth as the fission products provide seeds for the next round of aggregation. In this subsection, we consider two mechanistic models of the fission process associated with either proteolytic cleavage of aggregates or mechanically induced breakage. We acknowledge that other models are possible.

We assume that the fission time is small compared to the aggregation time, and thus work in the limit of ‘instantaneous fission’ in which breakage of an aggregate happens much more rapidly than aggregation. This implies a narrow distribution of fission sizes peaked, say, at aggregate size $`𝒜`$, and in this limit our results are expected to be independent of the width of the fission distribution.

We consider two extreme limiting models of fission: (i)Mechanical. In this case, once the aggregate reaches fission size $`𝒜`$, it splits into two fragments of equal size $`𝒜/2`$. This should approximately describe the situation in which aggregate size is limited by nerve cell curvature (e.g., aggregation favors flat planar or linear structures, but the curvature of the neuron tends favors curved structures). (ii) Physiological. In this case, the aggregate can break into all smaller lengths at the fission scale. This mimics the outcome of protease attack for which there is no obvious preferred site for breakage.

Taking a fixed background concentration of monomers, which should be reasonable for at least short times in the disease. The kinetic equation for the time evolution of aggregates with size $`n`$ (measuring the number of dimers present) and concentration $`[a_n]`$ is, for $`n<𝒜`$,

$$\frac{d[a_n]}{dt}=p_{n1}[a_{n1}]p_n[a_n]+p_{f,n}[a_𝒜]$$
(12)
and, for $`n=𝒜`$,

$$\frac{d[a_𝒜]}{dt}=p_{𝒜1}[a_{𝒜1}]p_{f,0}[a_𝒜].$$
(13)
Here, for 1d aggregation $`p_n=p_0`$, while for the 2d aggregation specified by our critical coordination three rules discussed in Sec. III, $`p_n=p+np^{}`$. For mechanical fission, $`p_{f,n}=2p_{f0}\delta _{n,𝒜/2}`$, while for the physiological fission, $`p_{f,n}=2p_{f0}/(𝒜1)`$. The instantaneous fission assumption requires $`p_{f,0}>>p,p^{}`$.

We can identify the doubling time from Eqns. 12 and 13 by the following procedure: (i) Laplace transform the set of coupled equations to obtain a matrix equation in transform space; (ii) identify the largest positive eigenvalue of the Laplace matrix. In all cases, we find but one positive eigenvalue. We have systematically varied the fission size $`𝒜`$ and studied the dependence of the exponential growth rate upon fission time. For fibrils, the mean time to aggregate to size $`𝒜`$ is $`t_m𝒜/p`$, while for the 2d aggregates, the mean aggregation time goes as $`t_m\mathrm{ln}(𝒜)/p^{}`$. In the one dimensional case, we find that for large $`𝒜`$, the doubling time $`t_2`$ tends to $`0.5(0.43)t_m`$ for mechanical(physiological) fission. Hence, there is but a factor of two difference between the aggregation time and the doubling time. Since the numerical difference between mechanical and physiological fission is not substantial, we have examined only the mechanical fission model for the 2d aggregate. In this case, we find that the largest eigenvalue of the Laplace matrix goes as $`0.4/p`$ independent of $`𝒜`$, while the aggregation time scales as $`\mathrm{ln}(𝒜)/p^{}`$. Hence, for sufficiently large $`𝒜`$ it is possible to make $`t_2/t_m<<1`$.

These results have the attractive feature of linking the aggregation time, which we associate with the lag phase, to the doubling time in fission. However, we acknowledge that other processes may properly describe fission. In particular, we cannot rule out continuous fissioning of fragments off of large aggregates which may lead to a very different result provided the fission rate is comparable to growth rates.

## VI Dose-Incubation curve

In this section, we will look at the total incubation time and how it varies as a function of the inoculated dose using the aggregation-time distributions derived in the previous sections, and explore the extent to which it provides a quantitative description of experimental dose incubation curves. As discussed before, a key advantage of the 2D growth models is that, with a suitable fissioning scenario, they lead to lag times much larger than doubling times. This is difficult to accomplish with the 1D growth models. However, in this section we will assume that the doubling time is an independent free parameter. This allows us to fit the experimental dose incubation curves by both 1D and 2D models. The constraints on relative values of lag times and doubling times will be brought up in our discussions in the next section.

We now proceed to calculate the incubation time as a function of inoculated dose within our model for both 1d and 2d aggregation models and compare with laboratory data. Consider first the DIC of the 263K hamster scrapie strain. Kimberlin et al  have determined the DIC along with independent measurements of the doubling time $`t_2`$ and the final infectivity for this strain. The doubling time $`t_2`$ can also be inferred from the DIC in the region where it shows a logarithmic dose-dependence. Besides this experimental data, we need to know the clearance ratio $`r_i`$ which gives the percentage of the number of infectious seeds in a given inoculum which are removed by rapid initial clearance. Let us illustrate this for the case of 1 $`\mathrm{L}D_{50}`$ unit : in our model this corresponds to having a 50$`\%`$ probability of attaching a single infectious seed. Correspondingly the initial inoculum has to contain, on average, $`\frac{50}{100r_i}`$ seeds (for $`r_i>50\%`$, typical values from experiment are $`r_i99\%`$ ). We treat $`r_i`$ as a parameter in our model to be determined by fitting the DIC. Using this, in conjunction with our results for the aggregation time distributions, we can generate theoretical DICs using the method outlined in Section II.

For a fixed aggregation size $`𝒜`$, probability of attachment $`p`$ and infectivity ratio $`r_i`$, we then calculate the mean-square deviation ($`S^2`$) (normalized by the experimental error estimates for each data point) between the theoretical and experimental DICs. Minimizing $`S^2`$ gives us the optimal parameters, $`p`$ and $`𝒜`$, for the particular strain for a given value of $`r_i`$ within our model.

We consider first the DIC for the 263K hamster scrapie strain. In carrying out the fitting, we have ignored the results at the highest doses since in this limit we are approaching the saturation of the incubation time. For our 2d growth model, we get a good fit to the experimental data as indicated in Fig. 4a. The optimal parameters in this case are $`𝒜=16`$ (with $`𝒜_i=10`$), $`p=0.025day^1`$ and $`r_i=0`$ (which seems unphysical) for which $`S^20.07`$. For a more realistic value of $`r_i88\%`$ we get $`𝒜=16`$, $`p=0.16day^1`$ with $`S^21.14`$. From the fitting, one can see that the clearance ratio $`r_i`$ cannot be much greater than $`r_i88\%`$ due to the constraints imposed by the experimental incubation and doubling times. It should be noted furthermore, that the above procedure does not uniquely determine these parameters since comparably good fits are obtained for higher values of $`𝒜`$ by correspondingly adjusting $`p`$. For the 1-d growth model, we also get a good fit with $`S^20.35`$ for $`𝒜=8`$ (with $`𝒜_i=4`$) and $`p=0.11day^1`$.

One of the unusual features of the 263K scrapie strain in hamsters is that at high doses the lag-time is negligibly small. This feature is clearly seen in the time-course measurements  of infectivity and also accords with the theoretical best-fit results described above. In order to test our procedure for a more representative strain we have also obtained a fit for the experimental DIC for the ME7 strain in C57Bl mice . The results obtained by using our 2d growth model are shown in Fig. 4b. In this case, the theoretical fit is not as good as that obtained for the 263K strain; the optimal parameters correspond to $`𝒜=140,p=0.033day^1`$ for $`r_i=99`$ which gives $`S^23.75`$. A comparably good fit was obtained by using the 1d growth model with the optimal parameters $`𝒜=40,p=0.23day^1`$ for $`r_i=99`$ which gave $`S^25.02`$. However in contrast with the 263K strain, the theoretical results for the ME7 strain give rise to a significant lag time of $`50`$ days for the highest dose inoculated. This is a testable prediction for time-course measurements of infectivity for the Me7 strain in C57Bl mice.

Finally, we have used the above procedure to analyze the Sc237 scrapie strain in hamsters. In many respects, this strain is similar to the 263K scrapie strain in hamsters, however an analysis of the respective DICs reveals some significant differences. Based on the theoretical best-fit formula for the DIC , we can estimate a doubling time $`t_2`$ for the Sc237 strain to be $`2.1`$ days whereas the doubling time for the 263K strain is $`3.9`$ days. This difference is also reflected in our theoretical best fit DIC curve for this strain which, for $`r_i=99`$, is given by the parameters $`𝒜=80`$ $`p=0.052day^1`$. For our 1d model the corresponding values are $`𝒜=62,p=0.35day^1`$. Note that the experimental error estimates were not available for dose-incubation data for this strain, so we assumed them to be the same as that for the 263K strain in determining the best fit. Using the above fits, we see that the theoretical prediction for the mean lag-time for an inoculated dose of 1 $`\mathrm{L}D_{50}`$ unit is $`50`$ days. This value seems to be in good agreement with the experimental results for low dose inoculations for this strain .

The key results from our fitting (using the 2d growth model) are summarized in Table I. Based on the above results, we conclude that while our model accounts for the features of the DIC and gives a good fit to the experimental data, the latter cannot be used to distinguish between the 1d and 2d growth morphologies or to ascertain the model parameters conclusively. However despite the ambiguity in the model parameters, there are some robust predictions we can make after fitting the DICs. We find that the experimental DIC, in conjunction with our model, can be used to make predictions for the time-course of infectivity: in particular we can predict the lag-time as a function of dose. While the duration of the lag-time that we calculate does depend on the clearance ratio $`r_i`$, we note that the trend is that increasing $`r_i`$ reduces the lag-time. Thus by measuring the lag-time at high doses we can determine the parameter $`r_i`$ in our model, which then yields testable predictions for the lag-time at low doses. In the case of the the Sc237 scrapie strain in hamsters, our calculated lag-time at low dose agrees well with experimental results for the same. For the other strains, the predictions for the lag-time at low dose are a key testable prediction of our model.

## VII Connection to Epidemiological Data

We have found no dose incubation time data available for large mammals in the literature, so to gain insight we have analyzed epidemiological data for BSE in cattle, a tabulation of incubation times from iatrogenic CJD associated with dura mater transplants, and the cluster of five victims of vCJD from the village of Queniborough in the United Kingdom. The goal is to produce approximate estimates of aggregate size and growth time scales for comparison to the small mammal data, and, potentially, to guide future epidemiological and public health studies.

In brief, our assumptions and methods are as follows:
1) Model Distributions We apply only the 2d model in the dilute dose limit (suitable for digested prions); comparable quality fits can be obtained from the 1d model, but the estimated $`t_2/t_m`$ ratio consistently and strongly violates our mechanistic model result from Sec. V, while the bound for the 2d case is satisfied.
2) Number of Doublings. For mice and hamsters, with 1g brains, 30 doublings to incubation is typical. Given that the mean cattle brain is 500g, and the mean human brain 1500g, we take $`n_2=40`$ doublings from infection to incubation for cattle and humans.
3) BSE fits. Ref.  provides a candidate incubation time distribution which best fits the epidemic time course, and yields a mean incubation time $`<t>`$=5 years and standard deviation of 1.3 years. The width fixes $`p^{}`$ uniquely for given $`n,l=1+p^{}/p`$, with $`p^{}`$ weakly dependent upon $`n`$. We readily calculate $`t_m(n,l,p^{})`$$``$ $`(\psi (n+l)\psi (l))/p^{}`$, and we take the difference $`<t>t_m`$ to be $`n_2t_2`$, the length of the exponential growth phase. With a minimum $`t_2`$=5 days typical for hamsters and mice (see Table II), the maximum $`n`$ value can be be found at given $`l`$. For small aggregates, $`n=6`$ is taken to be a plausible minimum value which bounds $`t_2`$ above as 15 days. The results for $`l=2,10`$ are summarized in Table II. Importantly, we find the maximum aggregate size at fission to be of order 80-100 monomers, approximately independent of $`l`$ for $`2l10`$.
4) Fits to iatrogenic CJD data. The data of Ref.  sharply constrain the size of the aggregate at fission, since the ratio of standard deviation ($`3`$ years) to mean ($`5.8`$ years) rules out $`n>4`$ for all $`l=1+p^{}/p`$ values. To fit the data we perform the following procedure: a) We produce a parameter free estimate of $`t_m`$ by multiplying $`\sigma `$ by the ratio of $`t_m/\sigma `$ for a given choice of $`n,l`$; this subtracted from $`<t>`$ yields the estimated doubling phase period (and $`t_2`$, assuming $`n_2=40`$ as above). (b) Using the estimated value of $`n_2t_2`$, we produce a one parameter fit in $`p^{}`$ to the cumulative distribution $`F^{(1)}(tn_2t_2)`$ for the given $`n,l`$ choice. We find for $`n=3`$, we can consistently fit for all $`l`$, while for $`n=4`$, comparable quality fits are obtained for $`l=2,3`$, but for $`l4`$ $`n=4`$ is no longer a viable choice. The resulting robust estimates of incubation parameters are $`t_2=5.4\pm 0.4`$ years, $`\sigma =3.2\pm 0.2`$ years, and $`t_2=7\pm 1`$ days. The attachment rates $`p^{}`$ vary systematically with $`l`$, but are at most 0.22$`\pm 0.042`$/year. (Error bars reflect 95% confidence intervals.)
5) Fits to Queniborough vCJD incubation time data. For vCJD, we have extracted an estimated onset time distribution from the Queniborough cluster, taken as incubation times, from the available scientific and journalistic literature. We assume that this is a single dose event, in the low dose limit. We determine the estimated mean $`<t>`$ from the five data points to be Feb., 1998. Using non-linear leas squares analysis, we fit the cluster data to the cumulative beta distribution $`F^{(1)}(t<t>+t_m(n,l,p^{}))`$, placing points at half steps of probability increment. We also estimated the upper and lower bounds to $`p^{}`$ corresponding to 95% confidence, which yields a corresponding range for $`t_m,\sigma `$. Unlike the iatrogenic CJD and BSE data, the overall incubation time is unknown, which prohibits an estimate of $`t_2`$. We thus take a reasonable estimate for the maximum doubling time at $`t_2^{max}30`$ days. Using the minimum estimate of 9.0 years for the incubation time from epidemiological studies, we thus estimate a minimum lag time of $`t_m^{min}=t_{inc}^{min}n_2t_2^{max}=5.7`$ years. The minimum aggregate size at fission which can exceed this $`t_m^{min}`$ value at the upper 95% confidence limit is $`n=4,l=2`$ (16 monomers). We can take the upper aggregation size limit from BSE as a reasonable bound ($`n=20,l=2`$ or $`n=10,l=10`$ giving 80 monomers). We summarize our results in Table III, for two values of $`l`$. We stress that our results here are at best a crude guide to expected model fits since: (1) there are only five data points to fit to and several model parameters, (2) the infection is likely to be characterized by slow heterologous protein attachment for short times and more rapid homologous attachment for long times, while we assume a single effective attachment probability, and (3) there is no guarantee that the postulated single event was in the low dose limit.

Given the simplicity of our model and the assumptions made for fitting, we shall emphasize the most robust results of our analysis. Our results are summarized in Tables II (BSE), III (vCJD). We find several robust features:(i) For BSE, iatrogenic CJD, and vCJD, we find that a wide range of $`l`$ produce comparable fits, corresponding to a different starting seed size. We have thus quoted values from $`l=2`$ and $`l=10`$, the former representing nearly minimal seed sizes, the latter maximal plausible ones. (ii) In all three cases, the standard deviation of the incubation time distribution is large (1.2 years for BSE, vCJD, and 3.2 years for iatrogenic CJD), supporting the assumption of small dose. (ii) The attachment rates are all comparable (BSE: 0.6/yr; iatrogenic CJD: $``$ 0.2/year; vCJD: $`0.66`$/yr) and are significantly smaller than for the small mammal analysis of the previous section ($``$ 11/year from Table I). (iii) For BSE and CJD, $`t_2`$ is comparable to values for small mammals (eg., $``$15 days for BSE, $`8`$ days for iatrogenic CJD). For all strains considered here with the exception of 263K for hamsters, ratios of $`t_2/t_m`$ are small. We shall discuss the large $`t_2`$ values for vCJD below. (iv) Aggregate sizes are comparable for large and small mammals ($`80100`$ monomers for BSE, for iatrogenic CJD with $`l=2,n=3`$ we get 16 monomers, and for vCJD apparently tens of monomers also), compared with 16-360 for the small mammals of Table I. Hence in each case, given a 2 nm characteristic protein size, the length scale is of order tens of nanometers. (If the aggregating entities are in fact oligomers as suggested by Ref. , we only double the length scale.)

Hence, our analysis strongly suggests interspecies variation in incubation times is dominated by $`t_m`$, given given similar doubling times and aggregate sizes. In our model, variation in $`t_m`$ is governed by $`p^{}`$, which then suggests a significantly higher homeostatic monomer concentration in small mammals compared to large ones, a point we discuss further below.

Our fitting of the Queniborough data also suggests that long ($`20`$ years) incubation times for vCJD are implausible. The maximum $`t_m`$ value found (at the upper 95% confidence limit) for $`n=l=10`$ is 9.3 years; even with the extreme $`t_2=30`$ day value, this gives a 12 year incubation time. We note that the most probable (best fit) $`p^{}`$ values yield larger and implausible $`t_2`$ values for the minimum epidemiological estimate of the incubation time. To the extent our model is applicable, this implies 9 years likely overestimates the incubation period. In consequence, in conjunction with epidemiological analyses, our work suggest that the number of vCJD infections is likely to number in the hundreds.

## VIII Discussion

Recall that our basic hypothesis is that incubation times are controlled by prion-aggregation around infectious external seeds on the neuronal surface. Furthermore, in our calculations, the distribution of aggregation times arises from the stochastic growth process from seeds of a given initial size. That only a narrow range of seed sizes is relavant here maybe motivated by size-sensitivity of (a) the blood-brain barrier (b) attachment probability and (c) transportability of the seeds. The lag phase corresponds to growth from initial seeds to a characteristic fissioning dimension $`𝒜`$, after which one gets a multiplication in seeding-centers and an exponential growth in infectivity.

That there is a long lag time despite external seeding by intracerebral inoculation , and a doubling time which is typically significantly shorter , both of which become sharply defined at high doses, seems to be a general feature of the prion diseases. Our 2d compact aggregate model, with the assumptions of the preceding paragraph, explains these facts. In particular, 2d compact aggregation generates a broad distribution of aggregation times for a single seed with a well defined sharp onset time ($`t_0`$) and mean aggregation time ($`t_m(1)`$). With increasing number of seeds ($`D_i`$), the distribution of times for the first seed to reach the typical fissioning size $`𝒜`$ will narrow. Correspondingly the lag time, determined by the first fissioning event, will become sharply defined and concentrate at the onset time, which only weakly depends upon $`D_i`$.

The doubling time ($`t_2`$) is defined by fissioning and subsequent growth from size $`𝒜/2`$ to $`𝒜`$. Since growth from different seeds is independent, self-averaging (‘law of large numbers’) gives a sharply defined $`t_2`$. Thus at high doses both the lag time and doubling time are sharply defined which accounts for one of the most striking features of prion diseases: the reproducibility of incubation times at high doses. Indeed, we can explain several features of the dose-incubation curve. Notably, above a saturation dose $`D_s`$, the incubation time does not decrease, while for $`D<D_s`$, the incubation time varies as $`\mathrm{log}(D)`$, showing deviations from the log only below a much smaller value $`D_{min}`$ . The total incubation time is the sum of the lag time and $`n_dt_2`$, where $`n_d`$ is the number of doubling steps. Assuming that the onset of clinical symptoms is related to the damage of a fixed number of neurons , the logarithmic dose dependence of the number of doubling steps follows from the fact that number of seeds grows exponentially in the fissioning stage. In the range $`D_{min}<D<D_s`$, the lag time does not change appreciably with dose, thereby giving rise to the logarithmic dose-dependence of the incubation time in this range. At low doses ($`D<D_{min}`$), the lag-time increases towards $`t_m(1)`$ (the mean aggregation time for a single seed) giving rise to a broad distribution of incubation times and the deviation from the logarithmic behavior in the DIC which is observed experimentally.

Furthermore, we note that in our mechanistic fissioning model, the doubling time ($`t_2`$) is bounded above by the time to grow from size $`𝒜/2`$ to $`𝒜`$. If the fission produces jagged fragments, these can be effectively filled by monomers which will accelerate the subsequent growth process. This is only possible for 2d compact aggregates and not for 1d fibrils, for which the exposed ends will always be limited to dimer growth. This possibility may account for the effective ‘1/c’ dependence in the incubation time observed for transgenic hamsters with multiple copies of the prion gene , noting that for hamsters the doubling phase appears to dominate incubation .

A key difference between the 1d and 2d morphology (shown in section IV) is that, within our model, in case of the latter (i) the lag-time can be an order of magnitude larger than the doubling time. If the total time in doubling-steps becomes large compared to the lag time, the overall distribution will be relatively narrow. Thus, only in the case of 2d growth can one get (ii) a wide distribution for the overall incubation time, with a width comparable to the mean. Thus assuming (i)-(ii) to result entirely from the growth processes discussed here, strongly points to a 2d (or 3d) morphology as controlling the incubation times.

The early growth morphology clearly deserves further experimental attention . Typically, the in vitro morphology of prion aggregates has been found to be fibrillar . Frequently large fibrillar aggregates are also observed post mortem in brain tissues. The morphology and size scale for aggregates that cause neuronal death and infection is not known. One could argue that the reason why in vitro aggregates are not infectious is because they do not have the proper morphology. We speculate that the attachment to lipid membranes could make a vital difference to the aggregation process, which is missing in in vitro experiments. It would be very interesting to carry out the in vitro studies of prion aggregation in presence of lipid membranes.

An important byproduct of our analysis is the ability to predict the time course of infectivity from the DIC. Provided we take the doubling time ($`t_d`$) as an independent experimental parameter, such predictions are very robust and do not rely on many details of the aggregation-fissioning model, including initial growth morphologies. Such predictions are particularly significant since experiments which measure the time-course of infectivity, and hence determine the lag phase, are extremely expensive and time consuming. This dose dependence of the lag phase may well be a significant factor in assessing the risk of infection.

Our results indicate that we can infer the (average) time-course of infectivity using the information supplied by the experimental DIC. Thus, based on our fits to the experimental DIC, we have made testable predictions for the time-course of infectivity, and in particular the lag time, as a function of dose. These predictions are in good agreement with the existing experimental results and their further experimental validation would prove very useful.

A factor which significantly affects the lag time is the probablity of dimer attachment $`p`$; lowering $`p`$ increases the lag time. This is relevant in understanding the species barrier effect in which there is a reduction of incubation times with multiple passages in inter-species infection . During first passage, the attachment of dimers is initially non-homologous but as the seed size increases it should change to homologous attachment. Since the non-homologous attachment probability should be smaller , the lag phase should be longer for first passage as compared to subsequent passages. Thus, in our picture, most of the difference in incubation times should come from the lag phase and the exponential growth phase should be similar between first and second passages. This has been observed experimentally for hamster scrapie passaged in mice . We note that the estimated dimer attachment rates for mice and hamsters range from 9-60 per year (c.f. Table I), while for humans and cattle they are maximally 0.6 per year.

Furthermore, we observe that the estimated aggregate sizes are comparable between large animals (humans and cattle) and small animals, all in the ballpark of tens of nm, which is precisely the estimated size of the lipid rafts on which prions are hypothesized to rest. The comparable sizes of aggregates at fission we find in our model between small and large mammals strongly suggests that the primary determinant of lag time is the dimer attachment rate which is regulated primarily by the concentration (for a given strain). This leads us to speculate that the concentration of normal prion proteins must vary dramatically between large and small species, which naturally leads us to envision a link to metabolic rate. Such a link is plausible if prion proteins play a functional role in relieving oxidative stress as has been proposed elsewhere. The hypothesis of enhanced homeostatic PrPc concentration in small animals relative to large ones is testable by direct examination of the brains of uninfected animals.

Finally, as noted in the preceding section, our model analysis yields 12 years as an extreme upper bound on the vCJD incubation time. We actually get reasonable doubling times only when we assume a total incubation below 9 years, which is the lowest epidemiological estimates. The important suggestion from this analysis is that the vCJD cumulative infection toll is likely to be several hundred, vs. several hundred thousand.

Acknowledgements. We acknowledge useful discussions with F. Cohen. We thank D.D. Cox for a critical reading and discussion of our fit to the Queniborough data. R.V.K. and D.L.C. acknowledge support from the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Research. A.S. is supported by Sandia which is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin company, for the United States Department of Energy under Contract No. DE-AC04-94AL85000. R.R.P.S. and D.L.C. have benefitted from discussions at workshops of the Institute for Complex Adaptive Matter. We are grateful for a grant of supercomputer time from the Lawrence Livermore National Laboratory.


text:
# First-principles study of band offsets in ferromagnetic semiconductor heterojunctions

## I Introduction

Spin “injection” refers to the induced transport of spin-polarized electrons or holes from one material into another. To date, two basic approaches to injection have been explored: photoinduced generation of spin-polarized carriers in one material, which then diffuse across an interface into a second material;  and direct electrical injection via application of a bias voltage across an interface. In both approaches, the heterojunction band offset represents an important quantity which may place fundamental constraints on the expected efficiency of the injection process. In this work we report a first-principles numerical study of band offsets for realistic heterojunctions between magnetic and nonmagnetic semiconductors.

We will refer to the magnetic and nonmagnetic sides of the heterojunction as the “source” and “sink,” respectively. Our choice of the source and sink materials are guided by a few simple principles. First, since spin lifetimes are generally longer for electrons and since electrons generally have higher mobilities and smaller spin-orbit interaction than holes, we focus on source materials which can supply spin-polarized electrons rather than holes. Second, we restrict our attention to semiconducting source materials so as to avoid the well-known “conductivity mismatch” problem that plagues direct injection from a metal into a semiconductor. Third, because lattice mismatch leads to dislocations (which act as spin scatterers), we consider sink materials with a reasonably close lattice match to the source material. A particularly interesting candidate source material is the chromium chalcogenide spinel CdCr2Se4, an n-type ferromagnetic semiconductor with a relatively high Curie temperature of 130 K. For sink materials, we consider Si and GaAs, both of which can be made n-type with appropriate doping and are well lattice matched to CdCr2Se4.

In this work we use density functional theory to calculate band alignments and offsets between CdCr2Se4 and Si or GaAs. We take the view that although the band alignment may ultimately place constraints on the injection efficiency, one may be able to “tailor” the interface during the growth process, thereby improving alignments that are otherwise unfavorable for injection. It is this variability in the alignments and offsets, with respect to various features of the interface, that we explore numerically in this work. Of the three possible types of alignment (see Fig. 1), the two most favorable for spin-injection of electrons are Type I or IIA. Less favorable is Type IIB alignment, which leads to Zener breakdown  (valence electrons of the source material tunnel into the conduction bands of the sink). Type IIB alignment also gives rise, through band bending, to confinement of holes in the source material in the layers adjacent to the interface, which would then act as strong spin-scattering centers and presumably degrade spin transport across the interface.

## II Theory of band offsets

In this work we evaluate the band offsets using methods pioneered by Van de Walle and Martin, here extended to include spin splitting. For completeness we review the essentials of that approach. We take the sink material to be on the left (L) side of the heterostructure, and the source material to be on the right (R). The valence band offset (VBO) is the relative alignment of the L and R valence band maxima (VBM) far from the interface, and is conveniently written as

$$E_{\mathrm{VBO}}=\mathrm{\Delta }E_{\mathrm{band}}+\mathrm{\Delta }V_{\mathrm{avg}}.$$
(1)
Here $`\mathrm{\Delta }E_{\mathrm{band}}`$ and $`\mathrm{\Delta }V_{\mathrm{avg}}`$ are generally referred to as the band lineup and potential lineup, respectively. $`\mathrm{\Delta }E_{\mathrm{band}}`$ is determined entirely by the electronic structure of the constituent bulk materials, while $`\mathrm{\Delta }V_{\mathrm{avg}}`$ accounts for the effects of the interface between those materials. For simplicity we do not show spin indices in Eq. 1; each spin channel is treated independently.

Although the band lineup, $`\mathrm{\Delta }E_{\mathrm{band}}`$, is a purely bulk quantity, its evaluation in general depends on the nature of the interface. For example, lattice mismatch between the substrate and overlayer will introduce tetragonal strain into an overlayer that has a nominally cubic lattice. This strain will in general lead to level shifts and band splittings that must be correctly included when computing the bulk band lineup. For the purposes of this work we consider the nonmagnetic sink materials to be unstrained substrates on which overlayers of CdCr2Se4 are grown epitaxially. Thus the overlayers experience an in-plane tensile strain (1.3% and 3.3% for Si and GaAs, respectively) that results in a contraction of the spacing between the overlayers. To fully account for these strain effects we consider bulk unit cells for CdCr2Se4 with a tetragonal distortion $`a_1=a_2=2a_{\mathrm{substrate}}a_3`$.

We also include fully the effects of the spin-orbit (SO) interaction on the band structure, by performing fully relativistic calculations for the bulk phases. With strain and spin-orbit effects thus included, the band lineup takes the form

$`\mathrm{\Delta }E_{\mathrm{band}}`$ $`=`$ $`(E_{\mathrm{VBM}}^\mathrm{L}V_{\mathrm{avg}}^\mathrm{L}+\mathrm{\Delta }E_{\mathrm{SO}}^\mathrm{L})`$ (2)
$``$ $`(E_{\mathrm{VBM}}^\mathrm{R}V_{\mathrm{avg}}^\mathrm{R}+\mathrm{\Delta }E_{\mathrm{SO}}^\mathrm{R}).`$ (3)
Here, $`E_{\mathrm{VBM}}`$ is the position of the valence-band maximum, and $`V_{\mathrm{avg}}`$ is the cell-averaged Kohn-Sham potential (which provides a convenient reference energy). $`\mathrm{\Delta }E_{\mathrm{SO}}`$ is the shift in the VBM associated with spin-orbit interaction. Figure 2 shows an an energy level-diagram including all the various terms defined above.

The second term in Eq. 1, the potential lineup $`\mathrm{\Delta }V_{\mathrm{avg}}`$, is determined by the electronic structure of a specific interface. We use periodic supercells to represent the heterojunction, with supercell periods sufficiently large to ensure that the interfaces do not interact. The potential lineup is the difference in the locally averaged potential on the two sides of the interface, far from the interface itself. To compute these local averages we use the slab-averaging technique discussed in Ref.. One first computes the planar average, $`\overline{V}(z)`$, of the total potential, $`V(𝐫)`$, by averaging across planes parallel to the interface. $`\overline{V}(z)`$ contains atomic-scale variations whose local periodicity on the two sides of the interface may in general may be different—even far from the interface—because of lattice mismatch. The correct local averages can be obtained by performing a “slab average” over a sliding window. This slab averaging is applied twice, in order to allow for different L and R windows:

$$\overline{\overline{V}}(z)\frac{1}{\lambda _L\lambda _R}_{z\lambda _L/2}^{z+\lambda _L/2}_{z^{\prime \prime }\lambda _R/2}^{z^{\prime \prime }+\lambda _R/2}\overline{V}(z^{})𝑑z^{}𝑑z^{\prime \prime }.$$
(4)
Here $`\lambda _L`$ and $`\lambda _R`$ denote the periodic repeat distances in the direction perpendicular to the interface for the L and R materials . Far from the interface, the slab-averaged potential, $`\overline{\overline{V}}(z)`$, is constant. The potential lineup is then given by

$$\mathrm{\Delta }V_{\mathrm{avg}}=\overline{\overline{V}}(z_L)\overline{\overline{V}}(z_R),$$
(5)
where $`z_L`$ and $`z_R`$ denote the midpoints of the L and R slabs. Finally, having computed the VBO directly from Eq. 1, the conduction band offset (CBO) is simply

$$E_{\mathrm{CBO}}=E_{\mathrm{VBO}}+E_{\mathrm{gap}}^LE_{\mathrm{gap}}^R,$$
(6)
where $`E_{\mathrm{gap}}^{L(R)}`$ is the band gap for L (R). For reasons discussed below (see Sec. V), we evaluate Eq. 6 using experimental band gaps rather than our theoretical values.

## III Interface Models

The crystal structure of bulk CdCr2Se4 is shown in Fig. 3. In this work we consider only (001) interfaces, for which tetrahedral bonding across the interface can be maintained in a particularly simple manner. In Fig. 3, different (001) basal planes can formed by truncating the spinel structure at different atomic planes. Only two chemically inequivalent planes can be so formed: one containing only Cd atoms, and one containing Se and Cr atoms in a 2:1 ratio. Such a representative plane containing Se and Cr atoms is indicated in Fig. 3. Within these planes, the atomic density in the Cd plane is only one fourth that of Si or GaAs (ignoring the lattice mismatch), which suggests that a 2$`\times `$2 reconstruction of the terminating Si layer would result for interfaces with the Cd plane. On the other hand, for interfaces with the mixed (Se,Cr) plane a simple unreconstructed interface is quite plausible, and so we limit ourselves to variants of this basic interface geometry.

For interfaces with Si, the atomic alignment at the interface is shown schematically in Fig. 4. For the case of Si we consider only atomically abrupt interfaces. At the interface boundary, every Si atom is tetrahedrally coordinated (to two Si and two Se atoms). The role of Cr atoms at the interface is not immediately obvious from these preliminary considerations, but will be analyzed in detail in Sect. VI.

For interfaces with GaAs we consider both abrupt and partially intermixed interfaces. In both cases we compute offsets for both Ga-terminated and As-terminated substrates. Fig. 5 shows an abrupt As-terminated substrate; Ga-terminated substrates are analogous.

Although abrupt interfaces are conceptually attractive, for a polar material like GaAs they may not be physically plausible. In particular, for interfaces between polar and non-polar materials, abrupt termination leads inevitably to the formation of long-range dipole fields. Such a dipole field is energetically unfavorable over large distances, and will typically lead to an atomic rearrangement near the interface in order to reduce such effects. For polar/nonpolar interfaces, a particularly simple way of eliminating the dipole field is to intermix the polar and nonpolar constituents at the interfacial layer in order to maintain local charge neutrality. We have followed this approach for interfaces with GaAs, by appropriately substituting Ga or As atoms for Se atoms at the interface in such a way as to maintain local neutrality. Fig. 6 shows an example of this substitution for an As-terminated interface.

## IV Computational Details

For both the interface and bulk calculations we use the plane-wave pseudopotential method with norm-conserving pseudopotentials, as implemented in the abinit code.  In both methods we use the local-spin-density approximation (LSDA) in the form of Perdew and Wang. We have carefully checked the accuracy of the pseudopotentials by comparing to results for bulk CdCr2Se4 obtained with the linearized augmented plane-wave (LAPW) method as implemented in the wien97 code. Table I summarizes the predictions of both methods for equilibrium lattice constant, bulk modulus, and magnetic moment, along with the available experimental values. The pseudopotential and LAPW lattice constants agree to better than 1%, the bulk moduli to better than 3%, and the magnetic moment essentially exactly. The agreement with experimental values is equally good with the exception of the bulk modulus (for which the spread in experimental values is very large). We have also computed the difference in total energy between the normal- and inverse-spinel forms of CdCr2Se4,

$$\mathrm{\Delta }E_s=E_{\mathrm{normal}}E_{\mathrm{inverse}},$$
(7)
and find that the spinel structure is energetically preferred (see Table I) within both PSPW and LAPW methods. Thus for the rest of this work we consider only the normal spinel phase.

The interface results we report here were obtained with supercells containing five atomic layers of Si or GaAs and five layers of CdCr2Se4. These thicknesses are sufficient to ensure that the interfaces do not interact, as confirmed by the convergence (discussed below) of the planar-averaged potentials to their bulk values within the first few layers of the interface. The coordinates of all the atoms were completely relaxed until all force components were smaller than $`10^3`$ Ha/bohr. By simultaneously optimizing the supercell dimension normal to the interface, we impose no artificial constraints on the interlayer spacing throughout the supercell. As described above, however, we do constrain the in-plane lattice constant to that of Si or GaAs. The resulting tensile strain experienced by the CdCr2Se4 overlayers gives rise to a decrease in the spacing between bulk-like CdCr2Se4 layers, consistent with total energy calculations for the bulk material with a tetragonal distortion. The cumulative effect of all these physical relaxations (relative to the ideal geometry) on the absolute band offsets is typically of the order of 0.1–0.3 eV, and so must be fully included for an accurate description of the alignments.

For all the pseudopotential calculations, a plane-wave cutoff of 70 Ry was used, which is sufficient to converge all band states relevant to the offset calculations to within 100 meV. A 4$`\times `$4 Monkhorst-Pack  sampling of the surface Brilloun zone was used for the supercell calculations, equivalent to 8$`\times `$8 sampling of the primitive Si or GaAs surface zone.

## V Bulk CdCr2Se4 Band Structure

The electronic structure of bulk CdCr2Se4 has been the subject of a number of experimental studies. Photoluminescence experiments give the fundamental gap —the energy difference between the highest p-like valence state and lowest s-like conduction state—to be 1.8 eV in the ferromagnetic phase. In contrast, optical studies reveal an absorption edge only 1.3 eV above the valence-band edge, suggesting the existence of d-like conduction states in the fundamental gap. This view is supported by transport measurements showing relatively low mobility, consistent with conduction taking place within narrow d-like bands.

Our results for the band structure of bulk CdCr2Se4 are show in Fig. 7. We find that the ground state is indeed a ferromagnetic semiconductor, with a magnetic moment of 6.0 $`\mu _B`$ per formula unit. Our results are in very good agreement with other recent LAPW calculations, but are slightly different from results obtained with the linearized-muffin-tin-orbital (LMTO) method, which gives semimetallic band structure in both spin channels and a slightly higher moment of 6.2 $`\mu _\mathrm{B}`$ per formula unit.

Figure 7 shows the density of states for bulk CdCr2Se4 decomposed into angular character and projected onto atomic sites. Each Cr site contributes a moment of $`3\mu _B`$ to the total moment, corresponding to four spin-up and one spin-down electron. Our predicted fundamental gaps, indicated in Fig. 7 by arrows, are found to be 2.42 and 3.05 eV for the majority and minority spin channels respectively; these results are in very good agreement with those reported in Ref..

For spin injection, the relevant energy gap is the band gap between valence and conductions states. Our results predict the band gap to be very small, of order 0.3 eV in both spin channels; this value is much smaller than the experimental optical absorption edge at 1.3 eV. Since there is experimental evidence suggesting that conduction takes place in the unoccupied d levels probed by optical absorption, we simply adopt the experimental value of 1.3 eV for the band gap of CdCr2Se4 appearing in Eq. 6. For similar reasons, we use experimental values of the band gaps of Si and GaAs.

## VI Calculated Band Offsets

We begin this section by illustrating the slab-averaging procedure used to determine the potential lineup across an interface. Figure 8 shows the planar-averaged and slab-averaged potentials for the case of CdCr2Se4/Si. The five layers of Si are in the middle portion of the figure, the atomic planes corresponding to minima in the planar-averaged potential; the three maxima in the CdCr2Se4 portion of the cell correspond to Cd planes. The slab-averaged potential is quite flat at the midpoints of the Si and CdCr2Se4 slabs, demonstrating that these slabs are sufficiently thick to converge the potential-lineup contribution to the band offset. Figure 8 also shows the planar-averaged potential for bulk Si. The rapid convergence of the planar-averaged potential in the supercell to its bulk value is readily apparent, and is similar to that found in previous studies.

Our calculated valence- and conduction-band offsets for the five interfaces described in Sec. III are summarized graphically in Fig. 9 and tabulated numerically in Table II. In Table II we show only offsets for the valence bands since only these were evaluated completely from first principles. For all five interfaces studied, we find that the CdCr2Se4 conduction-band edges for both spin channels sit well above the conduction-band edge of the substrate—independent of the substrate material (Si or GaAs), termination layer (Ga or As), and the atomic intermixing described in Sec. III. The band alignments are all Type IIB for the majority spin channel, and either IIA or weakly IIB for the minority channel. Thus, according to the criteria outlined in Sec. I, these band alignments are not favorable for efficient spin injection.

A major goal of this work is to identify specific interface features that strongly affect band offsets, and to evaluate whether such features can be externally controlled, for example, during the growth process. For several of the interfaces we have studied, we have identified one such feature: the concentration of Cr atoms at the interface. Our motivation for studying the effects of varying the interfacial Cr content comes from analyzing the environmental dependence of Cr removal energies. For example, the calculated formation energy for a Cr vacancy in bulk CdCr2Se4 is 3.27 eV, but for Cr atoms in the interface layer of a CdCr2Se4/Si interface the formation of a single vacancy is actually exothermic, with an energy gain of 1.60 eV/vacancy.

To explore the impact of interfacial Cr content on band offsets, we consider the limiting case of total depletion of Cr from the mixed (Se,Cr) layer at the interface. Since the Cr atoms do not participate directly in the tetrahedral bonding across the interface boundary, their removal leads to relatively minor rearrangements of the interface structure—which we nevertheless fully account for by re-optimizing the geometry. The energetics of removing this Cr layer vary somewhat from interface to interface, but even in the limiting case of a Cr-rich growth environment, the total depletion of Cr from the final interface layer is thermodynamically favorable for all five interfaces, as shown in Table III.

The band alignments and offsets for Cr-depleted interfaces are summarized in Fig. 10. Compared to the non-depleted interfaces, all of the offsets are reduced. The reduction is remarkably large for the interface with Si (about 1 eV), smaller for both abrupt and intermixed Ga-terminated GaAs interfaces (about 0.5 eV), and smallest for As-terminated GaAs interfaces (about 0.1 eV). In two cases (Si, and Ga-terminated GaAs) the reduction is actually sufficient to change the alignment type from IIB to IIA.

Finally, we emphasize that the question of whether the concentration of interfacial Cr can be externally controlled remains open. Although the energetics clearly favors Cr-depleted interfaces, kinetic barriers may be more important—and perhaps decisive—at the relatively low growth temperatures typical of molecular-beam epitaxy. Theoretical studies of complex interfaces are notoriously difficult, not the least because of uncertainties in the interface morphology. The present work clearly demonstrates the possible dangers of assuming an “idealized” interface structure, and should make clear that an understanding of the systematics of band offsets requires—at the very least—detailed microscopic calculations.

TABLE

## VII Summary

In summary we have applied first principles methods to determine the absolute valence band offsets for heterojunctions between CdCr2Se4 an n-type ferromagnetic semiconductor, and Si or GaAs. By using experimentally determined band gaps, we have also determined the type of band alignments in these heterojunctions, and find favorable alignment for electrical injection for Cr-depleted heterojunctions with Si and with Ga-terminated GaAs substrates. These findings suggest that CdCr2Se4 may be a promising candidate for the injection of spin-polarized electrons. Finally, we point out that although the type Type IIB alignments we have found are not favorable for spin injection, such alignments are ideal for spin-polarized resonant interband tunneling diodes in which CdCr2Se4 would function as a spin-selective quantum well.

## VIII Acknowledgements

One of the authors (J.M.S.) acknowledges the National Research Council for support during this work in the form of a postdoctoral fellowship. We thank Andre Petukhov, Steve Hellberg, Berry Jonker and Aubrey Hanbicki for many useful conversations during the course of this work. This work was funded in part by DARPA and ONR. Computational work was supported in part by a grant of HPC time from the DoD Major Shared Resource Center ASCWP.


text:
# Quenched Chiral Behavior of Hadrons with Overlap Fermions Presented by T. Draper at Lattice 2001, Berlin, Germany. This work is supported in part by the U.S. Department of Energy under grant numbers DE-FG05-84ER40154 and DE-FG02-95ER40907.

## 1 Simulation Details

Using a $`\beta =7.60`$ (tree-level tadpole-improved) Lüscher–Weisz  gauge action, we study the chiral properties of hadrons on a $`20^4`$ lattice with the overlap fermion . The massive Dirac operator  is defined so that the tree-level renormalization of mass and wavefunction is unity.

$`D(m_0)`$ $`=`$ $`(1{\displaystyle \frac{m_0a}{2\rho }})\rho D(\rho )+m_0a`$
$`=`$ $`(\rho +{\displaystyle \frac{m_0a}{2}})+(\rho {\displaystyle \frac{m_0a}{2}})\gamma _5ϵ(H)`$
$`ϵ(H)=H/\sqrt{H^2}`$, $`H=\gamma _5D_w`$, and $`D_w`$ is the usual Wilson fermion operator, except with a negative mass parameter $`\rho =1/2\kappa 4`$ in which $`\kappa _c<\kappa <0.25`$; we take $`\kappa =0.19`$ in our calculation which corresponds to $`\rho =1.368`$. See  for more details of the simulation.

## 2 Pion Mass and Chiral Logs

To search for chiral logs, we fit $`m_\pi ^2a^2`$ to

$`m_\pi ^2a^2=`$
$`Am_0a\{1\delta [\mathrm{ln}(Am_0a/\mathrm{\Lambda }_\chi ^2a^2)+1]\}+Bm_0^2a^2`$
The best fits which give stable values of A and B and with errors less than half of the fitted values of $`\delta `$ for a range of $`\mathrm{\Lambda }_\chi =0.6`$ GeV to 1.4 GeV are listed in the table.

We plot the fit with $`\mathrm{\Lambda }_\chi =0.8\mathrm{GeV}`$ as a solid line in Fig. 1. To check if there is indeed a quenched chiral log, we fit $`m_\pi ^2a^2`$ alternatively without the chiral log term and find

$`m_\pi ^2a^2=`$
$`1.88(10)m_0a2.3(11)m_0^2a^2+7.8(30)m_0^3a^3`$
with $`\chi ^2/NDF=0.46`$. A quadratic fit (shown in Fig. 1 as the dotted line) has $`\chi ^2/NDF=0.90`$. Thus the presence of a chiral log term is favored.

## 3 Quark Masses

To renormalize our estimates of the quark mass, we follow  to determine the renormalization constants $`Z_M^1=Z_S=Z_P`$ and thus $`m^{RGI}=Z_M(g_0)m_0(g_0)`$: calculate

$`Z_M(g_0)`$ $`=`$ $`U_M{\displaystyle \frac{1}{(r_0m_0)}}|_{(r_0m_\pi )^2=x_{\mathrm{ref}}}`$ (4)
using their tabulated values

$`U_M`$ $`=`$ $`\{\begin{array}{cc}0.181(6),\hfill & x_{\mathrm{ref}}=1.5736\hfill \\ 0.349(9),\hfill & x_{\mathrm{ref}}=3.0\hfill \\ 0.580(12),\hfill & x_{\mathrm{ref}}=5.0\hfill \end{array}`$ (8)
where $`U_M`$ is a universal factor which they estimated from published results of another regularization, namely, $`𝒪(a)`$-improved Wilson.

The value of the bare quark mass $`m_0a`$ that reproduces a given pseudoscalar mass (through $`x_{\mathrm{ref}}=(r_0m_\pi )^2`$, with $`r_0/a=4.05(2)`$ for our lattice) is obtained by interpolation of our data for the pion mass. From this we obtain $`Z_M^1=Z_S`$ at three values of quark mass (the lowest two quark masses are roughly half-strange and strange) as shown in Fig. 2.

The operators and action are $`𝒪(a)`$ improved. For the axial vector renormalization constant, $`Z_A`$, we find small $`(m_0a)^2`$ dependence as well . Thus we fit $`Z_M^1=Z_S`$ to a constant.

$$Z_M^1=Z_S=1.233(34)\chi ^2/NDF=0.86$$
(9)
From the $`(m_\pi a)^2`$ versus $`m_0a`$ figure, we obtain $`m_0a`$ at the physical pion mass. Using $`m^{RGI}=m_0/Z_S`$ and, from the 4-loop calculation with $`N_f=0`$ , $`Z_S^{\overline{MS}}=Z_S/0.72076`$, we obtain (very preliminary results) for $`(m_u+m_d)/2`$.

$`m^{RGI}`$ $`=`$ $`7.3(7)(10)\mathrm{MeV}`$ (10)
$`m^{\overline{MS}}(\mu =2\mathrm{GeV})`$ $`=`$ $`5.3(6)(7)\mathrm{MeV}`$ (11)
## 4 Chiral Loops in Quenched Propagators

Bardeen et al.  present evidence of a strong “quenched chiral loop” (QCL) effect, attributable to an $`\eta ^{}`$$`\pi `$ intermediate state, in the valence propagation of the scalar, isovector meson ($`a_0`$). The effect is manifested in a negative-metric contribution to the two-point correlation function which increases for decreasing quark mass. They were able to uncover the effect only after using their pole-shifting procedure, which mitigates the poor chiral properties of the Wilson fermions. We see this effect clearly with overlap fermions without the need for any such procedure.

Fig. 3 shows the two-point scalar-isovector local-local correlators for quark masses $`m_0a=0.0150`$, $`0.0164`$, $`0.0192`$, and $`0.0274`$; they are plotted on a linear scale to show the correlator becoming negative at early times ($`t/a12`$, the source is at $`t_0=3`$) and then asymptotically approaching zero from below.

We also see the effect for $`N^{}(\frac{1}{2}^{})`$. Fig. 4 shows the two-point $`N^{}`$ local-local correlator for quark masses $`m_0a=0.0150`$, $`0.0547`$, and $`0.1642`$. The correlation function becomes negative at successively earlier times for decreasing mass. The effect becomes very pronounced at the lighter masses. See  for more details.

## 5 Summary

We found evidence for chiral logs in the dependence of the pion mass on quark mass. (See  for more details.) We also presented, for the first time, our very preliminary results for the renormalized light quark mass and the effect of chiral loops on the $`a_0`$ and $`N^{}`$ propagators. In a companion proceeding , we study zero modes and their effect on the pion propagator, and show our results for the axial renormalization constant $`Z_A`$ and pion axial decay constant $`f_\pi `$, and the presence of chiral logs in the pseudoscalar decay constant, $`f_P`$. See  for complete details.