text:
# Int. J. Bifurcation and Chaos 5, 51–61 (1995) Entropy, Transinformation and Word Distribution of Information–Carrying Sequences

## I Introduction

In the manifold of structures which are used as information carriers in nature, culture and engineering, linear strings consisting of sequences of letters play a central role. This may be demonstrated by the following examples: The proteins and polynucleotids, the main information carriers in living systems are sequences of amino acids and/or nucleotides . Further most of the messages transporting information between human informational systems have the form of strings of letters. Examples are books or letters, music, computer programs etc. By using the methods of symbolic dynamics any trajectory of a dynamic system may be mapped to a string of letters on certain alphabet . Hence each sequence can be interpreted as the trajectory of a discrete dynamic system.

This work is devoted to the investigation of strings of the type characterized above, i.e. to objects (documents, programs etc.) which may be mapped to strings of letters. The main aim of the investigation is the analysis of long range correlations in information carrying strings. For several reasons we expect the existence of long range structures in these sequences . Especially we expect correlations which range from the beginning of a string up to its end. Let us discuss some of the reasons for this behavior:

1. Predictability: We know from our every-day experience and from scientific research that the identification of the first hundred or thousand letters of the string tells us already a lot about the continuation. Often we make a decision in a book shop after reading just one page. For example, if we find there several times the word love and tennis we expect to find them on the other pages again and again, but if we find first words like file and program we expect to remain in quite another field. Listening to the beginning of a Bach praeludium where the general theme is worked out we expect to hear it again in many variations up to the very end. Such expectations are only justified if there exist indeed long range correlations. This is the scientific expression of our intuitive expectations which are based on the experience that texts and music have certain inherent predictability.
2. Syntactical limitations: Another heuristic reason to expect long range correlations is the exponential explosion of the variety of possible subwords with increasing length for uncorrelated strings. Uncorrelated sequences generated on an alphabet of $`\lambda `$ letters have a variety of

$$N(n)=\lambda ^n=\mathrm{exp}\left(\mathrm{ln}\lambda n\right)$$
(1)
different subwords of length $`n`$. A subword is here any combination of letters including the space, punctuation marks etc. Strings of this type may be generated by stochastic processes of Bernoulli–type or by fully developed chaotic dynamics using alphabets of $`\lambda `$ letters. For $`n>100`$ the number $`N(\lambda )`$ is extremely large. In other words we need very sharp restrictions to select a meaningful subset out of it. Long range correlations provide such a selection criterion. Denoting the selected subset by $`N^{}(n)`$ we expect

$$\frac{N^{}(n)}{N(n)}0\mathrm{if}n\mathrm{}$$
(2)
Bounds of this kind are given by syntax and semantics. The syntactical rules do not allow for an arbitrary concatenation of words to sentences, most of them are forbidden. Furthermore we know that texts (and pieces of music) are formed by keywords (motifs) which are the raw material for the generation of a text (a piece of music). In fact all these rules lead to slower growth with $`n`$.

A rather sharp restriction on the growth with $`n`$ corresponds to the power law

$$N^{}(n)n^\alpha .$$
(3)
Symbolic strings generated by the logistic map in the Feigenbaum–point have this scaling property . The conjecture that several natural objects have this type of scaling has been made earlier . We will show here that pieces of music possibly belong to this class. Another growth law which is much faster is the stretched exponential law

$$N^{}(n)\mathrm{exp}\left(Cn^\alpha \right)\mathrm{with}\alpha <1.$$
(4)
This scaling which is typically for intermittent processes was observed for several texts . We mention that in the limit $`\alpha 0`$ this law corresponds to the scaling in eq. (3). The reduction due to the scaling rule (4) is not as strong as in the scaling rule (3), however, the reduction factor corresponding to rule (4) is still enormous for large $`n`$.
3. Evolution: The third reason to expect such behavior is the rather general idea that evolution operates in regions where long relaxation times, $`1/f`$-noise and other long range correlations are essential .

## II Entropy–like Measures of Sequence Structure

In the following section we will study entropy-like quantities as a measure of the long range correlations in sequences. In order to describe the structure of a given string of length $`L`$ using an alphabet of $`\lambda `$ letters $`\{A_1A_2\mathrm{}A_\lambda \}`$ we introduce the following notations :

Let $`A_1A_2\mathrm{}A_n`$ be the letters of a given substring of length $`nL`$. Furthermore let

$$p^{(n)}(A_1\mathrm{}A_n)$$
(5)
be the probability to find in a string a block of length $`n`$ (subword of length $`n`$) with the letters $`A_1\mathrm{}A_n`$. The probability of having a pair with $`(n2)`$ arbitrary letters in between we denote by

$$p^{(n)}(A_1,A_n)=p^{(n)}(A_1\mathrm{?}\mathrm{?}\mathrm{?}\mathrm{?}A_n).$$
(6)
We introduce the following quantities:

1. the mutual information, also called transinformation

$$I(n)=\underset{A_iA_j}{}p^{(n)}(A_i,A_j)\mathrm{log}\left(\frac{p^{(n)}(A_i,A_j)}{p^{(1)}(A_i)p^{(1)}(A_j)}\right),$$
(7)
2. the entropy per subword of length $`n`$

$$H_n=p^{(n)}(A_1\mathrm{}A_n)\mathrm{log}p^{(n)}(A_1\mathrm{}A_n),$$
(8)
3. the uncertainty of the letter following a block of length $`n`$

$$h_n=H_{n+1}H_n,$$
(9)
4. the entropy of the source (related to the Kolmogorov–Sinai entropy)

$$h=\underset{n\mathrm{}}{lim}h_n.$$
(10)

In an earlier paper  we assumed the scaling behavior

$$\begin{array}{c}H_n=nh+gn^{\mu _0}(\mathrm{log}n)^{\mu _1}+e\hfill \\ 0\mu _0<1\mathrm{or}\mu _0=1,\mu _1<0.\hfill \end{array}$$
(11)
Related assumptions have been made independently by several authors . From that point of view the long range order of strings may be well characterized by the asymptotic for large $`n`$. Chaotic and stochastic strings of the standard type have the property $`h>0`$. Special cases are Bernoulli processes or fully developed chaos with

$$\begin{array}{c}p^{(n)}(A_1\mathrm{}A_n)=a/\lambda ^n\hfill \\ H_n=n\mathrm{log}\lambda \hfill \\ h_n=\mathrm{log}\lambda .\hfill \end{array}$$
(12)
In the following we write all entropy measures in units of $`\mathrm{log}\lambda `$. For Markov processes with memory $`m`$ the uncertainty decreases during the first $`m`$ steps and remains then constant

$$h_{m+k}=h,k>0.$$
(13)
On the other hand any string with a finite period $`p`$ is characterized by

$$H_{p+k}=\mathrm{const}.,h_{p+k}=0\mathrm{if}k>0.$$
(14)
This means that any new letter added to a string does not increase the complexity of the sequence. Consequently we find for periodic strings $`h=0`$ and $`g=0`$. Of special interest to our further considerations are systems which are neither periodic nor chaotic and which are characterized by

$$g>0,h1.$$
(15)
Presently there are not enough data available to estimate the limit entropy $`h`$ for homogeneous texts (written by one author). We follow in this respect the seminal work by Claude E. Shannon who concluded in his pioneering paper: “From this analysis it appears that, in ordinary literary English, the long range statistical effects (up to 100 letters) reduce the entropy…”. Shannon gave an estimate of 40 bits for the entropy of $`n=100`$ letters, i.e. about 0.4 bit per letter. Transforming the bits to $`\lambda `$–units (which we use throughout this paper) we get $`H_{100}=8`$ and we find for the entropy per letter the value $`0.08`$. According to several general investigations  it is not likely that the uncertainty (entropy per letter) decreases still further for larger $`n`$n. Based on these considerations we assume in the following that the limit entropy (in $`\lambda `$–units) is in the region

$$0.01h0.1.$$
(16)
Since a reliable estimate seems to be impossible at present we simply neglect the contribution $`nh`$ to the block entropy in eq. (11). In this way we obtain for the intermediate region $`1n30`$ the formula which will be the basis for our further investigations

$$H_n=gn^{\mu _0}(\mathrm{ln}n)^{\mu _1}+e.$$
(17)
Special cases are:

* logarithmic tails

$$H_n=g(\mathrm{ln}n)^{\mu _1}+e,\mu _1<0$$
(18)
* power law tails

$$H_n=gn^{\mu _0}+e,0<\mu _0<1$$
(19)

The working hypothesis developed earlier  is that strings characterized by eqs. (18) or (19) being on the borderline between order and chaos might be prototypes of information carrying sequences.

Following a relation derived by McMillan and Khinchin the $`n`$–letter entropy and the mean number of subwords of length $`n`$ are related

$$N_n^{}=\lambda ^{H_n}.$$
(20)
In this way a logarithmic entropy scaling corresponds to a power law of the numbers of subwords and a power law scaling of the entropy corresponds to a stretched exponential growth of the number of subwords.

The mutual information (transinformation) defined by eq. (7) is not monotonically decreasing with increasing $`n`$. We define here long range tails as any non exponential decay or increase of the averaged transinformation $`I(n)`$. Periodic strings show correlations of infinite long range. Periodicity with the period $`p`$ implies that all conditional probabilities $`p^{(n)}(A_iA_j)/p^{(1)}(A_i)`$ for $`n>0`$ are also periodic. This leads to a periodic behavior of the transinformation. Therefore sequences with long range correlations show a fluctuating behavior at large scales .

## III Mutual Information and Word Distributions for Finite Information–Carrying Strings

Printed texts in natural languages and music written in the language of notes are examples of information–carrying strings. Other examples, which we do not consider here, are biosequences, where some evidence for the existence of long range correlations was found . Originally texts or pieces of music were generated by the writer or composer as a dynamical process in real time. Today we find in books or albums the frozen in results of this process in the form of a symbolic sequence. Certainly texts or pieces of music are symbolic sequences of high complexity. In contrast to other dynamical processes, writing and composing have developed during a long way of evolution being intended for communication between human beings. In spite of all these difficulties let us now follow the largely simplifying assumption due to Shannon, Mandelbrot and others that texts and pieces of music may be considered as the results of a stationary random process. Although this assumption is still controversial we will take it here as a basis for the further analysis. In our analysis we considered the following sequences:

1. Sonata for piano forte op. 31 No 2 by L. v. Beethoven ($`L4,920`$)
2. Moby Dick by H. Melville ($`L1,170,200`$)
3. Grimm’s Fairy Tales by the Brothers Grimm ($`L1,435,800`$)

Furthermore a few comparisons were made with the Praeludium in F–Major by J. S. Bach and with the sonata KV 311 by W. A. Mozart. In the case of texts we used an alphabet consisting of the 32 symbols

a b c d e f $`\mathrm{}`$ x y z , . ( ) # $``$

The last symbol $``$ stands for the empty space and # for any number. In the case of music we encoded the notes for $`2\frac{1}{2}`$ octaves beginning with the low $`A`$ and ending with the high $`D`$ on an alphabet with 32 symbols. The white keys on the piano forte where encoded by the 18 symbols

A H C D E F G a h c d e f g m o p r

and the black keys beginning with the lower Be and ending with the high Cis were encoded by the 12 symbols

B I J K L b i j k l n q

The pause was encoded by the score “\_” and holding a tone by the “–”. In order to get a better statistics we also used compressed alphabets consisting of 3 letters $`O`$, $`M`$, $`L`$ only. The letter $`O`$ codes for vowels in the case of texts or for a move down in the case of music, the letter $`M`$ stands for consonants or move ups, the letter $`L`$ codes for all other symbols, e.g. pauses (spaces), holding the tone and punctuation marks. For the analysis of the mutual information we have to count here the frequencies of pairs. Since the number of different pairs is $`32^2=1024`$ we have a rather good statistics if the length $`L`$ of the string satisfies the inequality $`L1024`$. As shown by several authors  the transinformation is a reliable measure of the correlations of letters at the distance $`n`$. Every peak at $`n`$ in the transinformation corresponds to a strong positive correlation. In fig. 1 we show the pair correlations of the Beethoven sonata and of music by Bach and Mozart . The peaks show that there exist strong correlations between two notes at characteristic distances. The interpretation of these peaks must be left to specialists. We further notice some similarity in the correlation structure of Bach’s and Beethoven’s music and a distinct different structure of Mozart’s music. A more detailed study of the differences between composers will be given in a separate paper .

In Figs. 2 and 3 the mutual information calculated for Moby Dick and for Grimm’s Tales ($`\lambda =32`$) are drawn. The results show that there are well expressed correlations in the range $`n=1\mathrm{}25`$ which are followed by a long slowly decaying tail. The results for the transinformation $`I(k)`$ become meaningless if the values are smaller than the level of natural fluctuations due to the finite length $`L`$ of the text which is

$$\delta I=\frac{\lambda ^22\lambda }{2L\mathrm{ln}\lambda }.$$
(21)
Although the fluctuation level decays with $`1/L`$ it has even for the rather long text Moby Dick with $`L=1,170,200`$ a value of about 0.00012. This means, as seen from Fig. 3 that any conclusions suggesting long range correlations beyond $`n=300`$ are rather problematic. However, the range of studies of $`I(k)`$ may still be extended by using length corrections . An alternative method is based on studies of the dependence of $`I(k)`$ on $`1/L`$ (which presumably is linear) and by extrapolations $`(1/L)\mathrm{}`$ .

As we see from Figs. 2 and 3 the mutual information decreases monotonously up to $`n300`$ and converges into the fluctuation level. There are no well expressed correlation peaks. Evidently long texts are in this respect less correlated than DNA sequences where well expressed long range pair correlations have been found .

The results obtained so far base only on the statistical distributions of pairs. In this case one can reach a rather good statistics by counting the probability of pairs. Let us study now the distribution of words of length $`n>2`$. Due to the fact that the number of different words of length $`n`$ using an alphabet consisting of 32 letters is

$$N_n^{}32^n=2^{5n}$$
(22)
there are for $`n=9`$ already $`2^{45}`$ words with approximately $`2^{50}`$ letters. This is much more than all texts stored in libraries which have been estimated to consist of about $`10^{15}2^{18}`$ letters . Therefore we could not do a statistic analysis if there were no additional constraints due to grammar and semantics. According to an earlier estimate  we expect that the growth law scales as

$$N_n^{}\mathrm{exp}\left(An^\alpha \right)$$
(23)
with $`\alpha \frac{1}{2}`$ for texts.

We have done the analysis with two long texts, namely Moby Dick and Grimm’s Tales. Figs. 4 and 5 show the rank ordered distribution of subwords of length $`n=4,9,16,25`$. The structures of both distributions are similar but the lists of words are quite different. For example among the most frequent subwords of length $`n=16`$ in the case of Moby Dick are “the\_sperm\_whale\_”, and “the\_quarter\_deck”. For Grimm’s Tales rather frequent subwords of length $`n=25`$ are e.g. “\_if\_i\_could\_but\_shudder.\_” and “\_princess,\_open\_the\_door\_”.

The shapes of the subword distributions are distinctly not Zipf–like, they do not follow a power law. The distribution tends to form a Fermi–like plateau. This tendency is based on the theorem of asymptotic equipartition derived by McMillan  and Khinchin . This theorem proves that for $`n\mathrm{}`$ the asymptotic form of the distribution is rectangular, i.e.

$$Z(j)=\{\begin{array}{cc}\hfill 1/N^{}(n)& ifi<N^{}(n)\hfill \\ \hfill 0& else.\hfill \end{array}$$
(24)
The effects due to finite $`n`$ and to the finite text length $`L`$ tend to smooth the edges of its distribution. Since our texts are rather long the finite size effects do not have too much influence to the distributions. Much more difficult is the analysis of relative short pieces of music. Here additional problems arise due to the short sample. For example Beethoven’s sonata consists of only 4,920 letters (notes). The importance of length corrections for estimating the frequencies of words was considered by several authors . For a deeper analysis of this problem we refer to . The method we used here was found by generalizing a method proposed in . According to this method the unknown distribution function of the words is guessed by a comparison of expected (generated) and observed distributions. Instead of the simple rectangular distributions in eq. (24) we applied a more realistic expression. For given word length $`n`$ we guess the true (not normalized) distribution, i.e. the distribution for $`L\mathrm{}`$, in the form

$$Z(j,n)=Z_0(j,n)+Z_1(j,n)+Z_2(j,n)$$
(25)
with

$`Z_0(j,n)`$ $`=`$ $`{\displaystyle \frac{z_0(n)}{1+\mathrm{exp}\left(b(jj_0(n))\right)}}`$ (26)
$`Z_1(j,n)`$ $`=`$ $`z_1(n)\mathrm{exp}\left({\displaystyle \frac{j1}{j_1(n)}}\right)`$ (27)
$`Z_2(j,n)`$ $`=`$ $`(z_2(n)z_1(n)z_0(n))`$ (29)
$`\mathrm{exp}\left({\displaystyle \frac{((j1)/j_2(n))^2}{2}}\right).`$
This distribution has 7 free parameters, one of them is fixed by the condition that the total number of words is given by the size of the text. For a string of length $`L`$ the total number of $`n`$–words is $`N=Ln+1`$. Each of the parameters has a simple meaning as:
$`z_2(n)`$ - frequency of the most frequent word of length n,  $`z_1(n)`$ - height of the exponential contribution at $`j=1`$,    $`z_0(n)`$ - height of the asymptotic “Fermi”-plateau,  $`j_2(n)`$ - number of frequent words,  $`j_1(n)`$ - number of relatively frequent words,  $`j_0(n)`$ - number of different words,  $`b(n)`$ - reciprocal “Fermi”-temperature of the plateau.

All these parameters characterize certain properties of the assumed generating process. For complex processes such as writing of novels or pieces of music the parameters are of course unknown. Due to their simple meaning, however, it is easy to guess the initial values of the parameters. The next step in an iterative process is then to generate a sample of $`n`$–words (with the help of the guessed probability distribution $`Z(j,n))`$ which has the same number of $`n`$–words as the text and transform it to a rank–ordered distribution. To overcome problems of low statistics we average over several generated distributions. The result of this procedure is called the expected distribution $`Z^{exp}(j,n,L)`$ according to the chosen set of parameters. The parameters are fitted by adaptation of the expected to the empirical (observed) distribution $`Z^{obs}(j,n,L)`$. The procedure is iterated up to a sufficient agreement between the observed and the expected distributions

$$\underset{j=1}{\overset{N}{}}\left[Z^{obs}(j,n,L)Z^{exp}(j,n,L)\right]^2min.$$
(30)
In the minimization procedure we used the simplex–method applying the program package MINUIT  developed for non–linear parameter optimization. For word length $`n=8`$ the optimization parameters are given in Tab. I. The corresponding distributions $`Z(j,n)`$, $`Z^{exp}(j,n,L)`$ and $`Z^{obs}(j,n,L)`$ are presented in figs. (6) and (7).

## IV Analysis of Entropy Scaling

As shown above the set of entropies of subwords is the basic measure for the investigation of the order of a string. The probabilities one need for the calculation of the entropies are in general unknown and have to be estimated from the frequencies of words. The theoretical value of the block entropy follows from the distribution discussed in the last section. For simplicity we will omit in the following the block length $`n`$. The theoretical entropy is then

$$H=\underset{i}{}\frac{Z(i)}{N}\mathrm{log}\left(\frac{Z(i)}{N}\right).$$
(31)
The calculation of the block entropy from the distribution function $`Z(i)`$ is rather difficult and time consuming. The direct observation of the entropy is based on the frequencies of the different words. Replacing the probabilities by the observed frequencies $`N_i/N`$ in the entropy definition we obtain the observed entropy

$$H^{obs}=\mathrm{log}(N)\frac{1}{N}\underset{i}{}N_i\mathrm{log}\left(N_i\right).$$
(32)
This is a random variable with the expectation value

$$H^{exp}=H^{obs}=\mathrm{log}(N)\frac{1}{N}\underset{i}{}N_i\mathrm{log}\left(N_i\right).$$
(33)
Let us assume now that the subwords are Bernoulli distributed. Then we find (omitting further the index $`n`$ which characterizes the length of the subwords) for the mean value

$`H^{obs}`$ $`=`$ $`\mathrm{log}N{\displaystyle \underset{i}{}}{\displaystyle \underset{N_i}{}}{\displaystyle \frac{(N1)!}{(N_i1)!(NN_i)!}}p_i^{N_i}`$ (35)
$`(1p_i)^{NN_i}\mathrm{log}N_i.`$
We will show now that the coefficients in this expansion are closely related to the Renyi entropies of order $`q`$ which are defined by

$$H^{(q)}=\frac{1}{1q}\mathrm{log}\left(\underset{i}{}p_i^q\right).$$
(36)
The Shannon entropies correspond to the limit $`q1`$ and the case $`q=0`$ is related to the total number of different subwords $`s`$:

$$H=H^{(1)},s=\lambda ^{H^{(0)}}.$$
(37)
In the limit of very long strings, i.e. $`Ns`$, the expected entropies are given by the approximation

$$H^{obs}=H\frac{s}{2N\mathrm{log}\lambda }.$$
(38)
In the opposite case, where $`Ns`$, the subwords may appear only very few times.  Therefore we obtain the series

$`H^{obs}`$ $`=`$ $`\mathrm{log}NN\mathrm{log}2{\displaystyle \underset{i}{}}p_{i}^{}{}_{}{}^{2}`$ (40)
$`{\displaystyle \frac{N^2}{2!}}\left(\mathrm{log}32\mathrm{log}2\right){\displaystyle \underset{i}{}}p_{i}^{}{}_{}{}^{3}+\mathrm{}`$
Using the definition of the Renyi entropies follows

$`H^{obs}`$ $`=`$ $`\mathrm{log}NN\mathrm{log}2\lambda ^{H^{(2)}}`$ (42)
$`{\displaystyle \frac{N^2}{2!}}\left(\mathrm{log}22\mathrm{log}2\right)\lambda ^{2H^{(3)}}\mathrm{}`$
The higher order Renyi entropies can be found by fitting the coefficients of $`N^k`$ in this series to the observed entropies. Since the entropies decrease monotonously with their order follows that the second order entropies ($`q=2`$) which can be obtained from eq.(42) represent a lower boundary of the first order Shannon entropies ($`q=1`$).

In this way we have obtained now three procedures to derive the entropies of human writings with a given length $`L`$:

1. Calculation of the usual entropy ($`q=1`$) from the observed distribution. Estimation of some finite length corrections using eq. (38). This is the standard procedure applied e.g. in .
2. Guessing first the “true” distribution $`Z(i,n)`$ on the basis of the observed distribution $`Z^{obs}(i,n,L)`$ and calculating then the entropies (for any $`q`$) from it. This method is essentially new. It was tested here but it still needs further elaboration.
3. Calculation of the higher order entropies (especially $`q=2`$) from the deviations between $`\mathrm{log}N`$ and the observed entropies. This method which seems to be new too is of limited accuracy and yields only rough estimates for the second and higher order entropies.

Let us consider now our examples: Beethoven’s Sonata is rather short ($`L=4920`$). Therefore the first approach is restricted to $`n<7`$ for the 3–symbol alphabet and to $`n<3`$ for the 32–symbol alphabet. We have obtained here as a new result the entropy for $`n=8`$ using method (2). The result is

$`H_8/\mathrm{log}(3)`$ $`=6.39`$ 3 symbols (43)
$`H_8/\mathrm{log}(32)`$ $`=2.49`$ $`\text{32 symbols}.`$ (44)
For $`9n26`$ the second order entropies were estimated by means of method (3). The result can be approximated by the fit formula

$$H_n/\mathrm{log}32.04\left(\mathrm{log}n+1\right),(\lambda =3).$$
(45)
This fit formula reminds the scaling of type eq. (18). The logarithmic function in eq. (45) yields a better fit than the power law postulated in an earlier paper .

Due to the rather large length of the text Moby Dick we may apply the method (1) up to much larger $`n`$–values. Taking into account that not all of the combinatorial possibilities are admitted we went up to $`n=16`$ for the 3–letter alphabet and up to $`n=10`$ for the 32–letter alphabet. The results are given in table II. We have checked these values also with a consistency test with method (2). A more detailed calculation by means of the distribution function method (2) will be discussed in . For larger $`n`$–values the second order entropies ($`q=2`$) were estimated from the differences between $`\mathrm{log}N`$ and the observed entropies (Fig. 8) using eq. (42). The result of our estimate is given by the fit formulae

$`H_n/\mathrm{log}(3)`$ $`=4.8\sqrt{n}7.6`$ $`(\lambda =3)`$ (46)
$`H_n/\mathrm{log}(32)`$ $`=0.9\sqrt{n}+1.7`$ $`(\lambda =32).`$ (47)
These fit formulae remind a scaling law according to eq. (19) with the exponent 0.5. The low accuracy of the data does not allow for a quantitative fit of the exponent, in fact any value between 0.4 and 0.6 gives a reasonable fit. A scaling law of square root type was first found by Hilberg  who fitted Shannon’s original data. This result was reproduced also for a textbook on selforganization .

Let us consider now the growth of the number of different subwords. According to relation (20) the number of different subwords in Moby Dick grows with $`n`$ according to a stretched exponential law

$$N_n^{}4.110^4\mathrm{exp}\left(5.2\sqrt{n}\right).$$
(48)
For the sonata the number of different subwords increases not as fast as for the English text, we find here a quadratic growth law

$$N_n^{}7.6n^2.$$
(49)
## V Conclusions

The present paper reports on results concerning information carrying strings such as texts and pieces of music. The results show that block entropies and the mutual information are appropriate measures of the correlations and the degree of order in strings. In agreement with earlier work  we have confirmed the hypothesis that strings produced by information processing sources are neither periodic nor chaotic but somehow in between. In the present work we have studied several examples (two books and one piece of music) to illustrate this hypothesis.

We have shown that there is some empirical evidence for short range and middle range correlations. This is to be seen in the statistics of pairs of letters and subwords. Another relevant information is contained in the block entropies. We developed methods how to guess these unknown quantities from relatively short samples. For subword length $`n100`$ we did not find any indication for the existence of pair or word correlations. Our empirical studies of the pair correlation function and the spectrum derived from this quantity  show a strong sensitivity on the length of the texts. In spite of all these difficulties we may state in conclusion that several distinct differences from chaotic strings are observed in texts or pieces of music.

The entropy scaling shows that texts and pieces of music resemble the sequences created by nonlinear dynamic systems near critical points through symbolic dynamics. However, more empirical data on long texts and pieces of music and a more detailed study of the block entropies of dynamical systems are needed to elaborate this point further. We state that with the present methods available so far it is not possible to calculate high order entropies for $`n>30`$ since there are no homogeneous texts which are significantly longer than several million letters.

###### Acknowledgements.

The authors thank H. Herzel and A. Schmitt for fruitful discussions, C. Frömmel for the permission to compare in Fig. 1 with an unpublished result on Mozart’s music and the Project Gutenberg Etext, Illinois Benedictine College, Lisle for providing the ASCII–text of “Moby Dick” and “Grimm’s Fairy Tales”.


format: text

text:
# Rashba spin-orbit interaction and shot noise for spin-polarized and entangled electrons

Permanent address: \]Department of Physics and Informatics, University of São Paulo at São Carlos, 13560-970 São Carlos/SP, Brazil.

## Abstract

We study shot noise for spin-polarized currents and entangled electron pairs in a four-probe (beam splitter) geometry with a local Rashba spin-orbit (s-o) interaction in the incoming leads. Within the scattering formalism we find that shot noise exhibits Rashba-induced oscillations with continuous bunching and antibunching. We show that entangled states and triplet states can be identified via their Rashba phase in noise measurements. For two-channel leads we find an additional spin rotation due to s-o induced interband coupling which enhances spin control. We show that the s-o interaction determines the Fano factor which provides a direct way to measure the Rashba coupling constant via noise.

Spin-related effects in transport form the basis of the emerging field of semiconductor spintronics springer . Moreover, the *electric* control of intrinsic *magnetic* degrees of freedom offers an important mechanism to manipulate and probe spin transport. For instance, the spin-transistor proposal of Datta and Das datta  highlights the relevance of a gate-controlled Rashba spin-orbit interaction as a means of spin rotating electron states in one-dimensional channels.

In this work, we investigate the transport properties of spin-polarized spin-pol ; egues  and spin-entangled entangler ; BLS  electrons for a beam-splitter configuration liu ; oliver ; henny  with a local Rashba spin-orbit (s-o) interaction, which acts within a finite region (length $`L`$) of the incoming leads, see Fig. 1 (left). Due to such a local s-o term the spinors of the incoming electrons can be varied continuously, which affects then the orbital symmetry of the wave function (via the Pauli principle) and thus transport properties such as current and noise. Within the scattering approach buttiker  we calculate current noise for leads with one and with two channels, see Fig. 1 (right). We find that shot noise for spin-polarized and entangled electrons strongly oscillates as a function of the Rashba coupling and the length $`L`$. In particular, singlet (triplet) pairs exhibit *intermediate* degrees of (anti-)bunching behavior BLS ; taddei . We show that entangled pairs as well as triplet states can be identified in noise experiments via their Rashba phase. We find that the s-o interaction determines the Fano factor, implying that the Rashba coupling can be measured via noise. Finally, we find an additional spin phase due to s-o induced interband coupling for leads with two channels. Since this modulation can be varied via the lateral confinement of the lead, this effect provides a new mechanism for electrical spin control.

*System.* We consider an experimentally feasible beam-splitter geometry liu ; oliver ; henny  with two incoming and two outgoing leads, Fig. 1. We assume that the local Rashba coupling in this lead can be externally controlled via a proper gating structure nitta ; grundler . The injected electrons (entangled or not) undergo a local spin evolution within an extension $`L`$ in lead 1. Note that we do not consider a Rashba interaction in lead 2; we are interested only in *phase differences* between leads 1 and 2. However, it is straightforward to extend our analysis to include s-o interaction in both incoming leads. We consider first single-channel leads (one occupied band), and then move on to the more involved case of two channels. For an incoming plane wave with wave vector *k* along the *x* direction the Rashba term is simply $`H_R=\alpha k\sigma _y`$ rashba ; molenkamp , where $`\alpha `$ denotes the s-o coupling constant.

*Spin transfer operator.* The rotation of the spin state can be described by a unitary transfer operator. This unitary operator can then be incorporated straightforwardly into the usual scattering matrix (Landauer-Büttiker) formalism for coherent transport buttiker  in order to calculate the current and the current correlators. We find that the transfer operator through the Rashba region is $`U_R=\mathrm{exp}\left(i\theta _R\sigma _y/2\right)`$, i.e., a rotation about the $`y`$ axis by an angle $`\theta _R=2\alpha m^{}L/\mathrm{}^2=2k_RL`$, $`m^{}`$ denotes the electron effective mass. This Rashba rotation is well known datta  and occurs only for incoming electrons having a spin component perpendicular to $`y`$. In the basis of the eigenstates of $`\sigma _z`$, we can write

$$U_R=\left(\begin{array}{cc}\mathrm{cos}\theta _R/2& \mathrm{sin}\theta _R/2\\ \mathrm{sin}\theta _R/2& \mathrm{cos}\theta _R/2\end{array}\right).$$
(1)
For instance, an incoming plane wave in lead 1 with spin up, $`|\psi _{\mathrm{in}}=|k`$, emerges at the other side of the Rashba region in the rotated state dynamical

$$x=L|\psi _{\mathrm{out}}=\left(\begin{array}{c}\mathrm{cos}\theta _R/2\\ \mathrm{sin}\theta _R/2\end{array}\right).$$
(2)
In the above we assume a *unity* transmission probability $`T_R`$ for electrons through the Rashba region. Indeed, $`T_R1`$ since there is no additional band offset due to different materials in our incoming lead mismatch . This implies that the Rashba interaction does not directly introduce noise in the lead – it simply rotates the incoming spin state. However, indirectly it does affect the noise characteristic of the entire system since it effectively changes the beam-splitter scattering matrix, as we shall see next.

*Scattering approach.* After leaving the Rashba region of length $`L`$ within lead 1 an electron is left in a linear superposition of spin-up and spin-down states with a phase $`\theta _R`$, Eq. (2). Both components of this superposition freely and independently move into the beam-splitter region where we assume they are partially transmitted to lead 3 or 4 . We can now combine the scattering matrices for the transmission of electrons through the Rashba region $`U_R`$ and for the spin-independent scattering of electrons from the incoming lead $`\beta `$ to the outgoing lead $`\alpha `$ at the beam splitter $`𝐬_{\alpha \beta }`$. The combined scattering matrices $`𝐬_{31}^R=𝐬_{31}U_R=rU_R`$ and $`𝐬_{41}^R=𝐬_{41}U_R=tU_R`$ describe both the Rashba evolution within lead 1 and the subsequent transmission into leads 3 and 4. Here, $`r`$ and $`t`$ denote the reflection and transmission amplitudes at the beamsplitter, respectively. Electrons from lead 2 are similarly transmitted into lead 3 or 4, but with no Rashba rotation. Here we have $`𝐬_{23}=𝐬_{14}`$ and $`𝐬_{24}=𝐬_{13}`$. Hence we define the total scattering matrix

$$𝐬=\left(\begin{array}{cccc}0& 0& 𝐬_{13}^R& 𝐬_{14}^R\\ 0& 0& 𝐬_{23}& 𝐬_{24}\\ 𝐬_{31}^R& 𝐬_{32}& 0& 0\\ 𝐬_{41}^R& 𝐬_{42}& 0& 0\end{array}\right),$$
(3)
relevant for calculating transport properties.

*Shot noise.* Shot noise blanter  is a non-equilibrium current fluctuation arising from the discrete nature of the charge flow (at zero temperature). At a time $`t`$, the current fluctuation about its average in lead $`\alpha `$ is $`\delta \widehat{I}_\alpha (t)=`$ $`\widehat{I}_\alpha (t)\widehat{I}_\alpha `$. As usual, shot noise is defined as the Fourier transform of the symmetrized current-current autocorrelation function between leads $`\gamma `$ and $`\mu `$

$$S_{\gamma \mu }(\omega )=\frac{1}{2}\delta \widehat{I}_\gamma (t)\delta \widehat{I}_\mu (t^{})+\delta \widehat{I}_\mu (t^{})\delta \widehat{I}_\gamma (t)e^{i\omega t}𝑑t.$$
(4)
The current in lead $`\gamma `$ in the scattering approach buttiker  is

$`\widehat{I}_\gamma (t)`$ $`=`$ $`{\displaystyle \frac{e}{h}}{\displaystyle \underset{\alpha \beta }{}}{\displaystyle 𝑑\epsilon 𝑑\epsilon ^{}e^{i(\epsilon \epsilon ^{})t/\mathrm{}}𝐚_\alpha ^{}(\epsilon )𝐀_{\alpha \beta }(\gamma ;\epsilon ,\epsilon ^{})𝐚_\beta (\epsilon ^{})},`$ (5)
$`𝐀_{\alpha \beta }(\gamma ;\epsilon ,\epsilon ^{})=\delta _{\gamma \alpha }\delta _{\gamma \beta }\mathrm{𝟏}𝐬_{\gamma \alpha }^{}(\epsilon )𝐬_{\gamma \beta }(\epsilon ^{}),`$
where $`𝐚_\alpha ^{}=(a_\alpha ^{},a_\alpha ^{})`$, and $`a_{\alpha \sigma }^{}(\epsilon )`$ $`[a_{\alpha \sigma }(\epsilon )]`$ denotes the creation (annihilation) fermionic operator for an electron with energy $`\epsilon `$ in lead $`\alpha `$; $`\sigma `$ is the spin component along a proper quantization direction ($`i=x,y,z`$). The *spin-dependent* s matrix is defined in Eq. (3). Below we determine explicit formulas for spin-polarized and spin-entangled electrons.

*Spin-polarized electrons*. For Fermi liquid leads we obtain the well-known noise formula buttiker  – but with a spin-dependent s matrix – after performing the ensemble average $`\mathrm{}`$ in Eq. (4). For spin-polarized electrons and a small bias $`eV`$ applied between the incoming (1,2) and outgoing (3,4) leads we find to linear order in $`eV`$ and at zero temperature

$$S_{33}^p=2eIT(1T)f_p,f_p=p\mathrm{sin}^2\frac{\theta _R}{2},$$
(6)
where $`T\left|t\right|^2`$ is the beam-splitter transmission, $`p`$ the degree of spin polarization in leads 1 and 2, and $`I=(2e^2/h)V/(1+p)`$ the mean current in lead 3.

*Spin-entangled electrons.* For electron pairs, the average in Eq. (4) is a quantum mechanical expectation value between two-electron states. We consider the following injected states delay  in leads 1 and 2

$`\begin{array}{c}|S\\ |Te_i\end{array}\}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[a_1^{}(\epsilon _1)a_2^{}(\epsilon _2)a_1^{}(\epsilon _1)a_2^{}(\epsilon _2)\right]|0\text{,}`$ (9)
$`|Tu_i`$ $`=`$ $`a_{1\sigma }^{}(\epsilon _1)a_{2\sigma }^{}(\epsilon _2)|0,\text{ }\sigma =,\text{,}`$ (10)
where $`|0`$ denotes the ground state (filled or not) Fermi sea of the leads. The states $`|S`$ and $`|Te_i`$ are the *entangled* singlet and triplets, respectively, while $`|Tu_i`$ are *unentangled* triplets. Here we have in mind an entangler entangler  attached to leads 1 and 2, Fig. 1. We assume that these pairs have *discrete* energies above $`\epsilon _F`$ BLS .

At zero temperature and applied voltage and $`\omega =0`$, the Fermi sea is completely inert and the noise in the system is solely due to the injected pairs above the Fermi surface BLS . For the singlet and triplets in Eq. (10) and the s matrix in Eq. (3), we find formula

$$S_{33}^X=\frac{2e^2}{h\nu }T(1T)f_X,$$
(11)
where the factor $`f_X`$, $`X=S`$, $`Te_i`$, $`Tu_i`$, $`i=x,y,z`$, depends on the Rashba phase $`\theta _R`$. The density of states $`\nu `$ in (11) arises because of the discrete levels. For the spin singlet

$$f_S=1+\mathrm{cos}\theta _R\delta _{\epsilon _1\epsilon _2},$$
(12)
where $`\epsilon _1`$ are $`\epsilon _2`$ denote the discrete energies of the paired electrons. For the triplet states in $`y`$ and $`z`$ directions

$`f_{Te_y}`$ $`=`$ $`1\mathrm{cos}\theta _R\delta _{\epsilon _1\epsilon _2},`$ (13)
$`f_{Tu_z}`$ $`=`$ $`1\mathrm{cos}^2(\theta _R/2)\delta _{\epsilon _1\epsilon _2},`$ (14)
$`f_{Tu_y}`$ $`=`$ $`f_{Te_z}=1\delta _{\epsilon _1\epsilon _2}.`$ (15)
*Two channels and s-o interband coupling.* So far we have considered a strictly 1D lead with a local Rashba interaction. Now we consider the case in which lead 1 has *two* transverse Rashba channels $`|a`$ and $`|b`$ moroz . We assume a weak s-o interband coupling which splits the bands near the crossing point $`k_c`$, Fig. 1 (right). To lowest order this splitting is $`2\alpha d`$, $`da|d/dy|b`$.

After traversing the Rashba region, a spin up (down) electron impinging at the band crossing is left in the state

$`e^{in\pi /2}\left(\begin{array}{c}e^{i\theta _R/2}\mathrm{cos}\left(\theta _d/2\right)\pm e^{i\theta _R/2}\\ i\left(e^{i\theta _R/2}\mathrm{cos}\left(\theta _d/2\right)e^{i\theta _R/2}\right)\end{array}\right){\displaystyle \frac{|a}{2}}`$ (18)
$`+e^{in\pi /2}\left(\begin{array}{c}ie^{i\theta _R/2}\mathrm{sin}\left(\theta _d/2\right)\\ e^{i\theta _R/2}\mathrm{sin}\left(\theta _d/2\right)\end{array}\right){\displaystyle \frac{|b}{2}},`$ (21)
where $`\theta _d=\theta _Rd/k_c`$ formula . To obtain (21) we have expanded an incoming spin-up $`n=0`$ (down, $`n=1`$) state in channel *a* in terms of the s-o interband-coupled states near the energy crossing at $`k_c`$, see Fig. 1. We describe these interband-coupled states perturbatively in analogy to the standard nearly-free electron model. Incoming electrons are now injected into linear combinations of unperturbed states (of channels *a* and *b*) near $`k_c`$ which satisfy proper boundary conditions for the velocity operator molenkamp .

Equation (21) clearly shows that impinging electrons with energies near the band crossing undergo further spin rotation $`\theta _d`$. This extra modulation arises because of channel mixing due to Rashba interband coupling. For $`\theta _d=0`$, Eq. (21) yields the state (2) with a single rotation $`\theta _R`$. An estimate of $`\theta _d`$ is readily obtained for infinite transverse confinement: assuming an energy (at the crossing)$`\epsilon _F=\epsilon (k_c)=24ϵ_R=24\mathrm{}^2k_R^2/2m^{}\alpha d/\epsilon (k_c)1/6`$ and $`d/k_c0.5`$; hence $`\theta _d=\theta _R/2`$ ($`\theta _R=\pi `$ for $`L=69`$ nm and $`\alpha =3.45\times 10^{11}`$ eVm nitta ; grundler ; for this $`\alpha `$ the lateral width of the channel is $`w=60`$ nm). Therefore even a “weak” interband coupling yields a sizable additional rotation $`\theta _d`$. This spin rotation for electrons injected at the band crossing produces an additional modulation of the transport properties. In particular, we find formula  for the spin-resolved charge current in lead 1

$$I_,1\pm \mathrm{cos}(\theta _d/2)\mathrm{cos}\theta _R,$$
(22)
which clearly shows the additional modulation $`\theta _d`$ hausler .

*Generalized Fano factors.* To determine shot noise in the presence of s-o interband coupling we proceed as before with the following extensions. For electron pairs, for instance, we consider the states in Eq. (10); here, however, the electron pair component in lead 1 evolves according to (21). After a somewhat lengthy calculation formula , we find that the Fano factors for the noise are now functions of both the Rashba angle $`\theta _R`$ and the interband mixing angle $`\theta _d`$,

$`f_p`$ $`=`$ $`{\displaystyle \frac{p}{2}}\left(1\mathrm{cos}{\displaystyle \frac{\theta _d}{2}}\mathrm{cos}\theta _R+{\displaystyle \frac{1}{2}}\mathrm{sin}^2{\displaystyle \frac{\theta _d}{2}}\right),`$ (23)
$`f_S`$ $`=`$ $`1+\left(\mathrm{cos}{\displaystyle \frac{\theta _d}{2}}\mathrm{cos}\theta _R\right)\delta _{\epsilon _1\epsilon _2},`$ (24)
$`f_{Te_z}`$ $`=`$ $`1{\displaystyle \frac{1}{2}}\left(\mathrm{cos}^2{\displaystyle \frac{\theta _d}{2}}+1\right)\delta _{\epsilon _1\epsilon _2},`$ (25)
$`f_{Tu_z}`$ $`=`$ $`1{\displaystyle \frac{1}{2}}(1+\mathrm{cos}{\displaystyle \frac{\theta _d}{2}}\mathrm{cos}\theta _R{\displaystyle \frac{1}{2}}\mathrm{sin}^2{\displaystyle \frac{\theta _d}{2}})\delta _{\epsilon _1\epsilon _2}.`$ (26)
The above equations reduce to the 1D case \[Eqs. (6) and (12)–(15)\] for $`\theta _d=0`$.

*Discussion.*  Figure 2 displays the “normalized” Fano factor $`fF/T(1T)`$, $`F=S/2eI`$, as a function of the Rashba phase $`\theta _R=2k_RL`$ for the (a) spin-polarized case \[Eq. (6) with $`p=1`$, here $`I=e^2V/h`$\], (b) injected singlet and triplet pairs \[Eqs. (12)–(15), here $`I=e/h\nu `$\]. In Fig. 2 we plot $`f`$ for two quantization directions: $`y`$ and $`z`$ ($`x`$ is equivalent to $`z`$), where “$`y`$” defines the Rashba rotation axis. In Fig. 2(a) only $`z`$-polarized electrons in leads 1 and 2 generate noise as $`\theta _R`$ is varied; $`y`$-polarized electrons are not affected by the Rashba rotation about $`y`$.

Because of the distinct symmetry of the orbital part of the pair wave function, shot noise for singlet and triplet states is *not* the same. As detailed in Ref. BLS  singlet pairs have a symmetric orbital wave function thus showing “bunching” behavior; triplets, on the other hand, show “antibunching” since their orbital wave function is antisymmetric. The Rashba phase modifies the symmetry of the spin part of the pair wavefunction; hence intermediate degrees of bunching or antibunching can be induced.

Figure 2(b) shows that shot noise for *entangled* singlet and triplet states display oscillatory bunching/antibunching behavior as a function of the Rashba phase (singlet and triplets differ by $`\pi `$). Via these oscillations it is possible to distinguish the entangled triplet states $`|Te_y`$ or $`|Te_z`$ from the respective unentangled ones $`|Tu_y`$ or $`|Tu_z`$. For $`\theta _R=0`$ the difference in noise vanishes, i.e., $`f_{Te_i}=f_{Tu_i}`$BLS .

Moreover, the oscillations in Fig. 2 suggest a direct way to obtain the s-o coupling constant $`\alpha `$ via measuring shot noise; for instance, from Eq. (6) ($`p=1`$) we find

$$\alpha =\frac{\mathrm{}^2}{m^{}L}\mathrm{arcsin}\sqrt{f_p}.$$
(27)
Figure 3 illustrates the effect of the additional spin-rotation $`\theta _d`$ (interband coupling) on the normalized Fano factor *f*; only the spin-polarized and the singlet cases are shown. This extra rotation can lead to a complete reversal of bunching/antibunching behavior for electrons near the band crossing \[see Fig. 1 (right)\]. Hence additional spin control is attained by varying the s-o interband coupling through the lateral width of the confining potential (e.g., via side gates).

*Conclusion.* Rashba s-o (interband) coupling strongly modulates current and shot noise for spin-polarized and entangled electrons in a beam-splitter geometry. This provides a means of probing spin properties in charge transport and offers a direct way to measure s-o coupling constants.

This work was supported by NCCR Nanoscience, the Swiss NSF, DARPA, and ARO. We acknowledge useful discussions with C. Schroll, H. Gassmann, and D. Saraga.


text:
# 1 Introduction

## 1 Introduction

In most situations in QCD, the question of whether a given observable can be calculated perturbatively reduces to two issues: firstly whether all the scales in the problem are large compared to $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ and secondly whether the observable is infrared and collinear safe.

One notable exception is the limit of high-energy scattering of two hard probes, for example quarkonia, or virtual photons. Fixed order analysis suggests that the scale of $`\alpha _\mathrm{s}`$ is simply determined by the inverse transverse size of the objects being scattered, independent of the centre of mass energy. But at very high centre-of-mass energies, $`\sqrt{s}`$, it becomes necessary to consider diagrams at all orders in the coupling. This resummation was initially considered for fixed coupling , but subsequently interest developed in its extension to the (more physical) running coupling case . One of the conclusions of these studies was that regardless of the hardness of the probes, there always exists a centre of mass energy beyond which the problem becomes entirely non-perturbative, see for example .

Other next-to-leading corrections , which however must be resummed, as is done for example in  or in other related  or less related  approaches, do not qualitatively change this picture. The transition to the non-perturbative regime takes place also in that case, although the diffusion effects can be different, due to the lower value of the intercept and the diffusion coefficient. This problem is now under detailed investigation .

For a long time, the view of the transition to the non-perturbative regime was based on the idea of diffusion. At lowest order in $`\alpha _\mathrm{s}`$ the dominant diagram at high energies is the exchange of a gluon between the two probes. The scale for $`\alpha _\mathrm{s}`$ is determined by the transverse momentum of the exchanged gluon, $`k^2`$, typically of order of transverse scale of the probes, $`Q^2`$ (in what follows we will generally refer to $`t=\mathrm{ln}Q^2/\mathrm{\Lambda }^2`$). At higher orders one needs to consider emissions of final-state gluons from the exchanged gluon — these emissions modify the transverse scale of the exchanged gluon, which undergoes a random walk in $`t^{}=\mathrm{ln}k^2/\mathrm{\Lambda }^2`$. In  this diffusion was illustrated graphically by plotting the median trajectory in $`t^{}`$ as a function of rapidity (the distance along the chain) as well as the dispersion around that median trajectory. A schematic version of such a plot is shown in figure 1, with the dark (blue) band representing the extent of the diffusion. Its characteristic shape led to it being named a *cigar* .

While at its ends the cigar is restricted to being narrow (by the presence of the probes), in the middle its width grows as a function of the centre of mass energy: $`\mathrm{\Delta }t^{}\sqrt{\overline{\alpha }_\mathrm{s}Y}`$ with $`Y=\mathrm{ln}s/Q^2`$. For the problem to remain perturbative, typical trajectories must remain well outside the non-perturbative infrared region and this leads to the requirement that the width of the cigar be less than $`t`$, i.e. $`\sqrt{\overline{\alpha }_\mathrm{s}Y}t`$. Using $`\overline{\alpha }_\mathrm{s}(t)=1/bt`$, with $`b=11/12`$ the coefficient of the $`1`$-loop $`\beta `$ function, this reduces to the condition $`Ybt^3`$. This limit was originally derived (with some additional numerical factors) in .

This picture of diffusion has its roots in the fixed-coupling approximation. It is tempting to consider the full running-coupling case as a small variation of the fixed-coupling situation. One then expects the differences to be that the centre of the cigar dips to $`t^{}`$ values below $`t`$, and that the cigar profile in $`t^{}`$ develops some asymmetry with respect to its centre (the cigar becomes a ‘banana’, figure 2a). This approach to the running-coupling case as a perturbation of the fixed-coupling case leads also to well-defined analytical predictions . In particular, the usual fixed-coupling exponential growth of the cross section, $`\mathrm{log}\sigma \alpha _\mathrm{s}Y`$, is supplemented by an extra term proportional to $`b^2\alpha _\mathrm{s}^5Y^3`$. This can be interpreted as due to the centre of the cigar being at $`t^{}`$ values of $`t𝒪\left(b\alpha _\mathrm{s}^3Y^2\right)`$, which for $`Y`$ of order $`bt^2`$ implies that the cigar reaches the non-perturbative domain. Indeed calculations of the ‘purely perturbative’ component of the cross section show unphysical (oscillating) behaviour for $`Ybt^2`$, which can be interpreted as meaning that the perturbative component is not a physically well-defined quantity in that region .

Recently however, it was discovered in the context of a simplified, collinear model for BFKL dynamics , that the running of the coupling can introduce further qualitative differences. Rather than the width of the cigar gradually increasing with $`Y`$ until it reaches the non-perturbative region at $`Ybt^3`$, or the centre of the cigar gradually moving in $`t^{}`$ down until it reaches the non-perturbative region at $`Ybt^2`$, it turns out that for $`Yt`$ an *abrupt* transition occurs and the dominant trajectory switches from being centred around $`t^{}=t`$ (fig. 1) to being centred around $`t^{}=\overline{t}`$, where $`\overline{t}0`$ is in the non-perturbative region (fig. 2b). For such a transition, which we named ‘tunneling’, there is no passage through any intermediate ‘banana’ configuration (fig. 2a).

In this article we wish to explain this phenomenon in more detail than was done previously and for the case of the full BFKL equation with running coupling. In section 2 we present analytical calculations which give the estimate of the nonperturbative Pomeron, while in section 3 we illustrate the phenomenon with a number of figures resulting from the numerical solution of the BFKL equation.

## 2 Analytical considerations

It is useful to recall the form of the BFKL equation as we shall use it here. We define a BFKL Green’s function $`G(Y,t,t_0)`$ through the following integro-differential equation

$$\frac{dG(Y,t,t_0)}{dY}=\frac{d^2\stackrel{}{q}}{\pi q^2}\overline{\alpha }_\mathrm{s}(\mathrm{ln}q^2)\left[\frac{k}{k^{}}G(Y,\mathrm{ln}k_{}^{}{}_{}{}^{2},t_0)G(Y,t,t_0)\mathrm{\Theta }(kq)\right],\begin{array}{c}k^2=e^t,\hfill \\ k^{}=|\stackrel{}{k}+\stackrel{}{q}|,\hfill \end{array}$$
(1)
with the initial condition

$$G(0,t,t_0)=\delta (tt_0),$$
(2)
and with $`\overline{\alpha }_\mathrm{s}=\alpha _\mathrm{s}N_C/\pi `$. All scales are in units of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. We shall take $`\overline{\alpha }_\mathrm{s}`$ as having its one-loop perturbative form, regularised by a cutoff at a non-perturbative scale $`\overline{t}`$:

$$\overline{\alpha }_\mathrm{s}(t)=\frac{\mathrm{\Theta }(t\overline{t})}{bt},b=\frac{11}{12}.$$
(3)
In some contexts we will use an equation similar to (1), but with $`\alpha _\mathrm{s}(\mathrm{ln}q^2)`$ replaced by $`\alpha _\mathrm{s}(t)`$.

It will be convenient to define some high-energy exponents: $`\omega _\mathrm{s}`$, the perturbative saddle-point exponent, determines the high-energy behaviour of the Green’s function when it is dominated by scales of order $`t`$. It is given by $`\omega _\mathrm{s}=(4\mathrm{ln}2)\overline{\alpha }_\mathrm{s}(t)+𝒪\left(\overline{\alpha }_\mathrm{s}^2\right)`$. We shall also need the non-perturbative pomeron exponent, $`\omega _{\mathrm{I}\mathrm{P}}`$, which governs the asymptotic high-energy behaviour of the Green’s function, $`G(Y)\mathrm{exp}(\omega _{\mathrm{I}\mathrm{P}}Y)`$, once the problem has become dominated by non-perturbative scales. It is given by $`\omega _{\mathrm{I}\mathrm{P}}=(4\mathrm{ln}2)\overline{\alpha }_\mathrm{s}(\overline{t})+𝒪\left(\overline{\alpha }_\mathrm{s}(\overline{t})^{5/3}\right)`$, though it is to be kept in mind that the expansion in powers of $`\overline{\alpha }_\mathrm{s}^{5/3}`$ converges very slowly.

For the purpose of our arguments we shall take the formal perturbative limit of $`t\overline{t}`$ and correspondingly $`\omega _\mathrm{s}\omega _{\mathrm{I}\mathrm{P}}`$.

### 2.1 Informal arguments for tunneling

The mechanism of tunneling is illustrated in figure 2b. To establish when tunneling takes place we have to equate the contribution from the tunneling configuration with that from the traditional perturbative configuration. Neglecting all prefactors, the tunneling configuration gives a contribution

$$G_{\mathrm{tunnel}}e^{(t\overline{t})/2}e^{\omega _{\mathrm{I}\mathrm{P}}Y}e^{(t\overline{t})/2},$$
(4)
corresponding (reading from left to right) to the ‘penalty’ for branching from $`t`$ to $`\overline{t}`$, the BFKL evolution at $`\overline{t}`$ and the penalty for branching back to $`t`$. In contrast, the perturbative contribution is proportional to

$$G_{\mathrm{pert}}e^{\omega _\mathrm{s}Y}.$$
(5)
Equating these two expressions gives

$$Y_{\mathrm{tunnel}}(t)=\frac{t\overline{t}}{\omega _{\mathrm{I}\mathrm{P}}\omega _\mathrm{s}(t)}.$$
(6)
Perhaps the most important characteristic of this equation is the fact that $`Y_{\mathrm{tunnel}}`$ is proportional to $`t`$ rather than $`t^3`$ as for the diffusion mechanism. One possible exception is the case in which the Pomeron is weak, that is $`\omega _{\mathrm{I}\mathrm{P}}1`$. Then $`\overline{t}`$ should be quite large, and including in Eq.(6) the value $`\omega _\mathrm{s}(t)=\chi _m/(bt)`$ and the rough estimate $`\omega _{\mathrm{I}\mathrm{P}}=\chi _m/(b\overline{t})`$ yields $`Y_{\mathrm{tunnel}}=(bt/\chi _m)\overline{t}`$, so that the perturbative regime is kept up to $`\omega _\mathrm{s}(t)Y`$ values of order $`\overline{t}`$.

### 2.2 Formal arguments

For an isolated Pomeron “bound state”, the tunneling probability is estimated on the basis of the approximate decomposition

$$G(Y;t,t_0)=G_{\mathrm{pert}}(Y;t,t_0)+e^{\omega _{\mathrm{I}\mathrm{P}}Y}\frac{_{\omega _{\mathrm{I}\mathrm{P}}}(t)_{\omega _{\mathrm{I}\mathrm{P}}}(t_0)}{𝑑t_{\omega _{\mathrm{I}\mathrm{P}}}(t)^2},$$
(7)
where $`_{\omega _{\mathrm{I}\mathrm{P}}}(t)`$ is the regular solution for $`t\mathrm{}`$ of the BFKL equation , satisfying the strong-coupling boundary condition also. The decomposition (7) separates the continuum from the lowest bound state contributions to the gluon’s Green functions  or, equivalently, the background integral from the leading Regge pole . The normalisation of the bound state factor in front of $`e^{\omega _{\mathrm{I}\mathrm{P}}Y}`$ is simply obtained by imposing that it be a projector by $`t_0`$ integration, or equivalently from the requirement that

$$G(Y;t,t_0)=𝑑t^{}G(Y;t,t^{})G(Yy;t^{},t_0).$$
(8)
Eq.(7) applies as it stands to a symmetrical BFKL kernel, as that with the scale choice $`\overline{\alpha }_\mathrm{s}(\mathrm{ln}q^2)`$ in Eq.(1). If instead the scale $`\overline{\alpha }_\mathrm{s}(t)`$ is chosen, the kernel is asymmetrical, the conjugated wave function is $`\overline{}_{\omega _{\mathrm{I}\mathrm{P}}}(t)=t_{\omega _{\mathrm{I}\mathrm{P}}}^{}(t)`$, and the bound state projectors differ by a $`t_0`$ factor and by the normalisation integral

$$G(Y;t,t_0)=G_{\mathrm{pert}}(Y;t,t_0)+e^{\omega _{\mathrm{I}\mathrm{P}}Y}\frac{_{\omega _{\mathrm{I}\mathrm{P}}}(t)_{\omega _{\mathrm{I}\mathrm{P}}}(t_0)t_0}{𝑑tt_{\omega _{\mathrm{I}\mathrm{P}}}(t)^2}.$$
(9)
In either case, the large $`t`$ behaviour of $`_{\omega _{\mathrm{I}\mathrm{P}}}(t)`$ is estimated in a semiclassical approach in which $`b\overline{\alpha }_\mathrm{s}(t),1/Y`$ and also $`b\omega _{\mathrm{I}\mathrm{P}}`$ are taken to be small parameters (weak Pomeron, or $`b`$ \- expansion ). In such a case the $`\gamma `$ \- representation of $`_{\omega _{\mathrm{I}\mathrm{P}}}`$ is valid, and takes the form

$$\begin{array}{c}_{\omega _{\mathrm{I}\mathrm{P}}}(t)=_{1/2i\mathrm{}}^{1/2+i\mathrm{}}\frac{d\gamma }{2\pi i}e^{(\gamma 1/2)t\frac{1}{b\omega _{\mathrm{I}\mathrm{P}}}X_{\omega _{\mathrm{I}\mathrm{P}}}(\gamma )}\hfill \\ \hfill \sqrt{\frac{b\omega _{\mathrm{I}\mathrm{P}}}{2\pi \chi ^{}[p(b\omega _{\mathrm{I}\mathrm{P}}t)]}}\mathrm{exp}\left[\frac{1}{b\omega _{\mathrm{I}\mathrm{P}}}_{\chi _m}^{b\omega _{\mathrm{I}\mathrm{P}}t}p(x)𝑑x\right],\end{array}$$
(10)
where the saddle point condition at $`\overline{\gamma }=\frac{1}{2}p`$

$$\chi \left(\frac{1}{2}p\right)=\frac{X(\gamma )}{\gamma }|_{\gamma =\frac{1}{2}p}=b\omega _{\mathrm{I}\mathrm{P}}t,$$
(11)
defines the semiclassical “momentum” $`p(b\omega _{\mathrm{I}\mathrm{P}}t)`$. The effective eigenvalue function in Eq.(11) is just the leading characteristic function $`\chi _0(\gamma )`$ for the $`\overline{\alpha }_\mathrm{s}(t)`$ choice, and differs by some corrections involving the parameter $`\omega _{\mathrm{I}\mathrm{P}}`$ in the symmetric scale choice $`\overline{\alpha }_\mathrm{s}(\mathrm{ln}q^2)`$ which essentially lead to a rescaling of the (right) regular solution by a constant factor $`1/\sqrt{\chi _m}`$ (see Appendix) .

The Pomeron suppression arising from Eq.(10) has a universal form $`e^{\left(\frac{1}{b\overline{\alpha }_\mathrm{s}(t)}g(\overline{\alpha }_\mathrm{s}(t))\right)}`$, where the function $`g(\overline{\alpha }_\mathrm{s})`$ ($`g(0)0`$) is easily found from Eq.(11), . On the other hand the perturbative part, corresponding to the continuum, is dominated by the saddle point exponent $`\omega _\mathrm{s}(t)`$

$$G_{\mathrm{pert}}(t,t_0;Y)=\frac{1}{\sqrt{2\pi \chi _m^{^{\prime \prime }}\overline{\alpha }_\mathrm{s}(\frac{t+t_0}{2})Y}}\mathrm{exp}\left(\omega _\mathrm{s}\left(\frac{t+t_0}{2}\right)Y+\mathrm{corrections}\right),$$
(12)
where the diffusion corrections are due to the running coupling .

From Eqs.(7,9,10) we estimate the $`t`$-dependence of the Pomeron contributions for the choices $`\overline{\alpha }_\mathrm{s}(k^2)`$ ($`\overline{\alpha }_\mathrm{s}(q^2)`$) of the running coupling scale respectively

$`G_{\mathrm{I}\mathrm{P}}^{(k)}`$ $`=`$ $`{\displaystyle \frac{N^{(k)}(\overline{t})}{2\pi b\omega _{\mathrm{I}\mathrm{P}}t}}\mathrm{exp}\left[t+{\displaystyle \frac{2}{\omega _{\mathrm{I}\mathrm{P}}}}\left(\mathrm{log}b\omega _{\mathrm{I}\mathrm{P}}t+1+\psi (1)+𝒪\left({\displaystyle \frac{1}{(b\omega _{\mathrm{I}\mathrm{P}}t)^3}}\right)\right)\right],`$
$`G_{\mathrm{I}\mathrm{P}}^{(q)}`$ $`=`$ $`{\displaystyle \frac{N^{(q)}(\overline{t})}{2\pi (b\omega _{\mathrm{I}\mathrm{P}}t)^2}}\left({\displaystyle \frac{b\omega _{\mathrm{I}\mathrm{P}}}{\chi _m}}\right)\mathrm{exp}\left[t+{\displaystyle \frac{2}{\omega _{\mathrm{I}\mathrm{P}}}}\left(\mathrm{log}b\omega _{\mathrm{I}\mathrm{P}}t+1+\psi (1)+𝒪\left({\displaystyle \frac{1}{(b\omega _{\mathrm{I}\mathrm{P}}t)^3}}\right)\right)\right].`$
Notice that they differ both in the $`t`$-dependence, because of the $`t`$ factor on Eq.(9) and for the $`\overline{t}`$-dependent normalisation integrals

$`N^{(k)}(\overline{t})=\left({\displaystyle _{\overline{t}}^{\mathrm{}}}𝑑tt_{\omega _{\mathrm{I}\mathrm{P}}}^2(t)\right)^1,`$
$`N^{(q)}(\overline{t})=\left({\displaystyle 𝑑t_{\omega _{\mathrm{I}\mathrm{P}}}^2(t)}\right)^1,`$ (14)
The latter quantities cannot be estimated in the semiclassical approximation, which breaks down at the turning points, and they need to be evaluated numerically. We can roughly estimate their size on the basis of the Airy diffusion model  obtained by the quadratic expansion of the characteristic function

$$\chi _\omega (\gamma )=\chi _m+\frac{1}{2}\chi _m^{\prime \prime }(\gamma \gamma _m)^2+\mathrm{}.$$
(15)
This approximation is expected to be a reliable one for the Pomeron dynamics provided the fluctuations $`(\gamma \gamma _m)^2\frac{2b\omega _{\mathrm{I}\mathrm{P}}}{\chi _m^{\prime \prime }}`$ are small compared to $`4!\chi _m^{\prime \prime }/(2\chi _m^{(4)})`$, and $`t`$ varies in a range for which

$$\xi _{\omega _{\mathrm{I}\mathrm{P}}}(t)=\left(\frac{2b\omega _{\mathrm{I}\mathrm{P}}}{\chi _m^{\prime \prime }}\right)^{\frac{1}{3}}\left(t\frac{\chi _m}{b\omega _{\mathrm{I}\mathrm{P}}}\right),$$
(16)
is a parameter of order unity. In this case the BFKL equation with the choice $`\overline{\alpha }_\mathrm{s}(k^2)`$ reduces, by Eq.(15) to a second-order differential equation in $`t`$ , and the cut-off condition for $`\overline{\alpha }_\mathrm{s}(t)`$ at $`t=\overline{t}`$ (3) implies that the regular solution for $`t\mathrm{}`$ must vanish at $`\overline{t}`$ , and for $`t`$ around $`\overline{t}`$ takes the following form

$$\left(\frac{\chi _m^{\prime \prime }}{2b\omega }\right)^{\frac{1}{3}}\stackrel{~}{}_\omega (t)=\pi \left[\mathrm{Bi}[\xi _\omega (t)]\frac{\mathrm{Bi}[\overline{\xi }_\omega (\overline{t})]}{\mathrm{Ai}[\overline{\xi }_\omega (\overline{t})]}\mathrm{Ai}[\xi _\omega (t)]\right],$$
(17)
while the regular solution for $`t\mathrm{}`$ becomes

$$\left(\frac{\chi _m^{\prime \prime }}{2b\omega }\right)^{\frac{1}{3}}_\omega (t)=\mathrm{Ai}\left[\xi _\omega (t)\right],$$
(18)
with relative normalisation fixed by the Wronskian

$$\mathrm{W}(_\omega ,\stackrel{~}{}_\omega )=\frac{2b\omega }{\chi _m^{\prime \prime }}.$$
(19)
Using now the Airy form for the regular solution (18) we can evaluate the normalisation integrals (14) to get

$`N^{(k)}(\overline{t})`$ $`=`$ $`{\displaystyle \frac{\chi _m}{b\omega _{\mathrm{I}\mathrm{P}}}}\left({\displaystyle \frac{2b\omega _{\mathrm{I}\mathrm{P}}}{\chi _m^{^{\prime \prime }}}}\right)^{\frac{1}{3}}{\displaystyle _{\xi _0}^{\mathrm{}}}𝑑\xi [\mathrm{Ai}(\xi )]^2+{\displaystyle _{\xi _0}^{\mathrm{}}}𝑑\xi \xi [\mathrm{Ai}(\xi )]^2=`$
$`=`$ $`\xi _0[\mathrm{Ai}^{^{}}(\xi _0)]^2{\displaystyle \frac{\frac{2}{3}\omega _\mathrm{s}(\overline{t})+\frac{1}{3}\omega _{\mathrm{I}\mathrm{P}}}{\omega _\mathrm{s}(\overline{t})\omega _{\mathrm{I}\mathrm{P}}}},`$
$`\chi _mN^{(q)}(\overline{t})`$ $`=`$ $`\left({\displaystyle \frac{2b\omega _{\mathrm{I}\mathrm{P}}}{\chi _m^{^{\prime \prime }}}}\right)^{\frac{1}{3}}{\displaystyle _{\xi _0}^{\mathrm{}}}𝑑\xi [\mathrm{Ai}(\xi )]^2=\xi _0{\displaystyle \frac{\omega _\mathrm{s}(\overline{t})}{\omega _\mathrm{s}(\overline{t})\omega _{\mathrm{I}\mathrm{P}}}}{\displaystyle \frac{b\omega _{\mathrm{I}\mathrm{P}}}{\chi _m}}[\mathrm{Ai}^{^{}}(\xi _0)]^2,`$ (20)
where we have used the Airy identities  and the assumption that the Pomeron intercept is given by the rightmost zero of the Airy $`\mathrm{Ai}(\xi )`$ function

$$\left(\frac{2b\omega _{\mathrm{I}\mathrm{P}}}{\chi _m^{\prime \prime }}\right)^{\frac{1}{3}}\left(\overline{t}\frac{\chi _m}{b\omega _{\mathrm{I}\mathrm{P}}}\right)=\xi _02.3381,\mathrm{Ai}(\xi _0)=0.$$
(21)
The additional factor $`\chi _m`$ which multiplies $`N^{(q)}(\overline{t})`$ comes from the rescaling of the (right)regular eigenfunctions due to a scale change $`_{\omega _{\mathrm{I}\mathrm{P}}}^{(q)}(t)\sqrt{\chi _m}_{\omega _{\mathrm{I}\mathrm{P}}}^{(k)}(t)`$.

The expressions (2.2) supplemented by the normalisation factors (2.2) can be compared with the naive estimate (6) by rewriting them in the form

$`G_{\mathrm{I}\mathrm{P}}^{(k)}`$ $`=`$ $`\overline{\alpha }_\mathrm{s}(t)\mathrm{exp}\left[(t\overline{t})+{\displaystyle \frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}}\mathrm{log}{\displaystyle \frac{t}{\overline{t}}}+\mathrm{\Delta }^{(k)}\right],`$
$`G_{\mathrm{I}\mathrm{P}}^{(q)}`$ $`=`$ $`\overline{\alpha }_\mathrm{s}^2(t)\mathrm{exp}\left[(t\overline{t})+{\displaystyle \frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}}\mathrm{log}{\displaystyle \frac{t}{\overline{t}}}+\mathrm{\Delta }^{(q)}\right],`$ (22)
where one sees the explicit $`\overline{\alpha }_\mathrm{s}(t)`$, ($`\overline{\alpha }_\mathrm{s}(t)^2`$) dependence for $`\overline{\alpha }_\mathrm{s}(k^2)`$ and $`\overline{\alpha }_\mathrm{s}(q^2)`$ scale choices respectively. The $`\mathrm{log}\frac{t}{\overline{t}}`$ in the exponent represents a collinear enhancement due to the evolution from scale $`t`$ to scale $`\overline{t}`$, and the $`\overline{t}`$ \- dependent factors

$`\mathrm{exp}(\mathrm{\Delta }^{(k)})`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \xi _0[\mathrm{Ai}^{}(\xi _0)]^2}}{\displaystyle \frac{[\omega _\mathrm{s}(\overline{t})\omega _{\mathrm{I}\mathrm{P}}]}{\frac{1}{3}\omega _{\mathrm{I}\mathrm{P}}+\frac{2}{3}\omega _\mathrm{s}(\overline{t})}}{\displaystyle \frac{1}{\omega _{\mathrm{I}\mathrm{P}}}}\mathrm{exp}\left[{\displaystyle \frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}}(1+\psi (1)+\mathrm{log}b\omega _{\mathrm{I}\mathrm{P}}\overline{t})\overline{t}\right],`$
$`\mathrm{exp}(\mathrm{\Delta }^{(q)})`$ $`=`$ $`{\displaystyle \frac{\chi _m}{2\pi \xi _0[\mathrm{Ai}^{}(\xi _0)]^2}}{\displaystyle \frac{[\omega _\mathrm{s}(\overline{t})\omega _{\mathrm{I}\mathrm{P}}]}{\omega _\mathrm{s}(\overline{t})\omega _{\mathrm{I}\mathrm{P}}^2}}\mathrm{exp}\left[{\displaystyle \frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}}(1+\psi (1)+\mathrm{log}b\omega _{\mathrm{I}\mathrm{P}}\overline{t})\overline{t}\right],`$ (23)
come from the normalisation factors (2.2) and of the semiclassical prefactors. By then comparing the result (2.2) with the perturbative Green’s function in Eq.(12) we obtain the estimate of the tunneling transition rapidity in both scale choice cases

$`(\omega _{\mathrm{I}\mathrm{P}}\omega _\mathrm{s}(t))Y_{\mathrm{tunnel}}^{(k)}`$ $`=t\overline{t}{\displaystyle \frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}}\mathrm{log}{\displaystyle \frac{t}{\overline{t}}}\mathrm{log}\overline{\alpha }_\mathrm{s}(t){\displaystyle \frac{1}{2}}\mathrm{log}(2\pi \chi _m^{\prime \prime }\overline{\alpha }_\mathrm{s}(t)Y_{\mathrm{tunnel}})\mathrm{\Delta }^{(k)}(\overline{t}),`$
$`(\omega _{\mathrm{I}\mathrm{P}}\omega _\mathrm{s}(t))Y_{\mathrm{tunnel}}^{(q)}`$ $`=t\overline{t}{\displaystyle \frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}}\mathrm{log}{\displaystyle \frac{t}{\overline{t}}}2\mathrm{log}\overline{\alpha }_\mathrm{s}(t){\displaystyle \frac{1}{2}}\mathrm{log}(2\pi \chi _m^{\prime \prime }\overline{\alpha }_\mathrm{s}(t)Y_{\mathrm{tunnel}})\mathrm{\Delta }^{(q)}(\overline{t}),`$
### 2.3 Limits of the tunneling picture

We have just seen that the Pomeron regime sets in, according to Eqs.(6) and (LABEL:eq:YTunnelExact) at $`\omega _\mathrm{s}(t)Y`$ values of order $`\overline{t}`$. This is to be compared with the validity boundary of the hard Pomeron behaviour  $`1\omega _\mathrm{s}(t)Yt/\sqrt{D}`$, where $`D=\frac{\chi _m^{\prime \prime }}{2\chi _m}`$ is the diffusion coefficient. Thus, for phenomenological values of $`t`$ and $`\overline{t}`$ there appears to be little room for the hard Pomeron (and its diffusion corrections) to be observable, at least for models (such as the leading BFKL equation) in which the hard Pomeron and the nonperturbative Pomeron are quite strong. The Pomeron contamination in such cases, is a serious problem even for a theoretical determination of perturbative diffusion corrections to the hard Pomeron. In fact, one has to use the $`b`$-expansion $`b0`$ limit with $`\overline{\alpha }_\mathrm{s}(t)`$ and $`\overline{\alpha }_\mathrm{s}(\overline{t})`$ fixed — in order to identify such corrections at the numerical level.

The situation may change, if subleading and unitarity corrections are taken into account. Subleading contributions to the kernel are large  and need to be resummed : they provide eventually a sizeable decrease of the hard Pomeron intercept and even more of the diffusion coefficient. This effect slows down both the nonperturbative effects and diffusion corrections, thus increasing the range of validity of the perturbative predictions. Furthermore, unitarity effects are expected to be important for large densities and more so in the strong coupling region. All this goes in the direction of a weaker asymptotic Pomeron — as is the actual physical Pomeron of soft physics  — thus increasing the effective $`\overline{t}`$ value, and the window of validity of the hard Pomeron regime.<sup>2</sup><sup>2</sup>2Indeed if the soft Pomeron has weaker energy dependence than the hard pomeron, then tunneling stops being a favoured contribution altogether since the non-perturbative Pomeron cannot overtake the hard pomeron before the oscillatory regime is reached.

It may thus be that in realistic small-$`x`$ models  with a weaker Pomeron, diffusion corrections to the hard Pomeron may set in before tunneling, so that the transition to the Pomeron regime is smoother. In any case, before diffusion corrections lead to a decreasing cross-section and to the asymptotic oscillatory regime the Pomeron is expected to take over.

## 3 Pictorial results

### 3.1 Cigars

Perhaps the most dramatic illustration of the phenomenon of tunneling can be obtained through pictorial representations analogous to those of figures 1 and 2, but obtained by exact numerical solution of the BFKL equation.

This involves studying the running-coupling dynamics of BFKL evolution and in particular the sensitivity of the Green’s function $`G(Y,t,t)`$ to evolution at some intermediate point $`y,t^{}`$. This can be done by examining:

$$f(Y,y,t,t^{})=\frac{G(y,t,t^{})G(Yy,t^{},t)}{G(Y,t,t)},$$
(25)
which, by virtue of the linearity of the BFKL equation, is normalised as follows

$$𝑑t^{}f(Y,y,t,t^{})=1.$$
(26)
From the point of view of graphical representation this quantity suffers from the problem that for $`y0`$ or $`yY`$ and $`t^{}t`$ one is sensitive to the $`\delta `$-function initial condition for $`G`$, eq. (2). To smooth-out this $`\delta `$-function we examine instead a smeared version of $`f`$:

$$\overline{f}(Y,y,t,t^{})=\frac{𝑑t_0𝑑t_1e^{(t_0t)^2/2(t_1t)^2/2}G(y,t_0,t^{})G(Yy,t^{},t)}{𝑑t_0𝑑t_1e^{(t_0t)^2/2(t_1t)^2/2}G(Y,t_0,t_1)}$$
(27)
Contour plots of $`\overline{f}(Y,y,t,t^{})`$ as a function of $`y`$ and $`t^{}`$ are shown in figure 3 for various values of $`Y`$. We have chosen a value for $`t`$ that is physically unrealistic, $`t=30`$. The reason for this choice is that it enhances the separation between the perturbative and non-perturbative regions, allowing us to illustrate more clearly the transition that takes place.<sup>3</sup><sup>3</sup>3For more phenomenologically relevant plots we note that in addition to taking more reasonable values for $`t`$ it would also be necessary to include the NLL corrections and collinear resummation. The coupling is regularised in the infrared, eq. (3), by a cutoff at $`\overline{t}=4`$.

At the lowest value of $`Y`$ in figure 3 the shape of the cigar is largely determined by the smearing of the initial condition, eq. (27). At $`Y=50`$ one can see a significant amount of diffusion at central values of $`y`$, as well as a slight asymmetry favouring lower $`t^{}`$ values, caused by the running of the coupling.

At $`Y=60`$ the tunneling transition starts to take place: at low values of $`t^{}t`$, an enhanced region appears which is *disconnected* from the main region of the cigar — it has not appeared through diffusion, but through tunneling. As $`Y`$ increases the perturbative region of $`t^{}t`$ becomes progressively less important, eventually being completely supplanted by the non-perturbative region $`t^{}\overline{t}`$. In this transition stage we take the liberty of referring to the contour plots as ‘BFKL aliens’.

Beyond $`Y=75`$, the evolution takes place entirely in the non-perturbative region, except in the wings, $`y0`$ and $`yY`$, which connect the initial conditions at $`t`$ to the infrared region. We refer to these contour plots as ‘bowls’.

We note that the position in $`Y`$ of the tunneling transition is quite consistent with the arguments of section 2. For our particular non-perturbative conditions, the non-perturbative power is $`\omega _{\mathrm{I}\mathrm{P}}0.32`$; substituting this into eq. (LABEL:eq:YTunnelExact) gives $`Y_{\mathrm{tunnel}}68`$.

### 3.2 Ducks

A convenient way of visualising the tunneling transition within a single plot is to consider $`\overline{f}(Y,Y/2,t,t^{})`$ as a function of $`Y`$ and $`t^{}`$, i.e. a contour plot built up of cross sections through the centres of the plots of figure 3, for a range of $`Y`$ values. This is shown in the upper plot of figure 4. At low $`Y`$ we see just the width of the initial condition, which gradually diffuses out while remaining centred at $`t^{}t`$. At $`Y60`$ a second region appears ‘out of nowhere’ at $`t^{}\overline{t}`$ and the original cigar region disappears completely just beyond. This, once again, is the tunneling transition. We refer to the resulting figure as a ‘BFKL duck’.

In the lower plot of figure 4 we show two quantities as a function of $`Y`$ (on the same $`Y`$ scale as the upper plot): the effective exponent $`\omega _{\mathrm{eff}}`$, defined as

$$\omega _{\mathrm{eff}}=\frac{d}{dY}\mathrm{ln}𝑑t_0𝑑t_1e^{(t_0t)^2/2(t_1t)^2/2}G(Y,t_0,t_1),$$
(28)
and the position of the maximum in $`t^{}`$ of $`\overline{f}(Y,Y/2,t,t^{})`$, which we refer to as $`t_{\mathrm{max}}`$. At the point where the tunneling transition occurs there is a jump in the value of $`t_{\mathrm{max}}`$ as it switches from the perturbative to the non-perturbative region. Simultaneously there is a rapid change in $`\omega _{\mathrm{eff}}`$ as it goes from the perturbative exponent $`\omega _\mathrm{s}`$ to the non-perturbative one $`\omega _{\mathrm{I}\mathrm{P}}`$.

We can use the point where the position of the maximum goes from being above to below $`t^{}=\frac{1}{2}(t+\overline{t})`$ as the definition of where tunneling occurs, $`Y_{\mathrm{tunnel}}`$ (in what follows, since we will no longer be faced with the problem of visualising an initial $`\delta `$-function, we will actually consider $`f(Y,Y/2,t,t^{})`$ rather than $`\overline{f}`$, though the difference in practice is minimal). Alternatively one can define the tunneling transition as occurring when the integrals of $`f(Y,Y/2,t,t^{})`$ above and below $`\frac{1}{2}(t+\overline{t})`$ are equal,

$$_{\mathrm{}}^{(t+\overline{t})/2}𝑑t^{}f(Y,Y/2,t,t^{})=_{(t+\overline{t})/2}^{\mathrm{}}𝑑t^{}f(Y,Y/2,t,t^{}).$$
(29)
These two definitions for $`Y_{\mathrm{tunnel}}`$ have been used to deduce, from our numerical solutions, the points with errors bars in fig. 5 (the lower error bar corresponds to the first definition, while the upper error bar corresponds to the second one), where we plot $`Y_{\mathrm{tunnel}}`$ as a function of $`t`$, for two values of $`\overline{t}`$. They are to be compared to the smooth lines, which correspond to the analytical predictions both in $`\overline{\alpha }_\mathrm{s}(t)`$ and $`\overline{\alpha }_\mathrm{s}(\mathrm{ln}q^2)`$ cases. We have tried the direct evaluation of the $`G_{\mathrm{I}\mathrm{P}}`$ from Eqs.(7), (9) by numerically calculating the $`\gamma `$-representation integrals (solid lines) and the analytical approximation based on the Airy model Eq. (LABEL:eq:YTunnelExact) (dashed lines).

One sees a remarkable agreement between the ‘numerically measured’ and analytical results, which would seem to be strong confirmation of even the fine details of our picture. We note in particular the roughly linear dependence of $`Y_{\mathrm{tunnel}}`$ on $`t`$ and the fact that its slope increases for larger $`\overline{t}`$ as a consequence of the smaller $`\omega _{\mathrm{I}\mathrm{P}}`$ and hence smaller denominator in eq. (6).

## 4 Conclusions

In this paper we have presented a detailed investigation of the abrupt tunneling transition from the perturbative regime to the non-perturbative ‘Pomeron’ regime for the case of leading-log BFKL evolution with a running coupling.

The tunneling transition takes place when $`\omega _\mathrm{s}Y\overline{t}`$, with $`\mathrm{exp}(\overline{t})`$ the effective non-perturbative scale of the problem (in units of $`\mathrm{\Lambda }_{QCD}^2`$). This is parametrically much earlier than estimates of the maximum perturbatively accessible $`Y`$ value that are based on diffusion arguments. Furthermore, in contrast to expectations from diffusion arguments, the transition does not involve evolution at scales intermediate between those of the hard probe and the non-perturbative region, as shown in figs. 3 and 4.

Since the purpose of this paper has been to give an illustration of the effect of tunneling, we have restricted ourselves to leading-logarithmic evolution and considered only very extreme scales for the hard-probe, where the effect is more dramatic because of the large ratio $`\omega _{\mathrm{I}\mathrm{P}}/\omega _\mathrm{s}(t)`$, or equivalently $`t/\overline{t}`$.

In a more phenomenologically relevant context there remain a number of other important issues: the separation between tunneling and diffusion may be less distinct, both because the accessible range of $`t`$ is smaller, and because for a given ratio of $`t/\overline{t}`$ higher-order corrections tend to reduce $`\omega _{\mathrm{I}\mathrm{P}}/\omega _\mathrm{s}(t)`$. Furthermore unitarity and saturation corrections may actually cause the effective value of $`\omega _{\mathrm{I}\mathrm{P}}`$ to be smaller than $`\omega _\mathrm{s}`$ which could eliminate tunneling altogether. Developing a better understanding of these and related questions remains one of the main challenges of small-$`x`$ physics.

## Acknowledgements

One of us (GPS) wishes to thank the Dipartimento di Fisica, Università di Firenze, for hospitality while part of this work was carried out.

## Appendix A: Pomeron weight with $`\overline{𝜶}_𝐬\mathbf{(}\mathrm{𝐥𝐧}𝒒^\mathrm{𝟐}\mathbf{)}`$

Taking $`\overline{\alpha }_\mathrm{s}(\mathrm{ln}q^2)`$ instead of $`\overline{\alpha }_\mathrm{s}(t)`$, as in Eq.(1) means adding a term $`\frac{b}{2}(\chi _0^{}+\chi _0^2)`$ to the NL kernel which emerges due to the evolution from scale $`t`$ to scale $`q^2`$. Correspondingly, in the $`\omega `$ \- expansion formulation  the effective characteristic function becomes

$$\chi _\omega =\chi _0(\gamma )+\omega \frac{1}{\chi _0}\left[\frac{b}{2}(\chi _0^{}+\chi _0^2)\right]=\chi _0(\gamma )\left(1\frac{b\omega }{2}\right)\frac{b\omega }{2}\frac{\chi _0^{}}{\chi _0},$$
(30)
so that

$$X_\omega (\gamma )=\left(1\frac{b\omega }{2}\right)\left[\mathrm{log}\left(\frac{\mathrm{\Gamma }(1\gamma )}{\mathrm{\Gamma }(1+\gamma )}\gamma \right)+2\psi (1)\left(\gamma \frac{1}{2}\right)\right]\frac{b\omega }{2}\mathrm{log}\frac{\chi _0(\gamma )}{\chi _0(\frac{1}{2})}.$$
(31)
As a consequence, the large $`t`$-anomalous dimension stays the same

$$b\omega t=\chi _\omega (\overline{\gamma })\frac{1}{\overline{\gamma }}.$$
(32)
but the saddle point exponent $`\mathrm{log}(\frac{G_{\mathrm{I}\mathrm{P}}^{(q)}}{N(\omega _{\mathrm{I}\mathrm{P}},\overline{t};t)})`$ in (2.2) changes

$$\begin{array}{c}t+\frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}\frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}\left\{\left[\mathrm{log}\left(\overline{\gamma }\frac{\mathrm{\Gamma }(1\overline{\gamma })}{\mathrm{\Gamma }(1+\gamma )}\right)+2\psi (1)(\overline{\gamma }\frac{1}{2})\right]\left(1\frac{b\omega _{\mathrm{I}\mathrm{P}}}{2}\right)\frac{b\omega _{\mathrm{I}\mathrm{P}}}{2}\mathrm{log}\frac{\chi _m(\overline{\gamma })}{\chi _m}\right\}\hfill \\ \hfill t+\frac{2}{b\omega _{\mathrm{I}\mathrm{P}}}\left[1+\psi (1)+\mathrm{log}b\omega _{\mathrm{I}\mathrm{P}}t\frac{b\omega _{\mathrm{I}\mathrm{P}}}{2}\mathrm{log}\chi _m+𝒪\left(\frac{1}{(b\omega _{\mathrm{I}\mathrm{P}}t)^2}\right)\right].\end{array}$$
(33)
It implies that approximately the eigenfunctions $`_{\omega _{\mathrm{I}\mathrm{P}}}^{(q)}(t)`$ get rescaled by the additional constant normalisation factor so that $`\sqrt{\chi _m}_{\omega _{\mathrm{I}\mathrm{P}}}^{(q)}_{\omega _{\mathrm{I}\mathrm{P}}}^{(k)}`$.


text:
# 1 Abstract

## 1 Abstract

Challenging stellar dynamical problems, such as the study of gravothermal oscillations in star clusters, have in the past initiated the very successful building of GRAPE special purpose computers. It is discussed, that present day tasks such as the formation and evolution of galactic nuclei with one or more massive black holes and the coupled stellar and gas dynamical processes in the formation of nuclei and star clusters, demand a new kind of hybrid architecture, using both GRAPE and a reconfigurable logics board called RACE. For such a system we have developed first implementations and floating point performance studies in the case of the SPH algorithm (smoothed particle hydrodynamics), which will be of great advantage for SPH modelling and also for direct $`N`$-body simulations using the more efficient Ahmad-Cohen neighbour schemes.

## 2 Historical Introduction

German-Japanese cooperation in stellar dynamics dates back for at least about two decades. One of its early highlights are the visit of D. Sugimoto as a visiting Gauß-Professor at the university observatory in Göttingen, Germany. At that collaboration, mainly with E. Bettwieser on the German side gravothermal oscillations of globular star clusters were detected using a gaseous model. The first author of this paper witnessed this as a freshman student in Göttingen. Subsequently a discussion arose whether such oscillations exist in real $`N`$-body systems; that question was a main motivation for the construction of GRAPE special purpose computers at the University of Tokyo in Japan. Their use made it possible to demonstrate that gravothermal oscillations are indeed present in real $`N`$-body systems. Subsequently GRAPE became an important computing equipment for many scientists working in stellar and galactic dynamics, also in Germany (see informations in www.astrogrape.org). In ongoing cooperative projects funded by the German and Japanese science foundations the first GRAPE computers began to work in Germany at the University of Kiel since 1993 (later moved to Astronomisches Rechen-Institut, Heidelberg) and at the two Max-Planck institutes (MPIA Heidelberg, MPA Garching). Presently, MPIA and ARI in Heidelberg, the University of Mannheim, Germany, and the University of Tokyo cooperate on the development of hardware and software suitable to tackle new challenges of modelling galactic nuclei, galaxy formation and evolution, and globular clusters. In the following a brief introduction is given into a selected sample of our astrophysical tasks. Finally some new ideas for the use of special purpose hardware are presented.

## 3 Gravothermal Star Clusters

Gravothermal Systems are those stellar systems in which the two-body relaxation has played an important role in their lifetime. It can approximately be modelled as a heat conducting gas. Astrophysical examples are dense (globular) star clusters and dense cusps in galactic nuclei around supermassive black holes. On the contrary many stellar systems such as galaxies as a whole are collisionless. Gravothermal systems are very difficult to study in direct numerical simulation, because most pairwise gravitational interactions have to be followed with high accuracy. The effect of gravothermal oscillations found by gas models could only 12 years later be seen in sufficiently large direct $`N`$-body simulations. The search for gravothermal oscillations was one of the driving ideas to build a fast special purpose computer called GRAPE (Gravity Pipe) . The dominant $`N^2`$ dependent part of the algorithm (pairwise force calculations) was mapped on it, using special arithmetics, pipelining and parallelisation. Still today we have not yet reached a final understanding of how globular clusters evolve; present studies try to improve theoretical modelling based mainly on the Fokker-Planck approximation and to compare its results with direct $`N`$-body models on rotating clusters, tidal fields and tidal shocks, a large number of primordial binaries, the influence of stellar evolution and the final fate of its remnants. The citation list is exemplary, but not exhaustive .

## 4 Galactic Nuclei with Black Holes

In the course of the merger of two galaxies their central black holes may ultimately coalesce, emitting in the very last phase strong gravitational radiation. During the early stages of the galaxy merger, the stellar component will form a dense core through violent relaxation within a rather short timescale. After that, the two supermassive black holes move through the stellar component with a velocity similar to the initial relative motion between the two galaxies. From this moment on, both massive bodies will feel dynamical friction. This friction leads the black holes to the newly-formed galactic center, while the frictional force becomes more efficient with increasing density. Through this process, the black holes must inevitably ‘find’ each other and form a binary system.

After being bound, the binary hardens (increases its binding energy) further through dynamical friction. At a certain binding energy dynamical friction between each black hole and the stellar system becomes inefficient. Close encounters and resonant three body interactions then provide further hardening for the binary. The most efficient process for binary hardening in this stage are those scatterings, after which the single stars gain very large velocities in a three body encounter with the black holes (superelastic scatterings). If the binary black hole centre of mass were fixed in the cluster, this process could evacuate the surroundings of the binary from suitable stars for further hardening, the hardening would stall. We observe in our simulations, however, strong recoils of the binary centre of mass due to the superelastic scatterings, so the hardening does not stall.

This is another challenging problem for our present-day software and hardware to simulate stellar systems. New hybrid codes including the recently developed parallel direct $`N`$-body integrators model a galactic nucleus containing two massive black holes using up to 128k single particles. With the recent GRAPE6 special purpose computer a model of three black holes in a nucleus of 512k particles has been simulated (Makino 2001, pers. communication). The shrinking of the black hole binary has to be followed until a phase where gravitational radiation sets in and a massive black hole merger occurs.

One of the most important parameters is the final eccentricity of the black hole binary, because the gravitational radiation induced merger depends critically on it. We present in Fig. 1 an example of the time evolution of the eccentricity of the black hole binary from the work of. A movie showing the black hole binary forming and moving in the galactic nucleus can be found under

ftp://ftp.ari.uni-heidelberg.de/pub/staff/marc/MPEG/simulation600.mpeg.

In the final hardening phase due to superelastic three body encounters the eccentricity of the black hole binary remains at a fairly large value (around 0.7), which is very interesting because it decreases the time scale for gravitational radiation merger of binary black holes dramatically, and thus increases our chances to detect gravitational radiation from such events via the planned gravitational radiation detector satellite LISA. In a system with three black holes Makino finds even larger eccentricities than 0.99.

## 5 Gas Dynamics, Galaxies and Star Clusters

While in certain evolutionary phases dynamics of star clusters and galactic nuclei is predominantly determined by stellar dynamical processes (they can be seen to first order as a system of gravitating point masses) their formation process is necessarily linked to a complex gas physics, including star formation, cooling, feedback processes and possibly even exotic events such as gas production by stellar collisions or star-gas interactions in dense star-gas systems. Again we have here for brevity selected a few exemplary references for the interested reader which cannot be exhaustive in any way. Many of the cited papers use a particle approach (smoothed particle hydrodynamics) to describe the dynamics of coupled gas-star systems; it is assumed that both stars and gas are represented by particles moving in a self-consistent gravitational potential, but the gas particles are subject to additional non-gravitational forces due to the gas physics.

## 6 Implementation Issues - GRACE

For astrophysical particle simulations including self-gravity, the determination of the full gravitational potential at each particles position is usually the most expensive step in terms of computational time required. This step, which is very efficiently done by the special hardware GRAPE, requires computational time $`T=\alpha N+\beta N^2`$, where $`\beta `$ is a time constant depending on hardware and on the algorithm (how many flops per individual pairwise force calculation). The first term linearly dependant on $`N`$ represents costs for advancing the particles and communication, if GRAPE is used. The idea of GRAPE is based on the strategy to reduce $`\beta `$ by hardware for the dominant term as much as possible. The most advanced $`N`$-body algorithms, however, use an Ahmad-Cohen scheme, where the short- and intermediate range force (neighbour or irregular force) is updated more often than the long range (regular) force. If on average for every $`\gamma `$ irregular force computations one regular force is computed the algorithm scales as $`T=\alpha N+\delta NN_n+\beta N^2/\gamma `$. Here $`N_n`$ is a typical neighbour number which is small compared to $`N`$. Such algorithm has been well implemented on general purpose parallel computers, but it does not work efficiently on the present GRAPE special purpose hardware. It is interesting to note that the second scaling also holds for high-accuracy codes using the SPH algorithm. The difference to pure $`N`$-body is that $`\delta `$ becomes larger, since SPH forces require of order 100 flops per pair of neighbour particles instead of only about 20 for the gravitational neighbour force only. Also in many existing codes the leading term is $`N\mathrm{log}N`$, instead of $`N^2`$, because a TREE algorithm is used for the long range force. In the challenging applications mentioned above the high accuracy $`N^2`$ algorithm will be required at least for a large subset of particles for accuracy reasons in gravothermal gas-star systems. Also the direct method is more advantageous for parallelisation and the use on GRAPE computers. Whether the TREE or direct force computations are used or not, the newly proposed hardware concept will be useful in both cases. The development of suitable algorithms for very large particle numbers and keeping a certain high degree of accuracy for the gravothermal systems is subject of present and future work.

Here we propose a new hardware architecture called GRACE (GRAPE and RACE), which combines the very fast GRAPE special purpose hardware with another hardware based on reconfigurable logics RACE (intermediate speed), both coupled via a common PCI interface to a host workstation.

The Fig. 2 plots a theoretical estimate to show, how a well balanced performance of the three components GRAPE, RACE and host workstation yields best performance. Just as an example we use for that picture an assumed performance of 1 Tflops for the GRAPE subsystem, 5 Gflops for the reconfigurable RACE subsystem, and 50 Mflop for the host workstation. In presently ongoing work (Lienhart, G., Wetzstein, M., et al., in preparation) a full loop of the SPH algorithm has been implemented on the available logic resources. The existing implementation works on the present hardware with roundabout 5 Gflops in floating point performance, as used in the plot. The interested reader can find more information on the reconfigurable RACE board under

http://www-li5.ti.uni-mannheim.de/fpga/?race/

and on the SPH implementation on it under

http://www-li5.ti.uni-mannheim.de/fpga/?astro/. Apart from the pure speed-up of such a combined GRACE system as compared to a pure GRAPE system such a new hybrid architecture has strong further potential advantages. First, it can overcome the host bottleneck for all kinds of $`N`$-body-SPH applications. Second, it can make the efficient use of GRAPE systems possible for the more intelligent neighbour schemes in direct $`N`$-body models. Third it allows flexibility where it is required: for the short- and intermediate range forces (softening, non-gravitational forces due to gas or molecular dynamics). While the general idea appears very promising, the implementation in detail has still to be carefully examined, in particular the way how the different hardware components communicate and store data and neighbour lists, for example. This work opens a path for more integrated GRAPE and RACE systems in the future.

## 7 Acknowledgements

Support and warm hospitality is gratefully acknowledged through collaborative German-Japanese grants 446 JAP 113/18/0 (1993-1998 with D. Sugimoto), 446 JAP 113/109/0 (1995-1997, with J. Makino, K. Nomoto, S.D.M. White) and 446 JAP 113/18/0-2 (2000-2002, with J. Makino, A. Burkert, R. Männer). This research has been performed within SFB 439 “Galaxies in the Young Universe” (funded by DFG) at the University of Heidelberg. Computational Resources of the HLRS Stuttgart and NIC Jülich have been used partly.