text:
# 1 Introduction

## 1 Introduction

Let $`(S(n),n_+)`$ be a simple symmetric random walk on $``$ with $`S(0)=0`$. Define

(1.1)
$$\xi (n,x):=\mathrm{\#}\{0kn:S(k)=x\},$$

which is referred to as the (site) local time of the random walk. For each $`n`$, consider

(1.2)
$$𝕍(n):=\{x:\xi (n,x)=\underset{y}{\mathrm{max}}\xi (n,y)\},$$

which is the set of sites which are the “most visited” by the random walk at step $`n`$. Following Erdős and Révész , an element of $`𝕍(n)`$ is called “favourite site”.

The study of $`𝕍(n)`$ was initiated by Erdős and Révész , and Bass and Griffin , and received much research interest from other mathematicians since. Somewhat surprisingly, many of the easy-to-formulate and innocent-looking open questions raised by Erdős and Révész in  and  remain unanswered so far. The present paper aims not only to have an overview upon known results in this field, but also to insist on unsolved problems in the hope that they will meet the interest of, and find solutions from, the reader.

To illustrate these solved or unsolved problems, we mention the following question: what is the probability that $`0𝕍(n)`$ for infinitely many $`n`$?

Since the random walk is symmetric, one would be tempted to think that this probability would be 1. However, the correct answer is “0”. In fact, Bass and Griffin  proved the following result:

(1.3)
$$\underset{n\mathrm{}}{lim}\underset{x𝕍(n)}{inf}|x|=+\mathrm{},\text{a.s.}$$

In words, the process of favourite sites is transient. Actually, Bass and Griffin  showed that, for any $`\epsilon >0`$, the distance of the set of favourite sites from the origin goes faster than $`\sqrt{n}/(\mathrm{log}n)^{11+\epsilon }`$ but slower than $`\sqrt{n}/(\mathrm{log}n)^{1\epsilon }`$, almost surely. We will discuss this result in more detail in Section 2.

Let us mention another innocent-looking question. It is trivial that $`\{\mathrm{\#}𝕍(n)=1,\text{i.o.}\}=1`$ (where “i.o.” stands for infinitely often). A little more thinking reveals that $`\{\mathrm{\#}𝕍(n)=2,\text{i.o.}\}=1`$. Erdős and Révész  asked:

$$\{\mathrm{\#}𝕍(n)3,\text{i.o.}\}=\mathrm{?}$$

This problem is still open. However, the following was recently proved by Tóth : $`\{\mathrm{\#}𝕍(n)4,\text{i.o.}\}=0`$. For more details, see Subsection 3.2.

The rest of the paper is splitted into four sections according to the natures of the problems invoked. Section 2 concerns the problems of how far and how close the favourite sites can be to the origin. For each of these two problems, we have some useful but incomplete information. In Section 3, we mention ten questions raised by Erdős and Révész in  and , and quoted as Open problems 1–10 in the book of Révész  (pages 130–131). Only a few of these questions have found solutions. Section 4 is devoted to some related problems. More precisely, we will discuss problems for rarely visited sites, favourite edges, and the location of favourite sites. Finally in Section 5, we briefly describe some known results for favourite sites of other stochastic processes such as Brownian motion and general Lévy processes.

## 2 Large and small values of the favourite sites

Throughout this section, we pick an arbitrary element $`V(n)𝕍(n)`$. According to (1.3), $`|V(n)|`$ goes to $`\mathrm{}`$ almost surely (when $`n\mathrm{}`$). The question is to determine the rate with which $`|V(n)|`$ goes to infinity. The answers are different for lower and upper limits.

### 2.1 Escape rates of favourite sites

As far as the lower limits are concerned, the best possible result available so far is due to Bass and Griffin .

Theorem A (). With probability one,

(2.1)
$$\underset{n\mathrm{}}{lim\; inf}\frac{|V(n)|}{n^{1/2}(\mathrm{log}n)^\gamma }=\{\begin{array}{cc}0\hfill & \text{if }\gamma <1\text{,}\hfill \\ \mathrm{}\hfill & \text{if }\gamma >11\text{.}\hfill \end{array}$$

Throughout the paper, when we state a limit result for $`V(n)`$ as in (2.1), it is to be understood that the convergence holds uniformly in all $`V(n)𝕍(n)`$.

Bass and Griffin  also asked about the exact rate of escape of the transient process of the favourite sites. This seems to be a very challenging problem. Here we formulate it in a weaker form.

###### Question 2.1

Find the value of the constant $`\gamma _0`$ such that with probability one,

$$\underset{n\mathrm{}}{lim\; inf}\frac{|V(n)|}{n^{1/2}(\mathrm{log}n)^\gamma }=\{\begin{array}{cc}0\hfill & \text{if }\gamma <\gamma _0\text{,}\hfill \\ \mathrm{}\hfill & \text{if }\gamma >\gamma _0\text{.}\hfill \end{array}$$

According to (2.1), we must have $`\gamma _0[1,11]`$. There is good reason to expect that $`\gamma _0`$ would lie in $`[1,2]`$.

### 2.2 Upper rates of favourite sites

If we are interested in the upper limits of $`|V(n)|`$, here is an answer in the form of the law of the iterated logarithm (LIL), which was discovered independently by Erdős and Révész  and Bass and Griffin .

Theorem B (, ). With probability one,

$$\underset{n\mathrm{}}{lim\; sup}\frac{V(n)}{(2n\mathrm{log}\mathrm{log}n)^{1/2}}=1,\text{a.s.}$$

Therefore, the favourite site $`V(n)`$ and the random walk $`S(n)`$ satisfy the same LIL. However, a closer inspection by Erdős and Révész  reveals that they have different Lévy upper functions:

Theorem C (). There exists a deterministic sequence $`(a_n)_{n1}`$ satisfying

$$\{S(n)a_n,\text{i.o.}\}=1,\{V(n)a_n,\text{i.o.}\}=0.$$

The upper functions of $`S(n)`$ are characterized by the Erdős–Feller–Kolmogorov–Petrowsky integral test (Révész , p. 35): if $`(a_n)_{n1}`$ is a non-decreasing sequence of positive numbers, then

(2.2)
$$\{S(n)n^{1/2}a_n,\text{i.o.}\}=\{\begin{array}{cc}0\hfill & \\ 1\hfill & \end{array}\underset{n}{}\frac{a_n}{n}\mathrm{exp}\left(\frac{a_n^2}{2}\right)\{\begin{array}{cc}<\mathrm{},\hfill & \\ =\mathrm{}.\hfill & \end{array}$$

Erdős and Révész  and Bass and Griffin  asked the following question:

###### Question 2.2

Find an integral test to decide whether $`\{V(n)a_n,\text{i.o.}\}=0`$.

We believe that a key to Question 2.2 would be to control the upper tail probability of the favourite site, formulated here for Brownian motion: let $`W`$ be a standard Brownian motion whose local time process is denoted by $`(L(t,x);t_+,x)`$, i.e., for any bounded Borel function $`f`$,

$$_0^tf(W(s))ds=_{}f(x)L(t,x)dx.$$

Let $`U`$ denote the (almost surely) unique favourite site at time 1: $`L(1,U)=sup_xL(1,x)`$. We pose the following

###### Conjecture 2.3

There exists a constant $`\nu >1`$ such that

(2.3)
$$0<\underset{x+\mathrm{}}{lim\; inf}x^\nu \mathrm{e}^{x^2/2}(U>x)\underset{x+\mathrm{}}{lim\; sup}x^\nu \mathrm{e}^{x^2/2}(U>x)<+\mathrm{}.$$

If (2.3) holds, then we think that we should be able to obtain an integral test characterizing the upper functions of $`V(n)`$:

###### Conjecture 2.4

Let $`\nu >1`$ be the constant satisfying $`(\text{2.3})`$. For any non-decreasing sequence $`(a_n)_{n1}`$ of positive numbers,

(2.4)
$$\{V(n)n^{1/2}a_n,\text{i.o.}\}=\{\begin{array}{cc}0\hfill & \\ 1\hfill & \end{array}\underset{n}{}\frac{a_n^{2\nu }}{n}\mathrm{exp}\left(\frac{a_n^2}{2}\right)\{\begin{array}{cc}<\mathrm{},\hfill & \\ =\mathrm{}.\hfill & \end{array}$$

Let us explain why we conjecture $`\nu >1`$. By the trivial inequality $`Usup_{0t1}W(t)`$ and the usual estimate for Gaussian tails, it is easily seen that if the first inequality in (2.3) holds, then $`\nu 1`$. On the other hand, if $`\nu `$ were equal to 1, then (2.4) would be the same as the Erdős–Feller–Kolmogorov–Petrowsky test in (2.2) which would contradict Theorem C. This leads us to the conjecture $`\nu >1`$.

More discussions upon the distribution of $`U`$ can be found in Subsection 3.5 below.

## 3 Ten Erdős–Révész questions

Erdős and Révész in  and  raised many questions about favourite sites of random walk, which we quote below. They correspond exactly to Questions 1–10 on pp. 130–131 of the book of Révész . Recall from (1.1) that $`\xi (n,x)`$ is the local time of the random walk, and from (1.2) that $`𝕍(n)`$ is the set of favourite sites. As before, $`V(n)`$ denotes an arbitrary element of $`𝕍(n)`$.

### 3.1 Joint behaviour of favourite site and maximum local time

Let $`\xi ^{}(n):=sup_x\xi (x,n)`$, which is the maximum local time. In their proof of Theorem B (cf. Subsection 2.2), Erdős and Révész  noticed that, for any $`\epsilon >0`$, almost surely there exist infinitely many $`n`$ such that simultaneously, $`V(n)(1\epsilon )(2n\mathrm{log}\mathrm{log}n)^{1/2}`$ and $`\xi ^{}(n)c(2n\mathrm{log}\mathrm{log}n)^{1/2}`$, for some positive constant $`c`$ depending on $`\epsilon `$. This led them to ask the following question: what can be said about the joint asymptotics of $`V(n)`$ and $`\xi ^{}(n)`$? In particular, if $`V(n)`$ is close to its maximal possible value, how large can $`\xi ^{}(n)`$ be?

If $`V(n)`$ and $`\xi ^{}(n)`$ were asymptotically independent, then one would expect that the limit set of $`\{V(n)/(2n\mathrm{log}\mathrm{log}n)^{1/2},\xi ^{}(n)/(2n\mathrm{log}\mathrm{log}n)^{1/2}\}`$ should be the half-disc $`\{(x,y):y0,x^2+y^21\}`$. However, the correct answer provided by Csáki et al.  shows that things do not go exactly like this.

Theorem D (). With probability one, the random sequence

$$(\frac{V(n)}{(2n\mathrm{log}\mathrm{log}n)^{1/2}},\frac{\xi ^{}(n)}{(2n\mathrm{log}\mathrm{log}n)^{1/2}})_{n3}$$

is relatively compact, whose limit set is identical to the simplex $`\{(x,y):y0,|x|+y1\}`$.

In particular, Theorem D implies Theorem B, and also the LIL for the maximum local time of random walk which was originally proved by Kesten .

The proof of Theorem D relies on an invariance principle for local times due to Révész  (recalled in (3.3) below) and on the Ray–Knight theorem for Brownian local time.

### 3.2 Many favourite sites

At each step (say, $`n(n+1)`$) of the random walk exactly one of the following three possibilities occurs:

* the currently occupied site is not favourite, $`S(n+1)𝕍(n+1)`$, and thus $`𝕍(n+1)=𝕍(n)`$ remains unchanged;
* the currently occupied site becomes a new favourite besides the favourites of the previous time $`n`$, thus $`𝕍(n)𝕍(n+1)`$ and $`𝕍(n+1)𝕍(n)=\{S(n+1)\}`$;
* the random walk revisits a site which was already favourite in the previous time $`n`$, and so this new site becomes the only new favourite $`𝕍(n+1)=\{S(n+1)\}𝕍(n)`$.

It follows that the number of favourite sites either remains unchanged, or increases by one, or drops down to 1. From the recurrence of the random walk it follows easily that $`\{\mathrm{\#}𝕍(n)=1,\text{ i.o.}\}=1`$ and $`\{\mathrm{\#}𝕍(n)=2,\text{ i.o.}\}=1`$. The question is

(3.1)
$$\{\mathrm{\#}𝕍(n)=r,\text{ i.o.}\}=\mathrm{?}r=3,4,\mathrm{}$$

Erdős and Révész  conjectured that this probability should be 0, for any $`r3`$.

Recently, Tóth  proved the following

Theorem E (). We have,

$$\{\mathrm{\#}𝕍(n)4,\text{i.o.}\}=0.$$

The main argument of the proof is based on some (quite natural) rearrangements of sums and the Ray–Knight representation of the local time process stopped at inverse local times. On the technical level, the proof of Theorem E relies on controlling the probability distribution of the number of global maxima of given height $`h1`$ of critical Galton–Watson processes.

### 3.3 Frequency of having many favourite sites

It is intuitively clear that in “most situations”, there is only one favourite site. Erdős and Révész  were interested in the question about how often there are at least two favourite sites. To formulate their question precisely, let

$$\nu _0:=0,\nu _{k+1}:=inf\{n>\nu _k:\mathrm{\#}𝕍(n)>1\},$$

that is: $`\nu _k`$ is the $`k`$-th time when there are more than one favourite sites. Or, alternatively, we can define

$$\kappa _n:=\mathrm{\#}\{0<jn:\mathrm{\#}𝕍(j)>1\}.$$

Now, the question is:

$$\text{ Is it true that }\underset{k\mathrm{}}{lim}\frac{\nu _k}{k}=\mathrm{},\text{ or, equivalently, }\underset{n\mathrm{}}{lim}\frac{\kappa _n}{n}=0,\text{ with probability one?}$$

The question remains open.

### 3.4 Big jumps of favourite sites

In order to formulate precisely questions 4, 8 and 9 of the Erdős–Révész list we define the last visited favourite site $`\mathrm{}(n)𝕍(n)`$ inductively, as follows:

(3.2)
$$\mathrm{}(0):=0,\mathrm{}(n+1):=\{\begin{array}{cc}\mathrm{}(n)\hfill & \text{if }S(n+1)𝕍(n+1)\text{,}\hfill \\ S(n+1)\hfill & \text{if }S(n+1)𝕍(n+1)\text{.}\hfill \end{array}$$

The fourth question concerns how large the jump sizes of favourite site can be. The answer is formulated in the following LIL.

Theorem F (). With probability one,

$$\underset{n\mathrm{}}{lim\; sup}\frac{|\mathrm{}(n+1)\mathrm{}(n)|}{(2n\mathrm{log}\mathrm{log}n)^{1/2}}=1.$$

In particular, Theorem F tells us that the extraordinarily large jumps of favourite site are asymptotically comparable to the size of the range of the random walk.

The main ingredient in the proof of Theorem F is the Ray–Knight theorem.

### 3.5 Limit law of favourite sites

The question here is: what is the limit distribution of $`V(n)/n^{1/2}`$ when $`n\mathrm{}`$?

Here is an argument to show that $`V(n)/n^{1/2}`$ has a non-degenerate limit distribution. Indeed, according to a theorem of Révész , possibly in an enlarged probability space, there exists a coupling for random walk $`(S(n),n_+)`$ and Brownian motion $`(W(t),t_+)`$ such that for any $`\epsilon >0`$,

(3.3)
$$\underset{x}{sup}|\xi (n,x)L(n,x)|=o\left(n^{1/4+\epsilon }\right),\text{a.s.},$$

where $`\xi `$ and $`L`$ denote the local times of $`S(n)`$ and $`W(t)`$, respectively. Thus, for any fixed $`a>0`$,

$`\left({\displaystyle \frac{V(n)}{n^{1/2}}}>a\right)`$ $``$ $`\left(\underset{x>n^{1/2}a}{sup}\xi (n,x)\underset{xn^{1/2}a}{sup}\xi (n,x)\right)`$
$``$ $`\left(\underset{x>n^{1/2}a}{sup}L(n,x)\underset{xn^{1/2}a}{sup}L(n,x)n^{1/4+\epsilon }\right)`$
$`=`$ $`\left(\underset{y>a}{sup}L(1,y)\underset{ya}{sup}L(1,x)n^{1/4+\epsilon }\right)`$
$``$ $`\left(\underset{y>a}{sup}L(1,y)\underset{ya}{sup}L(1,y)\right)=\left(U>a\right),`$
when $`n\mathrm{}`$, where $`U`$ denotes as in Subsection 2.2 the location of the maximum (on $``$) of $`xL(1,x)`$: $`L(1,U)=sup_xL(1,x)`$. Similarly,

$`\left({\displaystyle \frac{V(n)}{n^{1/2}}}>a\right)`$ $``$ $`\left(\underset{x>n^{1/2}a}{sup}\xi (n,x)>\underset{xn^{1/2}a}{sup}\xi (n,x)\right)`$
$``$ $`\left(\underset{x>n^{1/2}a}{sup}L(n,x)>\underset{xn^{1/2}a}{sup}L(n,x)+n^{1/4+\epsilon }\right)`$
$``$ $`\left(U>a\right),`$
so that,

$$\underset{n\mathrm{}}{lim}\frac{V(n)}{n^{1/2}}=U,\text{in distribution.}$$

We are grateful to Endre Csáki for having communicated to us this simple argument for the weak convergence.

The distribution of $`U`$ was characterized by Theorem 6.2 of Borodin , where a double Laplace transform of the limit distribution was computed, namely, he got a close-form expression for $`𝔼(\mathrm{e}^{a\sqrt{A}|U|})`$ for $`a>0`$ and an exponential variable $`A`$ which is independent of $`U`$. However, the expression for $`𝔼(\mathrm{e}^{a\sqrt{A}|U|})`$ found in  looks very complicated, involving ratios of Whittaker functions. We have not been able to invert the double Laplace transform, or even to get reasonably good information for the tail probability of $`U`$ which would be useful for the upper functions of $`V(n)`$, see Conjectures 2.3 and 2.4 in Subsection 2.2.

### 3.6 Total number of favourite sites

Let $`\alpha (n)`$ denote the number of all the different favourite sites up to step $`n`$, i.e., $`\alpha (n)=\mathrm{\#}(_{k=0}^n𝕍(k))`$. Is it true that with probability one, for all large $`n`$, $`\alpha (n)(\mathrm{log}n)^c`$ (for some constant $`c>0`$)?

The question is still open. Omer Adelman (personal communication) has a lower bound for $`\alpha (n)`$, proving that with probability one, for all large $`n`$, $`\alpha (n)c^{}\mathrm{log}n`$, for some constant $`c^{}>0`$.

### 3.7 Durations of favourite sites

The question is: how long a favourite site can stay favourite? More precisely, let $`\beta (n):=\mathrm{max}\{ji:\mathrm{\hspace{0.17em}0}ijn,_{k=i}^j𝕍(k)\mathrm{}\}`$. In words, $`\beta (n)`$ is the duration of the longest period (before $`n`$) during which a favourite site stays favourite. What can be said about the asymptotic behaviour of $`\beta (n)`$?

No answer available so far.

### 3.8 “Capricious” favourite sites

If $`x`$ is a favourite site at some stage, can it happen that the favourite site moves away from $`x`$ but later returns to $`x`$? More precisely, do infinite random sequences $`\mathrm{}<c_{n1}<a_n<b_n<c_n<a_{n+1}<\mathrm{}`$ of positive integers exist such that $`\mathrm{}(a_n)=\mathrm{}(c_n)\mathrm{}(b_n)`$, for $`n=1,2,\mathrm{}`$? (Recall the definition of the last visited favourite site, (3.2).)

The question remains open.

### 3.9 Small jumps of favourite sites

Consider the random increasing sequence of times, when the last visited favourite site changes value

$$\lambda _0:=0,\lambda _{k+1}:=inf\{n>\lambda _k:\mathrm{}(n)\mathrm{}(n1)\},$$

and the jump sizes of $`\mathrm{}(n)`$ at these times:

$$j_k:=|\mathrm{}(\lambda _{k+1})\mathrm{}(\lambda _k)|.$$

It seems likely that $`j_k\mathrm{}`$ almost surely, with $`k\mathrm{}`$. How to describe the limit behaviour of $`j_k`$? This question is a companion to the one in Subsection 3.4. While in Subsection 3.4 the upper behaviour of the jump size was determined with satisfying precision, we have not been able to get any non-trivial information about the lower behaviour.

### 3.10 Occupation times and favourite sites

The arcsine law says that with big probability, the random walk spends a long time on one half of the line (say, $`_+`$) and only a short time on the other half (in this case $`_{}`$). Is it true that the favourite site is located on the same side where the random walk spends the most time? For example, if $`(\mu _k)_{k1}`$ is a random sequence of integers satisfying $`\mu _k^1_{j=1}^{\mu _k}\mathrm{𝟏}\mathrm{l}_{\{S_j0\}}1`$ ($`k\mathrm{}`$), then Erdős and Révész  conjectured that $`V(\mu _k)+\mathrm{}`$, a.s.

The answer to the conjecture is no. Relying on the Ray–Knight theorem and careful analysis of the sample paths of Bessel processes, Csáki and Shi  prove the existence of a random sequence $`(\mu _k,k1)`$ such that $`\mu _k^1_{j=1}^{\mu _k}\mathrm{𝟏}\mathrm{l}_{\{S_j<0\}}c/(\mathrm{log}\mathrm{log}\mu _k)^2`$ (for some constant $`c>0`$ and all $`k`$), and yet $`V(\mu _k)<0`$ for all $`k`$.

## 4 Some related questions for simple random walk

### 4.1 Rarely visited sites

(1) Rarely visited sites. Let $`R(n):=\{S(0),S(1),\mathrm{},S(n)\}`$ be the range of the random walk up to step $`n`$. Erdős and Révész  conjectured that for any integer $`r3`$, the probability that there are infinitely many $`n`$ such that each of the sites of $`R(n)`$ has been visited at least $`r`$ times up to step $`n`$ is 0. This was disproved by Tóth , who showed that for any integer $`r`$, this probability actually equals 1. The proof of this result again relies on the Ray–Knight representation of local times. This time one has to control the probability of the event that the value of the critical Galton–Watson process drops down to zero from a given value $`r`$, uniformly in the initial condition.

Another aspect of rarely visited sites was studied by Major . He proved that, if $`Z(n)`$ denotes the number of sites of $`R(n)`$ which have been visited exactly once, then with probability one, $`lim\; sup_n\mathrm{}Z(n)/(\mathrm{log}n)^2=c`$, where $`c`$ is a finite and positive constant.

### 4.2 Favourite edges

(2) Favourite edges. Instead of looking at favourite sites (which by definition maximize the site local time of the random walk), we look at favourite edges which maximize the edge local time. We can ask a similar question as Question 3.2: are there infinitely many $`n`$ such that there are at least 3 favourite edges at time $`n`$? Tóth and Werner  proved that with probability one, there are at most finitely many $`n`$ such that there are at least 4 favourite edges at time $`n`$. This is a simpler predecessor of Theorem E, cited above.

### 4.3 Location of favourite sites

(3) Location of favourite sites. The question is whether it is possible for a favourite site to be close to the boundary of the range $`R(n)`$, i.e., close to $`\mathrm{max}_{0kn}S(k)`$ or to $`\mathrm{min}_{0kn}S(k)`$. (The question was originally due to Omer Adelman). Csáki and Shi  proved that the answer is no: with probability one, $`\mathrm{max}_{0kn}S(k)V(n)`$ (and by symmetry, $`V(n)\mathrm{min}_{0kn}S(k)`$) goes to infinity, and the rate of escape was determined.

## 5 Favourite sites of other processes

(1) Brownian motion. Many questions and results mentioned above can be formulated for the favourite sites of Brownian motion in the obvious way. Some further discussions and questions can be found in Leuridan .

(2) Symmetric stable processes. In Bass et al. , it was proved that the favourite site of a symmetric stable process is also transient. Eisenbaum  proved a collection of interesting results for favourite sites of symmetric stable processes. For example, she showed that for each given $`t0`$ there is almost surely a unique favourite site at time $`t`$, and that with probability one, for all $`t0`$, there are at most two favourite sites at time $`t`$.

(3) Lévy and Markov processes. Most of the results for symmetric stable processes can be extended to a larger class of Lévy and even symmetric Markov processes. See Marcus , Eisenbaum and Khoshnevisan .

(4) Two-dimensional random walk. Let $`V_2(n)`$ denote a favourite site, at time $`n`$, of a simple symmetric random walk on $`^2`$. Dembo et al.  proved that $`(\mathrm{log}V_2(n))/\mathrm{log}n`$ converges to $`1/2`$ with probability one, where $`x`$ denotes the Euclidean modulus in $`^2`$.

(5) Transient processes. If $`\{X(t),t0\}`$ is a transient process, having local time at $`t=\mathrm{}`$, denoted by $`L(\mathrm{},x)`$, then $`\{x[0,T]:L(\mathrm{},x)=sup_{y[0,T]}L(\mathrm{},y)\}`$ is the set of favourite sites of $`X`$ in $`[0,T]`$. Bertoin and Marsalle  studied the case when $`X`$ is a Brownian motion with a positive drift, and Hu and Shi  the case when $`X`$ is the modulus of a $`d`$-dimensional Brownian motion ($`d>2`$). They obtained respective rates of escape (when $`T\mathrm{}`$) of favourite sites.

(6) Poisson process. Khoshnevisan and Lewis  studied favourite sites of a Poisson process. They obtained several laws of the iterated logarithm.

(7) Random walk in random environment. The favourite sites can be defined for nearest neighbour random walk in random environment, exactly as for the usual random walk. Hu and Shi  considered the recurrent case and proved that the process of favourite sites is again transient, and the escape rate was characterized via an integral test. The latter question is still open for the usual simple random walk, see Question 2.1 above. The problem of escape rates of favourite sites is the only problem we are aware of, which is solved for random walk in random environment, but which is open for the usual random walk.

|  | Zhan Shi | Bálint Tóth |
| --- | --- | --- |
|  | Laboratoire de Probabilités UMR 7599 | Institute of Mathematics |
|  | Université Paris VI | Technical University Budapest |
|  | 4 place Jussieu | Egry J. u. 1 |
|  | F-75252 Paris Cedex 05 | H-1111 Budapest |
|  | France | Hungary |
|  | zhan@proba.jussieu.fr | balint@math.bme.hu |


format: text

text:
# Forecasting price increments using an artificial Neural Network

## 1 Introduction

Forecasting future values of an asset gives, besides the straightforward profit opportunities, indications to compute various interesting quantities such as the price of derivatives (complex financial products) or the probability for an adverse mode which is the essential information when assessing and managing the risk associated with a portfolio investment.

Forecasting the price of a certain asset (stock, index, foreign currency, etc.) on the ground of available historical data, corresponds to the well known problem in science and engineering of time series prediction. While many time series may be approximated with a high degree of confidence, financial time series are found among the most difficult to be analyzed and predicted. This is not surprising since the dynamics of the markets following at least the semi-strong EMH should destroy any easy method to estimate future activities using past informations.

Among the methods developed in *Econometrics* as well as other disciplines ∗\*∗\*$``$ see the vast bibliography with more than 800 entries at
www.stern.nyu.edu/ aweigend/Time-Series/Biblio/SFIbib.html reported from \[Weigend and Gershenfeld (1994)\], the artificial *Neural Networks* (NN) are being used by “non-orthodox” scientists as non-parametric regression methods \[Campbell, Lo and MacKinlay (1997), Moody and Neuneier + Zimmermann (1998)\]. They constitute an alternative to nonparametric regression methods like *kernel regression* \[Campbell, Lo and MacKinlay (1997)\]. The advantage of using a neural network as non linear function approximator is that it appears to be well suited in areas where the mathematical knowledge of the stochastic process underlying the analyzed time series is unknown and quite difficult to be rationalized. Besides, it is important to note that the lack of linear correlations in the financial price series and the already accepted evidence of an underlying process different from i.i.d. noise point out to the existence of higher-order correlations or non-linearities. It is this non-linear correlation that the neural net may eventually catch during its learning phase. If some macroscopic regularities, arising from the apparently chaotic behaviour of the large amount of components are present, then a well trained net could identify and “store” them in its distributed knowledge representation system made by units and synaptic weights \[Moody and Neuneier + Zimmermann (1998), Refenes, Burgess and Bentz (1997)\].

In the following we will see that a well suited NN for each of a set of price time series showing a “surprising” rate of success in predicting the *sign* of the price change on a daily base *can* be found. Not less interesting, we will see that the foretold regularities in the time series seem to be more present on larger time scale than on high frequency data, as the performance of the net degrades if we go from monthly to minutes data.

## 2 Multi-layer Perceptron

*Multi-layer perceptrons* (MLP) are the neural nets usually referred to as function approximators. A MLP is a generalization of Rosenblatt’s *perceptron* (1958); $`n_i`$ input units, $`n_h`$ hidden and $`n_o`$ output units with all feed forward connections between adjacent layers (no intra-layer connections or loops). Such net’s topology is specified as $`n_i`$-$`n_h`$-$`n_o`$.

A NN may perform various tasks connected to classification problems. Here we are mainly interested in exploiting what is called the *universal approximation property*, that is, the ability to approximate any nonlinear function to any arbitrary degree of accuracy with a suitable number of hidden units \[White (1992), Cybenko (1989)\].

The approximation is performed finding the set of weights connecting the units. This can be done with one of the available methods of non-parametric estimation techniques like *nonlinear least-squares*. In particular we choose *error back propagation* which is probably the most used algorithm to train MLPs \[(7), Rumelhart, Hinton and Williams (1986)\]. It is basically a gradient descent algorithm of the error computed on a suitable learning set. A variation of it use *bias*, *terms* and *momentum* as characteristic parameters. Moreover we fixed the learning rate $`\eta =0.05`$, the momentum $`\beta =0.5`$ and the usual sigmoidal as nonlinear activation function.

## 3 Detrending analysis

We have trained the neural nets on “detrended” time series. The detrending analysis was performed to mitigate the unbalance between the *learning set*, and the *test set*. In fact, subdividing the available data in learning set and testing set as specified in the following section (have a look at figure 1), we train the nets on a data set corresponding to a periods much back in time while we test the nets on data set corresponding to the most recent period of time. This problem is know in literature as *noise/nonstationarity tradeoff* \[Moody and Neuneier + Zimmermann (1998)\].

It is known that during the last ten years the American market has noticeably changed in that almost all the titles connected to the information technology have not only jumped to record values but also the fluctuations of price today are much stronger than before ††††††$`\mathrm{}`$ $`P_t`$ is what we use to train our nets. Considering $`\mathrm{log}(P_t)`$ instead of $`P_t`$ would mitigate the problem but it would introduce further nonlinearities . Ignoring this fact would lead to a mistake because the net would not learn the characteristics of the “actual situation”.

To detrend a time series we performed a nonlinear least squares fit using the *Marquardt-Levenberg* algorithm \[Campbell, Lo and MacKinlay (1997), Press, Teukolsky, Vetterling and Flannery (1994)\] with a polynomial of sixth degree. Then we just computed the difference of the series with the fitting curve. For each time series considered we ended up with a detrended series composed by 2024 points corresponding to the period from about January 1990 to February 2000. For example, the plot in figure 2 shows the detrended time series of the index S&P500 along with the original series and the polynomial fit.

We choose daily closing for 3 indexes and 14 assets historical series on the NYSE and Nasdaq. In particular the assets were chosen among the most active companies in the field of information technology.

## 4 Determining the net topology

One of the primary goals in training neural networks is to ensure that the network will perform well on data that it has not been trained on (called ”generalization”). The standard method of ensuring good generalization is to divide our training data into *multiple data sets*. The most common data sets are the *learning* $`L`$, *cross validation* $`V`$, and *testing* $`T`$ data sets. While the learning data set is the data that is actually used to train the network the usage of the other two may need some explanation.

Like the learning data set, the cross validation data set is also used by the network during training. Periodically, while training on the learning data set, the network is tested for performance on the cross validation set. During this testing, the weights are not trained, but the performance of the network on the cross validation set is saved and compared to past values. If the network is starting to overtrain on the training data, the cross validation performance will begin to degrade. Thus, the cross validation data set is used to determine when the network has been trained as well as possible without overtraining (e.g., maximum generalization).

Although the network is not trained with the cross validation set, it uses the cross validation set to choose a ”best” set of weights. Therefore, it is not truly an *out-of-sample* test of the network. For a true test of the performance of the network the testing data set $`T`$ is used. This data set is used to provide a true indication of how the network will perform on new data.

In figure 3, an example of MLP with $`n_i=3`$, $`n_h=7`$ and one output unit takes $`P_{t_0},P_{t_1},P_{t_2}`$ in input and gives the successive value $`P_{t_3}`$ as forecast. The number of free parameters is given by the number of connections between units $`(n_i+n_o)n_h`$.

While the choice of one output unit comes from the straightforward definition of the problem, a crucial question is “how many input and hidden units should we choose?”. In general there is no way to determine apriori a good network topology. It depends critically on the number of training examples and the complexity of the time series we are trying to learn. To face this problem a large number of methods are being developed (recurrent networks, model selection and pruning, sensitivity analysis \[Moody and Neuneier + Zimmermann (1998)\]), some of which follow the evolution’s paradigm (Evolutionary Strategies and Genetic Algorithm).

Because we have observed a critical dependence of the performance of the net from $`n_i`$ and $`n_h`$, and to avoid the great complexity of more powerful strategies \[Moody and Neuneier + Zimmermann (1998)\], we ended up with the decision to explore all the possible combinations of $`n_i`$-$`n_h`$ in a certain range of values. Our “brute force” procedure consists of training nets of different topologies (varying $`2n_i15`$ and $`2n_h25`$) and observe their performance. More precisely we select good nets on the basis of the mean square error (see eq(4.1)) computed on 200 points out of the sample set constituting the test set. Thus, besides the separation in Learning-Validation-Testing of our time series, we further distinguish a subset from the Testing set: the *Checking* $`C`$ (see fig. 1). The reason is that while we train the net to interpolate the time series (minimizing the mean square error) we finally extrapolate to forecast the *sign of the increments* (to be defined later).

To assess the efficiency of the learning and to discard bad trained nets during the search procedure we use the *mean square error* $`ϵ`$ defined as

$$ϵ=\frac{1}{\sigma }\frac{1}{|C|}\underset{tC}{}\left(G_tP_t\right)^2$$
(4.1)
where $`P_t`$ is the price value, $`G_t`$ is the forecasted value at time $`tC`$ and $`\sigma `$ is the standard deviation of the time series. For good forecasts we will have small positive values of $`ϵ`$ ($`1ϵ0`$).

We set the threshold 0.015 to discriminate good from bad nets. Only those nets for which $`ϵ0.015`$ are further tested for sign prediction.

In summary, first we learn on set $`L`$, and through validation $`V`$ we find when to stop learning; then through check on $`C`$ we see if the learning process worked well, and in case it did, we make predictions in the test phase on set $`T`$ for ”future” (i.e. previously unused) price changes and compare them with reality.

## 5 Stopping criteria

To avoid overfitting and/or very slow convergence of the training phase, the stopping criteria is determined by the following three conditions, one of which is sufficient to end the training phase (early stopping):

1. Stopping is assured within 5000 iterations of cross validation (see section 4);
2. during cross validation the mean square error on the validation set $`V`$ is computed as $`\epsilon _V=\frac{1}{2}_{tV}\left(G_tP_t\right)^2`$; during training $`\epsilon _V`$ should decrease, so a stopping condition is given if $`\epsilon _V`$ increase again more than 20% of the minimum value reached up to then;
3. learning is also stopped if $`\epsilon _V`$ reaches a plateau; this is tested during cross validation averaging 1000 successive values of $`\epsilon _V`$ and checking if the actual value is above this average.

## 6 Results

The plot in figure 4 compares the forecasted $`G_t`$ and the real $`P_t`$ values for the time series of Apple Corp. on the test set $`T`$. It also shows a linear fit for the points $`\{P_t,G_t\}`$. A raw measure of performance on the test set $`T`$ can be obtained by the slope of the fitting line (let’s call it $`\theta `$). It will be a value close to one if the fit corresponds to the $`y=x`$ line, i.e., if $`P_t=G_t`$. We obtained the following $`\theta `$’s for the time series in table 2 and 3: $`\theta _{\text{S\&P500}}=0.906`$, $`\theta _{\text{DJI}}=0.874`$, $`\theta _{\text{Nasdaq100}}=0.860`$. $`\theta _{\text{AAPL}}=0.976`$, $`\theta _\text{T}=0.921`$, $`\theta _{\text{AMD}}=0.914`$, $`\theta _{\text{STM}}=0.885`$, $`\theta _{\text{HON}}=0.885`$, $`\theta _{\text{INTC}}=0.874`$, $`\theta _{\text{CSCO}}=0.860`$, $`\theta _{\text{WCOM}}=0.847`$, $`\theta _{\text{IBM}}=0.842`$, $`\theta _{\text{ORCL}}=0.824`$, $`\theta _{\text{MSFT}}=0.803`$, $`\theta _{\text{SUNW}}=0.774`$, $`\theta _{\text{DELL}}=0.692`$, $`\theta _{\text{QCOM}}=0.488`$.

The final estimation of the performance in forecasting is made by means of the *one-step sign prediction rate* $`\zeta `$ defined on $`T`$ as follows

$$\zeta =\frac{1}{|T|}\underset{tT}{}HS(\mathrm{\Delta }P_t\mathrm{\Delta }G_t)+1HS(|\mathrm{\Delta }P_t|+|\mathrm{\Delta }G_t|)$$
(6.1)
where $`\mathrm{\Delta }P_t=P_tP_{t1}`$ the price change at time step $`tT`$ and $`\mathrm{\Delta }G_t=G_tP_{t1}`$ is the *guessed* price change at the same time step. Note that we assume to know the value of $`P_{t1}`$ to evaluate $`\mathrm{\Delta }G_t`$. $`HS`$ is a modified ‡‡‡‡‡‡$`\mathrm{}`$ The usual $`HS`$ function gives 1 in zero, i.e., $`HS(0)=1`$  Heaviside function $`HS(x)=1`$ for $`x>0`$ and 0 otherwise. The argument of the summation in eq(6.1) gives one only if $`\mathrm{\Delta }P_t`$ and $`\mathrm{\Delta }G_t`$ are non-zero and with same sign, or if $`\mathrm{\Delta }P_t`$ and $`\mathrm{\Delta }G_t`$ are both zero. In other words $`\zeta `$ is the probability of a correct guess on the sign of the price increment estimated on $`T`$.

In the lower-right inset of figure 4 it is shown $`P_tG_t`$ as function of $`P_t`$. One can see that the difference between the real and the forecasted values clusters for small $`P_t`$. Another way see it is to look at the histogram of $`\zeta `$ as function of $`\mathrm{\Delta }P_t`$. In other words the rate of correct guesses on the sign of the price increment relative to the magnitude of the fluctuation of the real price.

To obtain an unbiased histogram we have to normalize it dividing each bin by the corresponding value of the $`\mathrm{\Delta }P`$’s histogram (the limit of $`\mathrm{\Delta }P`$ follows a power law so that large fluctuations are much less probable). The resulting distribution is plotted in figure 5. It is now clearly visible that the net does not favor large increments over small ones or vice versa. In fact the probability to make a correct guess on the sign of the increment seems independent from the magnitude of the increment itself. This does not means that the net forecasts “rare events” (i.e., a profit opportunity) as easily as normal fluctuation, because the statistics here calculated are not significant with respect to extreme events.

To interpret the results that we are going to show we have to concentrate our attention on the way we select a good net to be used to make forecast. For each time series we have performed a search to determine the topology of a good net as specified in the last section. Once we get a pool of candidates the question is “how many of them give a sign prediction rate above fifty percent?”

This question is answered in table 1. There, *tot* indicates the number of nets such that $`ϵ0.015`$, that is, we judged as good nets, while *ok* is the number of them that gave $`\zeta 50`$. This ratio can be seen as an estimation of the confidence that the net will perform a “sufficient” forecast of price change, where sufficient means above fifty percent.

The value of $`\zeta `$ together with the specification of the number of units per layer of the best net is reported in table 2 and table 3 along with the dimension of the learning and validation set.

The sign prediction rates range from 50.29% to 54%. While the smallest values 50.29 may be questionable, the larger values above 54 seem a clear indication that the net is not behaving randomly. Instead it has captured some regularities in the nonlinearities of the series.

A quite direct test for randomness can be done computing the probability that such forecast rate can be obtained just by flipping a coin to decide the next price increment. For this purpose we use a random walk ($`pr(\text{up})=pr(\text{down})=1/2`$) as forecasting strategy $`G_{rw_t}`$ and observe how many, over 1000 different random walks, give a sign prediction rate $`\zeta _{rw}`$ defined in eq(6.1) above the value obtained with our net. Note that each random walk perform about 1000 time steps, the same as $`|T|`$ for that specified time series (see table 2 and 3). These values are reported in table 4. They indicate that except for QCOM the random walk assumption cannot give the same prediction rate as the neural net §§§§§§$`\mathrm{\S }`$ In other words, given a neural net which produce $`\zeta `$ as prediction rate over a certain time series $`P_t`$ we may compute the probability at which the *null hypothesis of randomness* is rejected. We use a random walk ($`pr(\text{up})=pr(\text{down})=1/2`$) as forecasting strategy $`G_{rw_t}`$ and then compute $`\zeta _{rw}`$ defined in eq(6.1) on the time series $`P_t`$. The random variable $`\zeta _{rw}`$ have mean 0.5 and standard deviation $`\sigma _{\zeta _{rw}}`$. By definition $`\zeta _{rw}`$ is the sample mean of $`T`$ i.i.d. Bernoullian random variables. Thus, assuming that $`\zeta _{rw}`$ converges to a Gaussian $`N(1/2,\sigma _{\zeta _{rw}})`$, we can estimate the unknown variance of $`\zeta _{rw}`$ as $`\widehat{\sigma }_{rw}^2=1/N_{i=1}^N(\zeta _{rw_i}1/2)^2`$. To have an estimation of $`\sigma _{\zeta _{rw}}`$ we ran $`N=1000`$ random walks each giving a value for $`\zeta _{rw}`$. Once we estimate $`\sigma _{rw}`$, the null hypothesis becomes “what is the probability $`P_{\zeta _{rw}}[x>\zeta ]`$ that the neural net is doing a random prediction on $`P_t`$ with rate $`\zeta `$ ?” or the other way around ”what is the probability $`P_{\zeta _{rw}}[x\zeta ]`$ that the net is not doing randomly?”. Formally $`P_{\zeta _{rw}}[x\zeta ]=_{\mathrm{}}^\zeta N(\frac{1}{2},\widehat{\sigma }_{rw})(x)𝑑x`$ where $`N(\frac{1}{2},\widehat{\sigma }_{rw})`$ is a Gaussian and $`\widehat{\sigma }_{rw}`$ is the estimation of the standard deviation $`\sigma _{rw}`$ of the random variable $`\zeta _{rw}`$. In summary, for every sign prediction rate $`\zeta `$ obtained with our neural net on a time series $`P_t`$, we first estimate $`\widehat{\sigma }_{rw}`$ as specified above, then we compute the probability $`P_{\zeta _{rw}}[x\zeta ]`$ at which the null hypothesis of randomness prediction is rejected. The results tell us that for some bad prediction values (like for QCOM or Nasdaq100) the randomess hyphothesis cannot be rejected but for the majority of the series the probability to reject the null hypothesis is something between 0.01 and 0.1. .

## 7 Weekly and intra-day data

It is interesting to ask if the MLP may exploit regularities in the time series of price sampled at a lower/higher rate than daily. Apart from the “scaling behaviour” observed empirically in real price series we are interested in the performance of our procedure (search plus learn) when we change the time scale on which we sample the price of the assets or the index at a stock market.

To answer this question we performed the same search for the good net on the IBM and AMD stock price sampled on weekly basis as well as taking intra-day data with the frequency of one minute. Both series consisted of 2024 points, the same as the daily price series.

The outcome is that intra-day data are much difficult to be forecasted with our MLPs. In fact for both the one-minute-delay data series the search did not succeeded to find a good net; all the good nets (few) have given a sign prediction rate $`\zeta <40\%`$.

On the other hand the forecast of weekly data gave a success rate comparable with that of daily series (e.g., a $`4`$-$`2`$-1 net performed $`\zeta =51.422764`$ with $`ϵ=0.004947`$).

## 8 Artificially generated price series

As last question, and to further test the correctness of our prediction, we tried to forecast the sign of price changes of an artificially generated time series. This was generated by the the Cont-Bouchaud herding model that seems one of the simplest one able to show fat tails in histogram of returns \[Cont and Bouchaud (1999)\]. This model shows the relation between the excess kurtosis observed in the distribution of returns and the tendency of market participants to imitate each other (what is called *herd behaviour*). The model consists of percolating clusters of agents \[Stauffer (2000)\]. At a given time step a certain number of coalitions (clusters) decide what to do: they buy with probability $`a`$, sell with probability $`a`$ or stay inactive with probability $`12a`$. The demand of a certain group of traders is proportional to its size and the total price change is proportional to the difference between supply and demand.

It is clear that such a model generates unpredictable time series, and our networks should not be able to make any predictions. Indeed, when our method was applied to this series it did not succeeded to find a good net as all the tried nets performed bad on the check set $`C`$, i.e., $`ϵ>0.015`$.

## 9 Discussion

We have shown that a suitable neural net able to forecast the sign of the price increments with a success rate slightly above 50 percent on a daily basis can be found. This can be an empirically demonstration that a good net exists but we do not have a mechanism to find it with “high probability”. In other words we cannot use this method as a profit opportunity because we do not know *a priori* which net to use. Perhaps a better algorithm to search for the good topology (model selection and pruning with sensitivity analysis \[Moody and Neuneier + Zimmermann (1998)\]) would give some help. The future work will likely undertake this direction.

As final remark we have found that intra-day data are much more difficult to be forecasted with our method than daily or weekly.

### Acknowledgements:

The author wishes to acknowledge D. Stauffer and G.H. Zimmermann for useful comments and hints.


text:
# Preprint DFPD 00/TH/29, TRI-PP-00-28 Spin observables for pion production from 𝑝𝑑 collisions

## 1 INTRODUCTION

Pion-production reactions have the potential to probe the nuclear phenomena at short distance since they typically involve processes transferring large momenta to the target nucleus. On the other hand, the pion is also the principal mediator of the nuclear force, hence, meson production (or absorption) plays a fundamental role in hadron dynamics and the study of these processes may reveal aspects of the meson-baryon coupling, and more generally of meson-exchange phenomena, which would remain otherwise hidden.

Reliable calculations on meson dynamics in few-nucleon systems are extremely difficult, and usually one has to assume that various reaction mechanisms compete in dominating the process. The description of the process is complicated because the treatment of the reaction mechanisms reveals the occurrence of many terms, and this forces one to make further assumptions in order to reduce the number of terms to a few, tractable ones. This reduction clearly introduces ambiguities making it more difficult to extract information about the nuclear wave function at short distances, or about the modifications of the hadron interactions because of the presence of other nucleons. A somewhat simplified situation is found when considering nucleon-induced production close to the pion threshold, since there the $`s`$-wave mechanisms of the elementary $`NN`$$`\pi `$$`NN`$ inelasticities dominate, while the $`p`$-wave mechanisms (including the isobar degrees of freedom) can be treated as corrections.

## 2 THEORETICAL MODEL

In a recent paper , the $`\pi `$-production reaction from $`pd`$ collisions has been calculated in the threshold region. The model includes meson-exchange rescattering diagrams in $`s`$ waves, with $`\pi N`$ contributions in the isovector and the isoscalar channels. The isovector contribution was generated by a $`\rho `$-mediated diagram, and the isoscalar interaction by means of a subtracted $`\sigma `$-exchange model which is suppressed on-shell but enhanced when the rescattered pion is off-mass shell . This off-shell enhancement has been previously advocated  to describe the production cross section at threshold for the simpler reaction $`pppp\pi ^o`$ and represents an alternative explanation for this process with respect to the one given by Lee and Riska  in terms of heavy-meson exchange (HME) effects. With the same parameterization suggested by Hamilton we found in Ref.  that the pion-production cross section at threshold can be reproduced, while without these off-shell effects the results are underestimated by one order of magnitude.

According to a previous treatment  centered around the $`\mathrm{\Delta }`$-resonance region, the effects of the $`\pi N`$ $`p`$-wave interaction have also been included, via the nonrelativistic $`\pi NN`$ vertex, as well as by means of the $`\pi N\mathrm{\Delta }`$ interaction. The intermediate $`\mathrm{\Delta }`$ propagation terminates through a $`\mathrm{\Delta }N`$ transition determined by the $`\pi `$ plus $`\rho `$ exchange diagrams, where the pseudoscalar meson provides the typical longitudinal structure to the transition potential, while the vector meson generates the transverse contribution.

## 3 RESULTS

With the same model of Ref. , we have calculated the spin observables at threshold. The formalism for the spin observables for this reaction has been given in Ref. . For the 3$`N`$ bound-state in the outgoing channel, we have used the wave functions which have been calculated in Ref. . As two-nucleon input for the three-nucleon equations we have used high rank separable representations  of the Paris and the Bonn B potentials, known as PEST and BBEST, respectively. These representations were originally constructed by the Graz group . The three-nucleon dynamics in the initial state (ISI) has been calculated by solving the Faddeev-type Alt-Grassberger-Sandhas equations . More details can be found in Refs. .

As an example, the angular distribution for the proton analyzing power $`A_{y0}`$ is shown for the pion center-of-mass momentum $`P_{c.m.}`$ $`0.25m_\pi `$, close to the $`\pi `$-production threshold (Fig. 2), and also at higher energy ($`P_{c.m.}1.36m_\pi `$), close to the $`\mathrm{\Delta }`$ resonance (Fig. 2). The two lines shown in Fig. 2 include all production mechanisms we have discussed, and differ because the solid line includes also the effects of ISI, while the dotted line refers to a plane-wave calculation. The ISI effects turn out to be quite important, while at higher energy their effects are smaller. In Fig. 2 we show the modifications obtained when adding the two $`s`$-wave diagrams (“$`\rho `$” and “$`\sigma `$”) on top of the $`p`$-wave mechanisms, at the $`\mathrm{\Delta }`$ resonance. Both these two diagrams are needed for the reproduction of $`A_{y0}`$ in the forward hemisphere, while at backward angles the situation is still controversial.

Acknowledgments

L. C. and G. P. acknowledge funds from the “Ministero dell’Universitá e della Ricerca Scientifica e Tecnologica”, within the Project “Fisica Teorica del Nucleo e dei Sistemi a Più Corpi”. W. Sch. and J. P. S. acknowledge support and hospitality from the Physics Department of the University of Padova, and from INFN, sez. di Padova. The work of W. Sch. and J. P. S. is also supported by the Natural Science and Engineering Research Council of Canada.


text:
# A note on the time evolution of generalized coherent states

## I Introduction

Coherent states were originally introduced into physics by Schrödinger in order to reconcile Heisenberg’s abstract solution of the quantum harmonic oscillator with the classical picture of a swinging pendulum. Schrödinger’s minimum-uncertainty wave-packets maintain their shape under the quantum time evolution and their “centre-of-mass” motion coincides with that of the pendulum. There are various generalizations of the original harmonic oscillator coherent states, but most of them exhibit the following classical-quantum correspondence:

* There is a family of states parameterized by some classically interpretable variables.
* For simple Hamiltonians, a state originally in this family remains in it.
* The parameters evolve according to the classical equations of motion.
* In addition to the parameter evolution, the states accumulate a phase equal to the exponential of the classical action.

The last property was pointed out by van Hove. If all four of these features are present, the classical dynamics completely determines the quantum time evolution.

One such extension of the Schrödinger states is the class of “generalized coherent states” introduced by Perelomov. Perelomov’s states are obtained by selecting a fiducial vector in some representation of a Lie group, and considering its orbit under the action of the group. For compact groups (and for some representations of non-compact groups) this construction guarantees the existence of a resolution of unity, and so suggests a path integral representation of the dynamics.

Usually the Perelomov coherent states are constructed by taking a greatest (or least) weight state as the fiducial vector. This generates minimum uncertainty states and, for Hamiltonians that are elements of the Lie algebra, the states evolve classically in the manner described above. In recent articles, however, M. Matsumoto has considered coherent states built on arbitrary elements of the representation space. His motivation is that such states are useful in quantum optics and that novel phase interference effects may be measurable.

The purpose of this note is to point out what must be well known—although I have not been able to find an explicit statement in the literature—that these “general” generalized coherent states do not necessarily inherit the classical-quantum correspondence property. If we insist on the traditional picture, then we must restrict our choice of fiducial vectors. The essential requirement is that two distinct notions of the isotropy group associated with the fiducial state must, in the end, define the same group.

In the next section I first review some basic facts about Perelomov coherent states and define the two isotropy groups associated with a fiducial vector. In the third section I show that only if these two isotropy groups coincide does the classical motion provide complete information about the quantum dynamics. Finally I make some brief remarks on the implication of these issues for path integrals.

## II Informative families

Suppose $`G`$ is a compact Lie group and $`|0`$ a vector in an irreducible representation space, $`V`$, of $`G`$. Since $`G`$ is compact, we may assume that $`V`$ comes equipped with an inner product with respect to which the representation is unitary. For $`gG`$ an element of the group we will write its action on a vector $`|vV`$ by $`|vg|v`$. Thus we make no notational distinction between the elements of the group and the corresponding operators $`g:VV`$ in the representation. We will also write $`|g=g|0`$. Following Perelomov, we will call the set $`\{|g\}`$ a family of generalized coherent states. The starting $`|0`$ is the fiducial vector  of the family.

If $`d\mu `$ denotes the invariant measure on the group then

$$B=𝑑\mu |gg|$$
(1)
commutes with all matrices $`g`$, and so, by Shur’s lemma, is proportional to the identity operator. Thus any $`|0`$ provides a resolution of the identity

$$𝐈=const.𝑑\mu |gg|.$$
(2)
This is probably the most important property of such sets of coherent states.

Although any state in the representation produces a family of states with resolution of the identity, the families are not all equivalent, nor are all equally useful. The most commonly seen are those built on highest (or lowest) weight vectors, such as the state $`|j,j`$ in the spin $`j`$ representation of $`SU(2)`$. As with the Schrödinger wave-packets, these families are composed of minimum uncertainty states, and have other nice mathematical properties. In particular they are naturally complex homogeneous manifolds with a Kähler structure. Here, however, we are interested in a broader class of fiducial vectors.

For any $`|0`$ consider the following sets:

* $`H_{|0}=\{hG|h|0=(\mathrm{phase})|0\}`$i.e. the set of elements of $`G`$ for which $`|0`$ is a common eigenvector.
* $`H_0=\{hG|h|\widehat{\lambda }|h=0|\widehat{\lambda }|0\widehat{\lambda }\mathrm{Lie}(G)\}`$i.e. the set of elements of $`G`$ which stabilizes $`f_0=0|\mathrm{}|0`$ considered as an element of $`\left(\mathrm{Lie}(G)\right)^{}`$.

Clearly both sets are subgroups and $`H_{|0}H_0`$.

The first subgroup, $`H_{|0}`$, is the isotropy group of the family of coherent states. In other words, the physically distinct states in the family are in one-to-one correspondence with the quotient space $`G/H_{|0}`$. States in any particular coset differ only by an overall phase.

The second group, $`H_0`$, is the isotropy group of the linear functional $`f_0:\mathrm{Lie}(G)𝐂`$, where $`f_0(\widehat{\lambda })=0|\widehat{\lambda }|0`$, under the co-adjoint action of $`G`$ i.e.

$$f_0(\widehat{\lambda })f_g(\widehat{\lambda })=f_0(g^1\widehat{\lambda }g).$$
(3)
As we will see, the degrees of freedom of the associated classical system live in the co-adjoint orbit $`G/H_0`$.

Since we wish the classical and quantum evolutions to be related, it would be nice if the two coset spaces $`G/H_{|0}`$ and $`G/H_0`$ were the same. Unfortunately the two subgroups $`H_{|0}`$ and $`H_0`$ do not necessarily coincide.

As an example, consider one of the families discussed by Matsumoto. We take as fiducial vector the state

$$|0=\sqrt{\frac{2}{3}}|1,1+\sqrt{\frac{1}{3}}|1,1$$
(4)
in the spin-$`1`$ representation of $`SU(2)`$. We have

$`0|\widehat{J}_1|0`$ $`=`$ $`0`$ (5)
$`0|\widehat{J}_2|0`$ $`=`$ $`0`$ (6)
$`0|\widehat{J}_3|0`$ $`=`$ $`{\displaystyle \frac{1}{3}}.`$ (7)
The isotropy group $`H_0`$ consists of rotations about the “$`3`$” axis. The group $`H_{|0}`$ must be a subgroup of this, but under such rotations

$$\sqrt{\frac{2}{3}}|1,1+\sqrt{\frac{1}{3}}|1,1e^{i\theta }\sqrt{\frac{2}{3}}|1,1+e^{i\theta }\sqrt{\frac{1}{3}}|1,1.$$
(8)
Thus $`H_0=\{e^{i\theta J_3}\}`$, while $`H_{|0}`$ contains only the identity element.

Suppose now that $`H_0`$ is equal to $`H_{|0}`$. In this case the cosets $`G/H_{|0}`$ and $`G/H_0`$ coincide and consequently $`|g`$ is determined (up to a phase) by the values of the expectations $`g|\widehat{\lambda }_i|g=\lambda _i`$ as $`\widehat{\lambda }_i`$ ranges over the Lie algebra of $`G`$. This latter property is a reasonable one to require, since then we can determine (up to a phase) the quantum state of the system from measurement of the these expectations. There may be a standard name for fiducial vectors for which $`H_0=H_{|0}`$, but I am unaware of it. For want of a suitable word, I will call such vectors, and the resulting family of coherent states, informative.

## III Dynamics

If we restrict our attention to Hamiltonians that are elements of the Lie algebra, then a state that starts as a coherent state will remain one under the quantum time evolution.

To associate a classical dynamical process with the quantum one, let $`\widehat{H}\mathrm{Lie}(G)`$ be the quantum Hamiltonian, $`|0`$ be our selected fiducial vector, and define the action

$$S=𝑑t\left\{i0|g^1\dot{g}|00|g^1\widehat{H}g|0\right\}.$$
(9)
The first term is the canonical (or Berry) phase which determines the symplectic structure on the phase space. The second term serves as the classical Hamiltonian.

The variation of $`S`$ is

$$\delta S=𝑑t\left\{0|[ig^1\dot{g}g^1\widehat{H}g,g^1\delta g]|0\right\}.$$
(10)
The corresponding equation of motion for $`g`$ is therefore

$$g^1\dot{g}=ig^1\widehat{H}g+i\widehat{\lambda }(t),$$
(11)
where $`\widehat{\lambda }(t)`$ is any element of the Lie algebra obeying $`0|[i\widehat{\lambda }(t),g^1\delta g]|0=0`$ for all $`g^1\delta g\mathrm{Lie}(G)`$. This condition means that $`\widehat{\lambda }(t)`$ lies in the Lie algebra of $`H_0`$. Since $`g^1\dot{g}`$ is indeterminate up to an an element of $`\mathrm{Lie}(H_0)`$, we must consider the classical trajectories as living in the co-adjoint orbit $`G/H_0`$, rather than in $`G`$. In general, knowing the classical trajectory in this space is not enough to determine the evolution of the quantum state. For example, the evolution of the fiducial state (4) with $`\widehat{H}=\widehat{J}_3`$ alters the relative phase of the two components as in (8). This phase-shift is invisible in the classical motion since the expectation values of the Lie algebra generators do not change.

We can make further progress, however, if we take $`|0`$ to be informative. This means $`\widehat{\lambda }(t)`$ is also in $`H_{|0}`$, so

$$\widehat{\lambda }(t)|0=\lambda (t)|0$$
(12)
for some (real) number $`\lambda (t)`$.

Given this information we note that

$$S=\lambda (t)𝑑t.$$
(13)
Next we observe that the solution of the equation of motion for $`g(t)`$ is

$$g(t)=T\left\{e^{i_0^t\widehat{H}𝑑t}\right\}g(0)\overline{T}\left\{e^{i_0^t\widehat{\lambda }𝑑t}\right\},$$
(14)
where $`\overline{T}`$ denotes anti-time-ordering. Therefore

$$|g(t)=g(t)|0=T\left\{e^{i_0^t\widehat{H}𝑑t}\right\}g(0)e^{i_0^t\lambda (t)𝑑t}|0.$$
(15)
In other words

$$T\left\{e^{i_0^t\widehat{H}𝑑t}\right\}|g(0)=e^{iS}|g(t).$$
(16)
The left hand side of this equation is the quantum time-evolved coherent state, while the right hand side is, up to a phase, the coherent state corresponding to the classically evolved variable $`g(t)`$.

At first sight the action $`S`$ appearing in the above expression is arbitrary. This is because in equation (11) $`\widehat{\lambda }`$ could have been any element of $`\mathrm{Lie}(H_0)`$. The ambiguity is removed however, when we select a specific representative, $`|g`$, from each ray in $`G/H_{|0}`$. This is what is normally done when we define a family of coherent states. For informative states, we have therefore recovered the traditional picture dating back to van Hove. This result, that the quantum evolution may be found by solving a purely classical problem, is essentially equivalent to the Wei-Norman disentangling procedure.

## IV Discussion

Because any fiducial vector gives rise to a resolution of unity, we can use the coherent states constructed on it to write down an exact discrete-time path integral for any transition amplitude. From this, by taking a formal limit of infinitely many intermediate steps, we may “derive” a continuous-time path integral representation of the quantum dynamics. The action appearing in this path integral is (9), together with some boundary terms that serve to make the initial and final value problem well defined. We might reasonably expect the classical paths, those with stationary variations, to play an important role in evaluating this path integral. Unfortunately we have seen that for the general fiducial vector these paths do not capture the full quantum dynamics. We must expect, therefore, substantial analytic difficulties in making rigorous the continuous-time limit of the discrete path integral.

As a symptom of these problems, consider the formal path integral that comes from taking (4) as fiducial vector. The canonical phase term in the classical action is then

$$\frac{1}{3}(\mathrm{cos}\theta 1)\dot{\varphi }𝑑t.$$
(17)
The coefficient, $`1/3`$, violates the condition required to make the “Dirac string” at the south pole invisible. Only integers and half integers are allowed as coefficients if the path integral is to be well-defined.

For the spin coherent states of $`SU(2)`$, it has recently been shown that, when the fiducial vector is taken to be a highest weight vector, the formal semi-classical expansion about the continuous-time classical paths does yield correct answers. It would be a salutary exercise to trace exactly what goes wrong with the continuous-time limit in the more general case.

###### Acknowledgements.

This work was supported by the National Science Foundation under grant DMR-98-17941. I would also like to thank the TCM group at the Cavendish Laboratory, Cambridge, England, for hospitality, and the EPSRC for funding my visit under grant number GR/N00364.


text:
# Magnetostructural effects and phase transition in Cr2O3 under pressure

## Abstract

We have successfully calculated the electronic and structural properties of chromia (Cr2O3) in the Local Spin Density Approximation (LSDA). We predict a transformation from the corundum to the Rh2O3 (II) structure around 15 GPa in the anti-ferromagnetic (AFM) phase as well as in the paramagnetic (PM) insulating state which occurs above the Néel temperature ($`T_N`$). This transition is relevant to interpreting the optical anomalies observed in the absorption spectrum of ruby under pressure. We have modeled the structural properties of the PM state using a Landau-like expansion of the magnetostriction energy. This treatment correctly describes the structural anomalies across $`T_N`$ in the corundum phase and indicates that the AFM and PM insulating states should have distinct compressive behaviors.

Oxides of the $`3d`$ transition metals are a fascinating class of materials with amazingly diverse physical properties. They have been a subject of intensive experimental and theoretical studies for many years . A wide variety of computational techniques — density-functional theory with linear augmented plane wave method , non-periodic embedded cluster approach , periodic unrestricted Hartree-Fock method  — have been used to perform ab initio calculations of these compounds. In this paper we report a successful application of the first principles pseudo-potential plane-wave (PPPW) approach  to compute structural properties of and predict a pressure induced structural transition in chromia ($`Cr_2O_3`$), a typical anti-ferromagnetic insulator in this class of materials. There were a few experimental papers on chromia  and only one ab initio calculation . The periodic unrestricted LCAO Hartree – Fock method was used  to calculate crystal parameters and elastic properties, as well as electronic structure of $`Cr_2O_3`$. Their results are in a good agreement with experimental data. The advantage of the PPPW approach is in the simplicity of the plane-wave basis set, which makes it easy to calculate ionic forces and lattice stresses. This allowed us to optimize dynamically the cell structure of chromia under pressure for various magnetic states.

Our motivation to study $`Cr_2O_3`$ is related to the high pressure behavior of ruby, i.e. $`Al_{2(1x)}Cr_{2x}O_3(x<0.05)`$. The pressure dependence of the fluorescence lines in ruby is widely used to determine pressure (the so-called ruby scale) in diamond-anvil-cell experiments. Alumina ($`Al_2O_3`$) and chromia exist in the corundum phase and form a completely isomorphous alloy system. Alumina has been shown to undergo a structural phase transition to the $`Rh_2O_3`$ (II) phase around 80-100 GPa . A recent theoretical study of the effect of this structural transition on the optical spectrum of ruby indicated that the neighborhood of the chromium site, a color center, might be undergoing a severe distortion around 30 GPa. This hypothesis is suggested by the behavior of the optical absorption lines which display a small discontinuity at 30 GPa and resemble more closely the transitions predicted in the high pressure $`Rh_2O_3`$ (II) phase beyond 30 GPa than those in the corundum phase. It was then anticipated that the cause of this distortion could be a similar phase transformation at lower pressures ($`<`$ 30 GPa) in chromia, the other end member of the alloy. At the moment there is no convincing experimental evidence for a phase transformation in $`Cr_2O_3`$ in this pressure range . However, this may be due to the fact that the $`Rh_2O_3`$ (II) phase has an X-ray diffraction pattern similar to corundum’s , which makes it difficult to observe.

In this paper we investigate this pressure induced transformation in chromia. However another interesting question arises: at room temperature chromia undergoes a change in magnetic phase under pressure. Its Néel temperature is $`T_N=308`$ K , with a pressure dependence of $`T_N/P=16K/GPa`$ . The paramagnetic (PM) state above $`T_N`$ is also insulating; therefore, the effect of the magnetic transition on the structural properties are not expected to be dramatic. Nevertheless, structural anomalies around the Néel temperature are well known  and a realistic prediction of a possible phase transition above 0.5 GPa should be carried out in the PM insulating phase with randomly oriented spins. Here we investigate from first principles this structural transition in the AFM, and in a hypothetical ferromagnetic (FM) phases. The structural properties of PM insulating state is then explored in relation to those of the AFM and FM phases using a phenomenological approach based on a Landau-like expansion of the magnetostriction energy. The predicted structural differences of this phase with respect to the AFM phase correlate well with the anomalies

observed around the Néel transition and are very different from those of the PM metallic phase predicted by a standard LDA calculation.

The crystal structure of $`Cr_2O_3`$ at ambient conditions is corundum-like. It can be described as a hexagonal closed packed array of oxygens with two thirds of the octahedral sites filled with chromium atoms. The unit cell is rhombohedral and contains two formula units. In $`Cr_2O_3`$ each chromium is left with three $`d`$-electrons, losing the other three to oxygens. The predominantly octahedral crystal field splits $`d`$-orbitals into (approximately) a $`t_{2g}`$-like triplet and an $`e_g`$-like doublet. The lower triplet accommodates three electrons. The $`Cr^{3+}`$ local magnetic moment is $`2.76\mu _B`$,  close to the spin only value of $`3\mu _B`$. The moments lie parallel to the $`c`$-axis, and are arranged in ferromagnetic layers perpendicular to the $`c`$-axis. The layers alternate up and down along the $`c`$-axis. In the insulating phase above $`T_N`$ chromia has randomly oriented spins.

To investigate the effect of magnetism on the structural properties we have performed three distinct calculations: 1) a standard spin-polarized LSDA calculation  in the AFM phase of $`Cr_2O_3`$ in the corundum structure; 2) same-type calculation in the FM phase, with the net magnetic moment of $`3\mu _B`$ per $`Cr`$ atom; 3) a standard non-spin-polarized LDA calculation in a PM phase. In all cases the lattice and internal degrees of freedom were dynamically relaxed under pressure . The zero-pressure structures obtained correspond to various local minima of the LSD functional.

The ground state at T = 0 K is the AFM state with a band gap of $``$ 1.5 eV (Fig. 1-a) and a local magnetic moment of $`3\mu _B`$ on chromium atoms (from straight band occupations). The overall band structure compares well with photo-emission data . Namely, the O2p and Cr3d band widths of 5 eV and 1 eV respectively, and O2p-Cr3d band centers separation of 4 eV are in a good agreement with experimental values. However, as expected, the band gap is underestimated with respect to the thermal gap of 3.3 eV . Zero-pressure equilibrium structural parameters presented in Table I are also in good agreement with experimental data. The cohesive energy of Ecoh = 6.1 eV/atom compares well with $`E_{coh}=5.55`$ eV/atom from experiments. The FM state is found to be insulating as well with a band gap of $``$ 0.9 eV (Fig.1-b). After structural relaxation this state is only $`35`$ meV/unit above the AFM ground state. Equilibrium lattice parameters (Table I) are quite different from those in the AFM state, indicating a substantial influence of the magnetic state on the structural properties. The standard paramagnetic non-spin-polarized LDA calculation stabilizes $`Cr_2O_3`$ in a metallic phase (Fig.1-c), as expected, with structural properties considerably different from the observed ones (see Table I). This state is 2.25 eV/unit above the AFM ground state and cannot properly account for the structural transition under consideration.

Our description of the structural properties of the PM phase is based on the Landau-like expansion of the crystal deformation energy at a certain fixed pressure, $`E=\lambda _{ikjl}u_{ik}u_{jl}+\beta _{ikjl}u_{ik}M_{1j}M_{2l}.`$ Here $`u_{ik}`$ is a strain tensor, $`𝐌_1`$, $`𝐌_2`$ are the magnetization vectors of the two AFM sublattices, and $`\lambda `$ and $`\beta `$ are constant tensors. The first term represents pure elastic deformation, while the second describes the magnetostriction energy, i. e., the coupling of the magnetic and structural degrees of freedom. Deformations conserving the corundum structure symmetry allow only for $`U=u_{33}`$, uniaxial strain, and $`V=u_{ii}`$, hydrostatic compression. We assume $`𝐌_1`$, $`𝐌_2`$ remain parallel to the $`z`$-axis. The above expression then simplifies to: $`E=\lambda _1V^2+\lambda _2U^2\lambda _{12}VU+(\beta _1V+\beta _2U)M_{1z}M_{2z}.`$ The equilibrium configurations are such that $`U=u_0M_{1z}M_{2z}`$, $`V=v_0M_{1z}M_{2z}`$, where $`u_0`$ and $`v_0`$ depend on $`\lambda `$ and $`\beta `$. For the AFM and FM phases $`M_{1z}M_{2z}=1`$ and $`M_{1z}M_{2z}=1`$ respectively. For the PM state the average product $`M_{1z}M_{2z}_{PM}`$ = 0. Therefore, the equilibrium lattice structural parameters of the PM phase are given by: $`U_{PM}=(U_{AFM}+U_{FM})/2`$, $`V_{PM}=(V_{AFM}+V_{FM})/2`$ where $`U_{AFM}`$, $`U_{FM}`$, $`V_{AFM}`$, and $`V_{FM}`$ have been determined from first principles. A similar procedure can be adopted for dealing with the Rh2O3 (II) phase. In this case, the deformation energy is expressed as: $`E=\lambda _{\gamma \delta }u_\gamma u_\delta +\beta _\gamma u_\gamma M_{1z}M_{2z}`$, where $`\gamma ,\delta =xx,yy,zz`$.

The predicted properties of the AFM and PM insulating states compare as follows: a) the difference in zero pressure lattice parameters (in the hexagonal cell description) between them, are similar to the anomalies observed around $`T_N`$. Throughout the Néel transition (AFM to PM), the calculated $`\mathrm{\Delta }a_H=+0.013\AA `$ and $`\mathrm{\Delta }c_H=0.11\AA `$ agree in sign and approximately in order of magnitude with experimental values of $`\mathrm{\Delta }a_H=+0.006\AA `$ and $`\mathrm{\Delta }c_H=0.018\AA `$ . We believe this is evidence of the satisfactory description of the PM insulating state. b) The calculated compressive behavior of our PM insulating state compares well with the experimental behavior under pressure above TN , while our AFM calculation is in better agreement with the experiment of Lewis and Drickamer, which we suspect may have been carried out in the AFM phase.

A summary of experimental data and our results is presented in Fig. 2. In Ref.  a substantial decrease of the rhombohedral angle with pressure was found (hexagonal $`c`$-axis less compressive than $`ab`$ axes). In contrast, in Ref.  a slight increase in this angle was observed ($`c`$-axis more compressive than $`ab`$ axes). This discrepancy has been attributed to non-hydrostatic stresses in Ref. . Our results suggest this discrepancy could be real if somehow in Ref.  chromia was kept in the AFM state. We predict here that the magnetic state affects noticeably the compressive behavior of chromia despite the insulating nature of both phases. Compression experiments well below and well above $`T_N`$ could help to clarify this situation.

Now we deal with the structural transition. The high pressure $`Rh_2O_3`$ (II) phase has $`Pbna`$ space group with an orthorhombic unit cell containing 20 atoms (four Cr2O3 units). This phase is structurally similar to the corundum structure and may be described as having a different stacking of similar octahedral layers with a periodicity which is twice that of the corundum along the ’$`c`$’ direction. In the corundum structure the CrO6 octahedra share three edges while in the Rh2O3 (II) structure they share only two. The relative stability of these structures in both magnetic states is shown in Fig. 3-a. We predict the corundum to $`Rh_2O_3`$ (II) transformation in chromia to take place at 14 GPa and 16 GPa in the AFM and in the PM insulating phases respectively, with fractional volume changes at the transition of $`2\%`$. This transition would also take place in the FM phase if it were somehow stabilized.

Fig. 3-b displays the average radius of the first coordination shell around chromium in AFM chromia and in ruby (from Ref. ). The average Cr-O bond-length in chromia increases across the corundum to Rh2O3 (II) transition by an amount similar to that required to explain the optical anomalies in ruby under pressure, which can be explained by a decrease in crystal field splitting. This verifies that this presumable rearrangement could arise from a preference of chromia for the Rh2O3 (II) phase above $``$ 15 GPa. The structural constraint imposed by the alumina host structure should naturally hinder the atomic rearrangement around the color centers until higher pressures, for instance 30 GPa.

These results should stimulate further experimental and theoretical work. The prediction of distinct compressive behaviors in chromia above and below TN and the structural phase transformation near 15 GPa await experimental confirmations. The latter, if verified, makes ruby an interesting study case: an isomorphous alloy in which both end members undergo the same structural transition but at very different pressures. Intermediate compositions should undergo similar transitions at intermediate pressures. However, before the transformation manifests macroscopically it could be nucleating around one of the components, even in the impurity limit. The possibility of investigating this phenomenon in ruby by EXAFS or anomalous X-ray scattering is fascinating.

Acknowledgments \- We thank Phil Allen and Stefano Baroni for helpful comments on the manuscript. This work was supported by the University of Minnesota MRSEC (AYuD), and NSF grant EAR-9973130 (RMW and WD).


text:
# Time series model based on global structure of complete genome

## 1 Introduction

The nucleotide sequences stored in GenBank have exceeded hundreds of millions of bases and they increase by ten times every five years. A great deal of information concerning origin of life, evolution of species, development of individuals, and expression and regulation of genes, exist in these sequences$`^{\text{[1]}}`$. In the past decade or so there has been an enormous interest in unravelling the mysteries of DNA. It has become very important to improve on new theoretical methods to do DNA sequence analysis. Statistical analysis of DNA sequences<sup>\[1-9\]</sup> using modern statistical measures is proven to be particularly fruitful. There is another approach to research DNA, namely nonlinear scales method, such as fractal dimension$`^{\text{[10, 11, 12, 13]}}`$, complexity$`^{\text{[14, 16]}}`$. The correlation properties of coding and noncoding DNA sequences was first studied by Stanley and coworkers$`^{\text{[5]}}`$ in their “fractal landscape or DNA walk” model. The DNA walk defined in  is that the walker steps “up” if a pyrimidine ($`C`$ or $`T`$) occurs at position $`i`$ along the DNA chain, while the walker steps “down” if a purine ($`A`$ or $`G`$) occurs at position $`i`$. Stanley and coworkers$`^{\text{[5]}}`$ discovered there exists long-range correlation in noncoding DNA sequences while the coding sequences correspond to regular random walk. But if one considers more details by distinguishing $`C`$ from $`T`$ in pyrimidine, and $`A`$ from $`G`$ in purine (such as two or three dimensional DNA walk model$`^{\text{[1]}}`$ and maps given in ), then the presence of base correlation has been found even in coding region. However, DNA sequences are more complicated than those these types of analysis can describe. Therefore, it is crucial to develop new tools for analysis with a view toward uncovering mechanisms used to code other types of information.

Since the first complete genome of the free-living bacterium Mycoplasma genitalium was sequenced in 1995$`^{\text{[17]}}`$, an ever-growing number of complete genomes has been deposited in public databases. The availability of complete genomes opens the possibility to ask some global questions on these sequences. The avoided and under-represented strings in some bacterial complete genomes have been discussed in . A time series model of CDS in complete genome has also been proposed in .

One can ignore the composition of the four kind of bases in coding and noncoding segments and only consider the roughly structure of the complete genome or long DNA sequences. Provata and Almirantis $`^{\text{[20]}}`$ proposed a fractal Cantor pattern of DNA. They map coding segments to filled regions and noncoding segments to empty regions of random Cantor set and then calculate the fractal dimension of the random fractal set. They found that the coding/noncoding partition in DNA sequences of lower organisms is homogeneous-like, while in the higher eucariotes the partition is fractal. This result is interesting and reasonable, but it seems too rough to distinguish bacteria because the fractal dimensions of bacteria they gave out are all the same. The classification and evolution relationship of bacteria is one of the most important problem in DNA research. In this paper, we propose a time series model based on the global structure of the complete genome and we find that one can get more information from this model than that of the fractal Cantor pattern. We have found some new results to the problem of classification and evolution relationship of bacteria.

A DNA sequence is a sequence over the alphabet $`\{A,C,G,T\}`$ representing the four bases from which DNA is assembled, namely adenine, cytosine, guanine, and thymine. But from views of the level of structure, the complete genome of organism is made up of coding and noncoding segments. Here the length of a coding/noncoding segment means the number of its bases. First we simply count out the lengths of coding/noncoding segments in the complete genome. Then we can get three kinds of integer sequences by the following ways.

i) First we order all lengths of coding and noncoding segments according to the order of coding and noncoding segments in the complete genome, then replace the lengths of noncoding segments by their negative numbers. So that we can distinguish lengths of coding and noncoding segments. This integer sequence is named whole length sequence.

ii) We order all lengths of coding segments according to the order of coding segments in the complete genome. We name this integer sequence coding length sequence. For some examples, we plot the distribution of coding length sequences of three bacteria genome and the 4th chromosome of Saccharomyces cerevisiae (yeast) in Figure 1.

iii) We order all lengths of noncoding segments according to the order of noncoding segments in the complete genome. This integer sequence is named noncoding length sequence.

We can now view these three kinds of integer sequences as time series. We want to calculate their correlation dimensions and Hurst exponents.

## 2 Correlation dimension and Hurst exponent

The notion of correlation dimension, introduced by Grassberger and Procaccia$`^{\text{[21, 22]}}`$, suits well experimental situations, when only a single time series is available. It is now being used widely in many branches of physical science. Consider a sequence of data from a computer or laboratory experiment:

$$x_1,x_2,x_3,\mathrm{},x_N,$$
(1)
where $`N`$ is a large enough number. These numbers are usually sampled at an equal time interval $`\mathrm{\Delta }\tau `$. We embed the time series into $`𝐑^m`$, choose a time delay $`\tau =p\mathrm{\Delta }\tau `$, then obtain

$$𝐲_i=(x_i,x_{i+p},x_{i+2p},\mathrm{},x_{i+(m1)p}),i=1,2,\mathrm{},N_m,$$
(2)
where

$$N_m=N(m1)p.$$
(3)
In this way we get $`N_m`$ vectors in the embedding space $`𝐑^m`$.

For any $`𝐲_i,𝐲_j`$, we define the distance as

$$r_{ij}=d(𝐲_i,𝐲_j)=\underset{l=0}{\overset{m1}{}}|x_{i+lp}x_{j+lp}|.$$
(4)
If the distance is less than a given number $`r`$, we say that these two vectors are correlated. The correlation integral is defined as

$$C_m(r)=\frac{1}{N_m^2}\underset{i,j=1}{\overset{N_m}{}}H(rr_{ij}),$$
(5)
where $`H`$ is the Heaviside function

$$H(x)=\{\begin{array}{c}1,\mathrm{if}x>0,\hfill \\ 0,\mathrm{if}x0.\hfill \end{array}$$
(6)
For a proper choice of $`m`$ and not too big a value of $`r`$, it has been shown by Grassberger and Procaccia$`^{\text{[22]}}`$ that the correlation integral $`C_m(r)`$ behaves like

$$C_m(r)r^{D_2(m)}.$$
(7)
Thus one can define the correlation dimension as

$$D_2=\underset{m\mathrm{}}{lim}D_2(m)=\underset{m\mathrm{}}{lim}\underset{r0}{lim}\frac{\mathrm{ln}C_m(r)}{\mathrm{ln}r}.$$
(8)
For more details on $`D_2`$, the reader can refer to .

To deal with practical problems, one usually choose $`p=1`$. If we choose a sequence $`\{r_i:1in\}`$ such that $`r_1<r_2<r_3<\mathrm{}<r_n`$, then a scaling region can be found in the $`\mathrm{ln}r\mathrm{ln}C_m(r)`$ plane, see , p.346. Then the slop of the scaling region is $`D_2(m)`$. When $`D_2(m)`$ does not change with $`m`$ increasing, we can take this $`D_2(m_0)`$ as the estimate value of $`D_2`$. We calculate the correlation dimensions of three kinds of length sequences of the complete genome using the method introduced above. From the $`\mathrm{ln}r\mathrm{ln}C_m(r)`$ figures of these sequences of different values of embedding dimension $`m`$, we find that it is suitable to choose $`m=7`$. For example, we give the $`\mathrm{ln}r\mathrm{ln}C_m(r)`$ figure of whole length sequence of A. fulgidus when $`m=6,7`$ (Figure 2). We take the region from the third point to the 20th point (from left to right) as the scaling region.

Hurst$`^{\text{[24]}}`$ invented the now famous statistical method — the rescaled range analysis ($`R/S`$ analysis) to study the long-range dependence in time series. Later on, B. B. Mandelbrot$`^{\text{[25]}}`$ and J. Feder $`^{\text{[26]}}`$ brought $`R/S`$ analysis into fractal analysis. For any time series $`x=\{x_k\}_{k=1}^N`$ and any $`2nN`$, one can define

$$<x>_n=\frac{1}{n}\underset{i=1}{\overset{n}{}}x_i$$
(9)
$$X(i,n)=\underset{u=1}{\overset{i}{}}[x_u<x>_n]$$
(10)
$$R(n)=\underset{1in}{\mathrm{max}}X(i,n)\underset{1in}{\mathrm{min}}X(i,n)$$
(11)
$$S(n)=[\frac{1}{n}\underset{i=1}{\overset{n}{}}(x_i<x>_n)^2]^{1/2}.$$
(12)
Hurst found that

$$R(n)/S(n)(\frac{n}{2})^H.$$
(13)
$`H`$ is called the Hurst exponent.

As $`n`$ changes from 2 to $`N`$, we obtain $`N1`$ points in the $`\mathrm{ln}(n)`$ v.s. $`\mathrm{ln}(R(n)/S(n))`$ plane. Then we can calculate the Hurst exponent $`H`$ of the length sequence of organisms using the least-squares linear fit. As an example, we plot the graph of $`R/S`$ analysis of the whole length sequence of A. fulgidus in Figure 3.

The Hurst exponent is usually used as a measure of complexity. The trajectory of the record is a curve with fractal dimension $`D=2H`$ (,p.149). Hence a smaller $`H`$ means a more complex system. When applied to fractional Brownian motion, the system is said to be persistent if $`H>1/2`$, which means that if for a given time period $`t`$, the motion is along one direction, then in a succeeding time, it is more likely that the motion will follow the same direction. For $`H<1/2`$, the opposite holds, that is, the system is antipersistent. But when $`H=1/2`$, the system is a Brownian motion, and is random.

## 3  Data and results.

More than 21 bacterial complete genomes are now available in public databases . There are five Archaebacteria: Archaeoglobus fulgidus, Pyrococcus abyssi, Methanococcus jannaschii, Aeropyrum pernix and Methanobacterium thermoautotrophicum; four Gram-positive Eubacteria: Mycobacterium tuberculosis, Mycoplasma pneumoniae, Mycoplasma genitalium, and Bacillus subtilis. The others are Gram-negative Eubacteria. These consist of two Hyperthermophilic bacteria: Aquifex aeolicus and Thermotoga maritima; six proteobacteria: Rhizobium sp. NGR234, Escherichia coli, Haemophilus influenzae, Helicobacter pylori J99, Helicobacter pylori 26695 and Rockettsia prowazekii; two chlamydia Chlamydia trachomatis and Chlamydia pneumoniae, and two Spirochete: Borrelia burgdorferi and Treponema pallidum.

We calculate the correlation dimensions and Hurst exponents of three kinds of length sequences of the above 21 bacteria. The estimated results are given in Table 1 ( we denote by $`D_{2,whole}`$, $`D_{2,cod}`$ and $`D_{2,noncod}`$ the correlation dimensions of whole, coding and noncoding length sequences, from top to bottom, in the increasing order of the value of $`D_{2,cod}`$ ) and Table 2 ( we denote by $`H_{whole}`$, $`H_{cod}`$ and $`H_{noncod}`$ the Hurst exponents of whole, coding and noncoding length sequences, from top to bottom, in the increasing order of the value of $`H_{cod}`$ ).

## 4 Discussion and conclusions

Although the existence of the archaebacterial urkingdom has been accepted by many biologists, the classification of bacteria is still a matter of controversy$`^{\text{[27]}}`$. The evolutionary relationship of the three primary kingdoms (i.e. archeabacteria, eubacteria and eukaryote) is another crucial problem that remains unresolved$`^{\text{[27]}}`$.

From Table 1, we can roughly divide bacteria into two classes, one class with $`D_{2,cod}`$ less than 2.40, and the other with $`D_{2,cod}`$ greater than 2.40. We observe that the classification of bacteria using $`D_{2,cod}`$ almost coincides with the traditional classification of bacteria. All Archaebacteria belong to the same class. All Proteobacteria belong to the same class except Rhizobium sp. NGR234, in particular, the closest Proteobacteria Helicobacter pylori 26695 and Helicobacter pylori J99 group with each other. Two Spirochete group with each other. Two Chlamydia gather with each other. Gram-positive bacteria is divided into two sub-categories: Mycoplasma genitalium and Mycoplasma pneumoniae belong to one class and gather with each other, Mycobacterium tuberculosis and Bacillus subtilis belong to another class and almost gather with each other.

If one classifies bacteria using $`D_{2,whole}`$, with the $`D_{2,whole}`$ of one subclass less than 2.55, that of the other larger than 2.55, almost the same results hold as those using $`D_{2,cod}`$. But when one classifies bacteria using $`D_{2,noncod}`$, the results are quite different. This is quite reasonable because the coding segments occupy the main part of space of the DNA chain of bacteria.

A surprising feature shown in Table 1 is that the Hyperthermophilic bacteria (including Aquifex aeolicus and Thermotoga maritima) are linked closely with the Archaebacteria if we only consider the length sequences of noncoding segments. But when we consider the length sequences of coding segments, they are linked closely with eubacteria. We notice that Aquifex, like most Archaebacteria, is hyperthermophilic. Hence it seems that their hyperthermophilicity property is possibly controlled by the noncoding part of the genome, contrary to the traditional view resulting from classification based on the coding part of the genome. It has previously been shown that Aquifex has close relationship with Archaebacteria from the gene comparison of an enzyme needed for the synthesis of the amino acid trytophan$`^{\text{[28]}}`$. Such strong correlation on the level of complete genome between Aquifex and Archaebacteria is not easily accounted for by lateral transfer and other accidental events$`^{\text{[28]}}`$. Our result is based on different levels of the genome from that used by the authors of

From Table 1, one can also see the $`D_{2,cod}`$ values are almost larger than the $`D_{2,noncod}`$ values. Hence the coding length sequences are more complex than the noncoding length sequences.

From Table 2, we can also roughly divide bacteria into two classes, one class with $`H_{cod}`$ less than 0.60, and the other with $`H_{cod}`$ greater than 0.60. One can see all Archeabacteria belong to the same class except Pyrococcus abyssi. All Gram-positive Eubacteria belong to the same class except Bacillus subtilis. All Proteobacteria belong to the same class except Haemophilus influenzae. Two Spirochete group with each other. Two Chlamydia almost group with each other.

We also find the $`H_{noncod}`$ values of all Archeabacteria except Pyrococcus abyssi, two Hyperthermophilic bacteria, and Mycoplasma genitalium and E. coli are less than $`1/2`$, while those of other bacteria are greater than $`1/2`$. Hence Hyperthermophilic bacteria have some common information with Archaebacteria in noncoding segments.

We calculate $`D_{2,whole}`$, $`D_{2,cod}`$, $`D_{2,noncod}`$ $`H_{whole}`$, $`H_{cod}`$ and $`H_{noncod}`$ of the 4th chromosome of Saccharomyces cerevisiae (yeast). They are 2.5603, 2.1064, 2.5013, 0.5517, 0.6255 and 0.5482 respectively. From Tables 1 and 2, if we consider $`D_{2,whole}`$, $`H_{whole}`$ and $`H_{cod}`$, we can see that Archaebacteria and Chlamydia are linked more closely with yeast which belongs to eukaryote than other categories of bacteria. There are several reports (such as ) that, in some RNA and protein species, archeabacteria are much more similar in sequences to eukaryotes than to eubacteria. Our present result supports this point of view.

In , we find that the Hurst exponent is a good tool to distinguish different functional regions. But now considering more global structure of the genome, we find the correlation dimension a better exponent to use for classification of bacteria than the Hurst exponent in this level.

## ACKNOWLEDGEMENTS

One of the authors Zu-Guo Yu would like to express his thanks to Prof. Bai-lin Hao of Institute of Theoretical Physics of Chinese Academy of Science for introducing him into this field and continuous encouragement. He also wants to thank Dr. Bin Wang of ITP and Dr. Fawang Liu at QUT for useful discussions about computer programmes. This project is supported by Postdoctoral Research Support Grant No. 9900658 of QUT.