text:
# References

Self-similar approach to market analysis
V.I. Yukalov

Bogolubov Laboratory of Theoretical Physics
Joint Institute for Nuclear Research, Dubna 141980, Russia

Abstract. A novel approach to analyzing time series generated by complex systems, such as markets, is presented. The basic idea of the approach is the Law of Self-Similar Evolution, according to which any complex system develops self-similarly. There always exist some internal laws governing the evolution of a system, say of a market, so that each of such systems possesses its own character regulating its behaviour. The problem is how to discover these hidden internal laws defining the system character. This problem can be solved by employing the Self-Similar Approximation Theory, which supplies the mathematical foundation for the Law of Self-Similar Evolution. In this report, the theoretical basis of the new approach to analyzing time series is formulated, with an accurate explanation of its principal points.

PACS. 01.75.+m Science and society - 02.30. Lt Sequences, series, and summability - 02.50.-r Probability theory, stochastic processes, and statistics

1. Law of self-similar evolution

There are two sides of the approach I am going to present here: philosophical and mathematical. The first side gives the general idea of why should we think in these or those terms, and the second side is supposed to answer how concretely could we realize the idea. It is natural to start with the formulation of the general idea.

The idea to be formulated is general since it may concern any complex system, for instance, a market, a firm, a society, a nation, species, a person, a brain, a heart, and so on. Shortly speaking, a complex system is such that cannot, at least at the present time, be completely described by any finite number of given equations. A complex system exists as an empirical organism, with its specific character and behaviour.

To understand the character of a system means to notice some rules governing its behaviour and, therefore, to be able to predict the latter. If the system is complex, such rules cannot be described by fixed equations. But one may always characterize the behaviour in terms of self-similarity, implying that there exist some specific features that, to some extent, are repeated during the evolution of the considered system. Thus, characterizing a person, we name the corresponding features that give us impression of what can be expected from him or her.

The evolution of each complex system is always governed by some laws. This is why there is a similarity in the behavioral facts of this system. Such a behavioural self-similarity is easily noticeable in the life of particular persons as well as in the history of societies, nations, and biological species. This observation suggests that there exists the Law of Self-Similar Evolution, according to which the main features of a complex system are preserved in the course of its evolution. And if one has grasped these features, he or she should have the feeling of what could be expected in future.

When the behavioral facts, related to the evolution of a complex system, can be associated with some quantities, one gets time series. These can be market prices or indices, measurable signals from a brain or from a heart, some data characterizing societies or species etc. The examples of time series are numerous, their consideration being usually based on statistical analysis and the construction of dynamic regression models \[1–5\]. These models provide a rather reasonable description of sufficiently stable dynamics with repeated events, like seasonal variations, but they fail in treating such unstable systems as markets. Financial markets are very difficult to predict. Perhaps the most central question in finance is under what circumstances is prediction possible at all?

The widespread opinion is that market prices fluctuate absolutely randomly and, thus, are unforecastable. This point of view is what is called the Efficient Markets Hypothesis, which can be traced back to Samuelson . According to this hypothesis, in an informationally efficient market, properly anticipated prices fluctuate randomly. This includes the assumption that prices fully incorporate the expectations and information of all market participants, that all market traders are identical in the sense that all of them possess the same and the whole information on a market, all traders having the same physical and mental abilities. Such a hypothesis that all investors are fully rational agents that instantaneously and correctly process all available information is clearly unrealistic. Moreover, people, thanks God, are never identical neither in their abilities nor in their wishes. They are not elementary quantum particles, but each of them is a very complex system by its own. There is also increasing empiric evidence that even the most competitive markets are not strictly efficient \[7–9\].

A word of caution is called for with respect to the meaning of the term ”random”, which is often confused with ”chaotic”. These two notions are rather different. A realistic process can be chaotic, at the same time, possessing properties of both regularity and randomness in different proportions, so that it is, in principle, possible to distinguish chaos from randomness, and even control or erase chaos and make predictions .

Dynamics of all realistic complex systems always exhibits some part of randomness, either due to internal reasons, specific for nonlinear dynamical systems, or caused by external stochastic noise. Therefore randomness is inavoidable to this or that extent and, for some time intervals, it can mask the existence of underlying tendencies and persisting trends. Nevertheless, such trends are to be noticeable on average, whether this concerns the evolution of financial markets or biological species. Persisting features retained during this evolution imply that the complex system evolves self-similarly. It is worth emphasizing that the Law of Self-Similar Evolution does not mean and in no sense requires that there should be displayed some rigidly fixed properties of the considered system, but it rather tells that there must exist some trends in evolution. The evolutional self-similarity is not a static notion requiring the occurrence of stationary fixed properties but it is a dynamic concept stating the persistence of underlying trends.

In this section, the Law of Self-Similar Evolution has been described as a philosophical category, in rather general and perhaps vague words. The corresponding mathematical realization will be presented in the following sections. The suggested approach is mathematically based on the Self-Similar Approximation Theory whose application to the analysis of asymptotic series is expounded in Sec. 2. The self-similar extrapolation of asymptotic series can be directly reformulated as a self-similar forecasting for time series, which is explained in Sec. 3. A method for evaluating the probabilities of self-similar patterns is developed in Sec. 4. These sections present, for the first time, the complete mathematical foundation of the self-similar approach to analyzing arbitrary time series. Specifications related to market time series are also discussed. Principal points of the approach are summarized in Sec. 5.

2. Self-similarity in asymptotic series

Assume that we are interested in finding a function $`f(x)`$ of a real variable $`x`$. Let this function be defined by so complicated equations that we are able to extract from them only asymptotic expansions in the vicinity of some point $`x_0`$. Without the loss of generality, the expansion point may be taken as zero, $`x_0=0`$. Let us have several such expansions,

$$f(x)f_k(x)(x0),$$
(1)
enumerated by the index $`k=0,1,2,\mathrm{}`$. The basic problem is what can be said about the value of $`f(x)`$ at finite $`x`$ if all we know are asymptotic expansions $`f_k(x)`$ in the vicinity of $`x0`$? This problem is constantly met in physics and applied mathematics, where it is often called the problem of function reconstruction or the problem of summation of asymptotic series .

An original general approach, named the Self-Similar Approximation Theory, for reconstructing functions from a set of their approximate expressions has been developed and successfully applied to various problems \[13–23\]. The name comes from the basic idea of the approach to present the passage between subsequent approximations as a self-similar transformation. More precisely, it was shown \[15–19\] that for an approximation sequence $`\{f_k(x)\}`$ it is possible to construct a cascade, that is a dynamical system with discrete time, whose trajectory is bijective to the sequence $`\{f_k\}`$, so that the sought function $`f(x)`$ corresponds to a fixed point of the cascade. If we treat the given complicated equations, together with a calculational algorithm, as a complex system generating a sequence $`\{f_k\}`$, then the latter is nothing but a prototype of a time series, the approximation number $`k`$ playing the role of time. Hence the self-similar approximation theory is the mathematical realization of the law of self-similar evolution.

Asymptotic expansions are usually presented as power-law series

$$f_k(x)=\underset{n=0}{\overset{k}{}}a_nx^{\alpha _n},$$
(2)
where $`\alpha _n`$ are arbitrary real numbers arranged in ascending order,

$$\alpha _n<\alpha _{n+1}(n=0,1,\mathrm{},k).$$
(3)
For the purpose of self-similar analysis, the presentation of the series (2) has to satisfy several general properties:

(i) Analysis should not depend on the choice of units for the variable. Hence the latter is to be taken in a dimensionless form.

(ii) Asymptotic expansions are to be reduced to a scale-invariant form. To this end, it is always possible to factor out the term

$$f_0(x)=a_0x^{\alpha _0}(a_00)$$
(4)
and to introduce the scale-invariant function

$$\phi _k(x)\frac{f_k(x)}{f_0(x)},$$
(5)
where $`x`$ is assumed to be dimensionless. Evidently, function (5) does not depend on the change of scales for $`f_k`$.

(iii) If the variable $`x`$ pertains to a finite interval, the latter is to be normalized to the unitary interval. So that everywhere in what follows, it is assumed that $`x[0,1]`$.

Note that if the variable pertains to an infinite interval, one may pursue different ways depending on what additional information on the behaviour of $`f(x)`$ at $`x\mathrm{}`$ is available. If there is no such information, it is possible to transform the infinite interval to the unitary one by means of the change of variables $`x^{}=x/(1+x),x=x^{}/(1x^{})`$, so that $`x^{}[0,1]`$.

For the scale-invariant function (5), we have

$$\phi _k(x)=\underset{n=0}{\overset{k}{}}b_nx^{\beta _n},$$
(6)
where $`x[0,1]`$ and

$$b_n\frac{a_n}{a_0},\beta _n\alpha _n\alpha _0\beta _0=0.$$
(7)
Clearly, $`\phi _0(x)=1`$. Also, all powers $`\beta _n`$ are positive, even if some $`\alpha _n`$ are negative, which follows from the ascending order (3). Recall that, by construction, the series (6) are assumed to be asymptotic, having sense only for $`x0`$, while for any finite $`x`$ sequence $`\{\phi _k(x)\}_{k=0}^{\mathrm{}}`$ diverges. It is possible to say that the latter sequence converges just at one point $`x=0`$.

To proceed further, we have to transform divergent series to a form that would have sense for finite $`x[0,1)`$. This can be done with the help of control functions , which can be introduced in different ways \[13–23\]. Dealing with asymptotic series, it is convenient to invoke the multiplicative power-law transformation \[20–22\] defined as

$$\mathrm{\Phi }_k(x,s)x^s\phi _k(x),$$
(8)
with the inverse transformation

$$\phi _k(x)=x^s\mathrm{\Phi }_k(x,s).$$
(9)
Since power laws are common in describing fractal objects, equation (8) may be called the fractal transformation. The transform (8), according to equation (6), is

$$\mathrm{\Phi }_k(x,s)=\underset{n=0}{\overset{k}{}}b_nx^{s+\beta _n}.$$
(10)
Here $`s=s(x)`$ is a control function, whose role is to make the series (10) meaningful for finite $`x`$. These series can be considered as an expansion in powers of the new variable $`x^s`$. As is evident, such series are asymptotic with respect to $`x^s0`$. The latter limit can be achieved if, instead of forcing $`x`$ to zero, we keep $`|x|<1`$ and setting $`s\mathrm{}`$. Thus, the series (10) can be treated as asymptotic with respect to $`s\mathrm{}`$ for all $`|x|<1`$. Now we may say that the sequence $`\{\mathrm{\Phi }_k(x,s)\}_{k=0}^{\mathrm{}}`$ converges for all $`|x|<1`$, provided that $`s\mathrm{}`$. In this way, we come to the natural choice of the control function $`s\mathrm{}`$. Recall that we are considering the case when no additional constraints are imposed on the behaviour of the sought function and all we know are its asymptotic expansions (1). In the intermediate expressions, the value of $`s`$ is assumed to be asymptotically large, and the actual limit $`s\mathrm{}`$ is to be taken after the inverse transformation (9).

Since our consideration here concerns functions, we need to define the property of functional self-similarity, which should not be confused with geometric self-similarity describing fractals. The notion of geometric self-similarity  is connected with the scaling of a variable, which is only a particular kind of the more general notion of the functional or group self-similarity \[15–19\]. To correctly define the latter, we need to introduce some notation. We define the expansion function $`x(\phi ,s)`$ by the equation

$$\mathrm{\Phi }_0(x,s)=\phi ,x=x(\phi ,s).$$
(11)
With the form (10), this gives $`x(\phi ,s)=\phi ^{1/s}`$ and $`\phi =x^s`$. Introduce the mapping

$$y_k(\phi ,s)\mathrm{\Phi }_k(x(\phi ,s),s).$$
(12)
Let $`\mathrm{\Phi }_k(x,s)`$ be real for all $`k=0,1,2,\mathrm{}`$ and all $`x[0,1]`$. From $`\mathrm{\Phi }_k(x,s)𝐑`$ it follows that $`y_k(\phi ,s)𝐑`$ for all $`k=0,1,2,\mathrm{}`$ and $`\phi 𝐑`$. Therefore the mapping (12) is an endomorphism on $`𝐑`$. It is this endomorphism that allows us to formulate the property of group self-similarity we need.

Our aim is to present the change of the endomorphism (12), when varying the approximation number $`k`$, as the evolution of $`y_k`$ with respect to the discrete time $`k`$. From the point of view of group theory, self-similarity is nothing but a semigroup property. The latter, for the evolution of $`y_k`$ with respect to $`k`$, reads $`y_{k+p}=y_ky_p`$. As follows from the definition (11), the unit element is $`y_0`$, since $`y_0(\phi ,s)=\phi `$. The family of endomorphisms, $`\{y_k|k=0,1,2,\mathrm{}\}`$, with the semigroup property forms a dynamical system in discrete time, called the cascade. The semigroup property, in terms of the notation (12), takes the form

$$y_{k+p}(\phi ,s)=y_k(y_p(\phi ,s),s).$$
(13)
Since in the accepted interpretation, we treat the sequence $`\{y_k(\phi ,s)\}`$ as a trajectory resulting from the evolution of the cascade $`\{y_k|k=0,1,2,\mathrm{}\}`$, the relation (13) may be called the evolutional self-similarity. As far as the corresponding semigroup property is natural for dynamical systems, the equation (13) may also be termed the dynamic self-similarity. This equation (13) is a necessary condition for the fastest-convergence criterion ; the cascade fixed point representing the sought function.

Following the general theory \[14–19\], the cascade $`\{y_k|k=0,1,2,\mathrm{}\}`$ can be embedded into a flow $`\{y_\tau |\tau 0\}`$, which is a dynamical system in continuous time. For the latter, one may write the Lie equation which is a differential equation of motion. The flow velocity field is defined, by means of the Euler discretization, as the cascade velocity

$$v_n(\phi ,s)y_n(\phi ,s)y_{n1}(\phi ,s)(n=1,2,\mathrm{},k).$$
(14)
From equations (10) - (12), it follows that

$$y_k(\phi ,s)=\underset{n=0}{\overset{k}{}}b_n\phi ^{1+\beta _n/s},$$
(15)
which results in

$$v_n(\phi ,s)=b_n\phi ^{1+\beta _n/s}.$$
(16)
The differential equation of motion can be integrated. The integration over the effective time goes from $`\tau =n`$ to $`\tau =n+\tau _n`$, with $`\tau _n`$ being the effective time required for reaching a fixed point after the $`n`$-th step. In this way,

$$_n^{n+\tau _n}𝑑t=\tau _n.$$
(17)
Therefore, the evolution integral \[15–19\] acquires the form

$$_{y_n}^{y_n^{}}\frac{d\phi }{v_n(\phi ,s)}=\tau _n.$$
(18)
This, with the cascade velocity (16), yields

$$\mathrm{\Phi }_n^{}=\left(\mathrm{\Phi }_{n1}^{\beta _n/s}\frac{\beta _n}{s}b_n\tau _n\right)^{s/\beta _n},$$

where $`\mathrm{\Phi }_n^{}=\mathrm{\Phi }_n^{}(x,s)y_n^{}(x^s,s)`$ and $`\mathrm{\Phi }_n=\mathrm{\Phi }_n(x,s)`$. In particular,

$$\mathrm{\Phi }_1^{}(x,s)=x^s\left(1\frac{\beta _1}{s}b_1\tau _1x^{\beta _1}\right)^{s/\beta _1}.$$

Returning back to the function $`\phi _k^{}(x)`$ by means of the inverse transformation (9), we have to take the limit $`s\mathrm{}`$. For example,

$$\phi _1^{}(x)\underset{s\mathrm{}}{lim}x^s\mathrm{\Phi }_1^{}(x,s),$$

which results in the first-order self-similar approximant

$$\phi _1^{}(x)=\mathrm{exp}\left(b_1\tau _1x^{\beta _1}\right).$$

To obtain higher approximations, we can accomplish such a renormalization procedure $`2k`$ times for each $`\mathrm{\Phi }_k^{}(x,s)`$, which was called  self-similar bootstrap. However the same result can be reached twice faster, by accomplishing $`k`$ renormalizations, in the following way. We may present the function (6) as

$$\phi _k(x)=1+b_1x^{\beta _1}\left(1+\frac{b_2}{b_1}x^{\beta _2\beta _1}\left(1+\frac{b_3}{b_2}x^{\beta _3\beta _2}\left(1+\mathrm{}\right)\right)\mathrm{}\right).$$

The latter can be written in the form $`\phi _k(x)=1+x_1`$, in which $`x_1`$ is expressed through $`x_2`$, and $`x_2`$ through $`x_3`$, and so on according to the rule

$$x_n=\frac{b_n}{b_{n1}}x^{\beta _n\beta _{n1}}\left(1+x_{n+1}\right),$$

where $`n=1,2,\mathrm{},k`$. Considering each $`x_n`$ as a small parameter, we need to accomplish $`k`$ times the first-order renormalization procedure described above. As a result, introducing the notation

$$c_n\frac{a_n}{a_{n1}}\tau _n,\nu _n\alpha _n\alpha _{n1}(n=1,2,\mathrm{},k),$$
(19)
we come to the $`k`$-order self-similar exponential approximant

$$\phi _k^{}(x)=\mathrm{exp}\left(c_1x^{\nu _1}\mathrm{exp}\left(c_2x^{\nu _2}\mathrm{}\mathrm{exp}\left(c_kx^{\nu _k}\right)\right)\mathrm{}\right).$$
(20)
This, for short, can also be named the $`k`$-order superexponential.

It is worth emphasizing that although the form (20) reminds the Euler nested exponentials , it is principally different from the latter. First of all, the Euler superexponentials are defined only for integer powers $`\alpha _n`$. Second, when one tries to sum power series by means of such continued exponentials, one fits the coefficients in the latter so that to reproduce those in the power series, as a result of which the constructed superexponentials have the same radius of convergence as the related power series on the real axis .

In our case, the self-similar exponential (20) contains the coefficients $`c_n`$ given by equation (19), where the effective renormalization time $`\tau _n`$ is yet undefined. The latter plays the role of a control function that is to be determined from additional conditions. One way could be to find $`\tau _n`$ from the fixed-point conditions  having the form of the minimal-difference criterion . Such fixed-point equations not always possess solutions. Here we suggest another, and more general way of defining the control functions $`\tau _n`$. This method is commonly employed for defining control functions in the optimal control theory. The idea is to construct a cost functional whose minimization yields the control functions of interest.

The quantity $`\tau _n`$ appears in the evolution integral (18), where it has the meaning of an effective time required for reaching, after the $`n`$-th step, a fixed point representing the sought function. One, clearly, would like to reach the answer as fast as possible. The minimal number of renormalization steps is, obviously, one. Therefore one would like that the total effective time $`n\tau _n`$ be also close to one. Hence, one should look for $`\tau _n`$ being close to $`1/n`$. At the same time, how fast one reaches the fixed point depends on the distance of the latter from the starting point. The distance that is passed during the time $`\tau _n`$, with a velocity $`v_n`$, can be evaluated as $`v_n\tau _n`$. Thus, we need to find a minimal time $`\tau _n`$, being close to $`1/n`$, corresponding to the minimal distance $`v_n\tau _n`$. This suggests us to construct the time-distance cost functional

$$F=\frac{1}{2}\underset{n}{}\left[\left(\tau _n\frac{1}{n}\right)^2+\left(v_n\tau _n\right)^2\right],$$
(21)
whose minimization with respect to the control function $`\tau _n`$ yields

$$\tau _n=\frac{1}{n(1+v_n^2)}.$$
(22)
The velocity $`v_n=v_n(x)`$ is to be understood as the image, in the domain of $`x`$, of the cascade velocity (16), with taking account of the inverse transformation (9), which gives $`v_n(x)x^sv_n(x^s,s)`$. From here, one has

$$v_n(x)=\phi _n(x)\phi _{n1}(x)=b_nx^{\beta _n},$$
(23)
which defines the control function $`\tau _n=\tau _n(x)`$ according to the expression (22).

In this way, we find the controllers

$$c_n(x)=\frac{a_n}{a_{n1}}\tau _n(x),\nu _n=\alpha _n\alpha _{n1},$$

$$\tau _n(x)=\frac{1}{n[1+v_n^2(x)]},v_n(x)=\frac{a_n}{a_0}x^{\alpha _n\alpha _0}.$$
(24)
Combining equations (5) and (20), we obtain the $`k`$-order self-similar exponential approximant

$$f_k^{}(x)=f_0(x)\mathrm{exp}\left(c_1x^{\nu _1}\mathrm{exp}\left(c_2x^{\nu _2}\mathrm{}\mathrm{exp}\left(c_kx^{\nu _k}\right)\right)\mathrm{}\right),$$
(25)
in which $`c_n=c_n(x)`$, and $`f_0(x)`$ is given by formula (4). The approximant (25) extrapolates the asymptotic series (2), valid only for $`x0`$, to the region of finite $`x[0,1)`$. The value of $`f(x)`$ at the point $`x=1`$ can be defined as the limit from the left, as $`x10`$. Therefore the superexponential (25) extrapolates the sought function $`f(x)`$ from asymptotically small $`x0`$ to the whole unitary interval $`0x1`$.

It is worth recalling that in the asymptotic series (2) the powers $`\alpha _n`$ were assumed to be arbitrary real numbers, with the sole requirement that they are arranged in the ascending order (3). Some, or even all, of these powers could be negative. If so, then the initial term (4), with $`\alpha _0<0`$, has the power-law form that has been so much discussed in literature. The asymptotic existence of power laws is well known in critical phenomena. The relevance of power laws has repeatedly been claimed to describe many natural phenomena, ranging from earthquakes  to different economic and financial distributions . Since such power laws are practically always asymptotic, their more general form should include corrections leading to the power-law series (2). Extrapolating these asymptotic series in the described self-similar way, one should come to the self-similar exponentials. The first-order approximation then results in a kind of a stretched exponential. The stretched exponential distributions describe many phenomena in nature and economy either not worse or even better than power-law distribution functions . More generally, the extrapolation of power laws should lead to the self-similar nested exponentials (25), whose structure evidently demonstrates the existence of many scales.

3. Self-similarity in time series

The technique of the self-similar extrapolation for asymptotic series can be reformulated as the method of forecasting for time series. It is necessary to call attention to the existence of several principal points in this reformulation. Overlooking these points would essentially restrict the applicability of the method. The correct general way of self-similar forecasting for time series is advanced below.

First of all, the same basic requirements that were imposed on asymptotic series are compulsory for time series:

(i) The measured quantity is to be presented in a scale-invariant form. This can be easily done by normalizing the given data to the value $`f_0`$ of the measured quantity at the initial time, which is analogous to introducing the scale-invariant function (5).

(ii) The time variable has to be normalized to a dimensionless form, such that the prediction time would pertain to the unitary interval $`[0,1]`$. This normalization eliminates the ambiguity in defining the power of the power-law transformation (8), requiring that $`s0`$, which results in self-similar superexponentials (20).

In addition to these requirements, common both for asymptotic as well as for time series, there arises a question specific for time series: how the latter should be presented, as a backward or forward recursion? The answer to this question follows from the comparison of asymptotic series with time series. In order that the extrapolation of asymptotic series could be directly extended to forecasting for time series, the prediction horizon is to be a unitary interval and the available information from the past has to provide approximate forecasts for asymptotically small time $`t+0`$. This requirement can be strictly accomplished only for backward recursion. Hence, we have one more restriction, specific for time series:

(iii) For the correct usage of the self-similar analysis, time series are to be arranged as a backward recursion. This means that, if we are given a set $`\{f_n\}`$ of data $`f_n`$ corresponding to a quantity of interest, measured at the times $`t_n`$, where $`n=0,1,2,\mathrm{},k`$, the moments of time are to be ordered so that

$$t_{n+1}<t_n(n=0,1,2,\mathrm{},k),$$
(26)
with the initial time $`t_0=0`$. The past-history data base of $`k`$-th order is the set

$$𝐃_k=\left\{f_k,f_{k1},\mathrm{},f_0\right|t_k<t_{k1}<\mathrm{}<0\}.$$
(27)
To start the procedure of self-similar forecasting for future times $`t[0,1]`$, we need to possess a sequence of approximations valid for asymptotically small $`t+0`$. The role of such asymptotic forecasts can naturally be played by functions $`f_k(t)`$ interpolating the given data base (27) for the past time horizon $`t_kt0`$, so that

$$f_k(t_n)=f_n(n=0,1,\mathrm{},k).$$
(28)
This interpolation can be uniquely defined by the Lagrange interpolation formula  presenting the interpolation function

$$f_k(t)=\underset{n=0}{\overset{k}{}}f_nl_n^k(t)(k1)$$
(29)
as a series over the Lagrange polynomials

$$l_n^k(t)\underset{m(n)}{\overset{k}{}}\frac{tt_m}{t_nt_m}(nk).$$
(30)
Because of the property $`l_n^k(t_m)=\delta _{mn}`$ of the Lagrange polynomials, condition (28) is automatically satisfied. Note that the interpolation form (29) is evidently scale invariant with respect to time and the ratio $`f_k(t)/f_0`$ is scale invariant with respect to units of $`f_k`$. Also, one may remark that the time moments $`t_n`$ are not necessarily equidistant, but can be chosen arbitrarily.

The interpolation formula (29) can be rewritten as the algebraic polynomial

$$f_k(t)=\underset{n=0}{\overset{k}{}}a_nt^n(a_0=f_0),$$
(31)
with the coefficients $`a_na_{nk}`$ immediately following from equations (29) and (30). The values of these coefficients depend, of course, on the given data base (27), but for the simplicity of notation, we shall omit in what follows the additional index $`k`$. Consider the subsequence $`\{f_n(t)\}_{n=0}^k`$ of the terms

$$f_n(t)\underset{m=0}{\overset{n}{}}a_mt^m,$$
(32)
which, for a given $`k`$, tends to the polynomial (31) as $`nk`$ and for which

$$f_n(0)=a_0=f_0(n=0,1,\mathrm{},k).$$
(33)
Then the subsequence $`\{f_n(t)\}_{n=0}^k`$ can be treated as a sequence of approximations asymptotically valid for $`t+0`$. Employing the self-similar extrapolation, described in the previous section, we may construct a forecast for the time interval $`0t1`$. The result is given by formulas (24) and (25), where we need to change the variable $`x`$ by the time $`t`$ and to put $`\alpha _n=n`$, $`\nu _n=1`$. Instead of equations (24), we have the controllers

$$c_n(t)=\frac{a_n}{a_{n1}}\tau _n(t)(n=1,2,\mathrm{},k)$$

$$\tau _n=\frac{1}{n[1+v_n^2(t)]},v_n(t)=\frac{a_n}{f_0}t^n.$$
(34)
And the self-similar exponential (25) takes the form

$$f_k^{}(t)=f_0\mathrm{exp}\left(c_1t\mathrm{exp}\left(c_2t\mathrm{}\mathrm{exp}\left(c_kt\right)\right)\mathrm{}\right),$$
(35)
in which $`c_n=c_n(t)`$. Expression (35) is the self-similar forecast for the future time interval $`0t1`$, predicted on the grounds of the past-history data base (27).

4. Ensemble of possible scenarios

For each given data base (27), one may construct the self-similar forecast (35). But one can take several different data bases by varying either the data-base order $`k`$, or by changing the time intervals between the moments of time $`t_n`$, where $`n=0,1,\mathrm{},k`$, or by varying both the data-base order as well as the data-base scale. So that, in general, one can consider an ensemble of different data bases. Each of the latter has to be labelled by two indices,

$$𝐃_k(j)=\left\{f_k^{(j)},f_{k1}^{(j)},\mathrm{},f_0\right|t_k^{(j)}<t_{k1}^{(j)}<\mathrm{}<0\},$$
(36)
one index, $`k`$, defining the data base order, and another, $`j`$, specifying the chosen time scale of the past. For each data base (36), one obtains a forecast $`f_k^{}(j,t)`$ according to the rule (35), but with different values of $`c_n=c_{nk}`$. Hence there exists an ensemble $`\{f_k^{}(j,t)\}`$ of possible forecasts, or admissible scenarios. Which of these possible forecasts should one trust?

Life teaches us that in the majority of cases nothing can be trusted for hundred percent. But, when there can happen several different events, they can be classified by estimating their probabilities. Hence we need to define a probability measure on the ensemble of scenarios $`\{f_k^{}(j,t)\}`$.

This problem is analogous to the problem of pattern selection occurring for nonlinear differential equations in partial derivatives. Such equations sometimes possess a set of solutions corresponding to different spatio-temporal structures, or patterns . A general approach for treating the problem of pattern selection has been suggested . This approach can be directly applied for weighting possible scenarios from the given ensemble of self-similar forecasts. For this purpose, we may consider the passage from $`f_k^{}(j,t)`$ to $`f_{k+1}^{}(j,t)`$ as the motion with respect to $`k`$. Then the map $`\{f_k^{}(j,t)|k=1,2\mathrm{}\}`$ is to be treated as the image of a dynamical system with discrete time $`k`$. The probability of a scenario $`f_k^{}(j,t)`$ can be defined as

$$p_k(j,t)=\frac{1}{Z_k(t)}\mathrm{exp}\left\{\mathrm{\Delta }S_k(j,t)\right\},$$
(37)
where $`Z_k(t)`$ is a normalization factor, being the sum

$$Z_k(t)\underset{j}{}\mathrm{exp}\left\{\mathrm{\Delta }S_k(j,t)\right\}$$

over the pattern indices $`j`$, and the entropy variation

$$\mathrm{\Delta }S_k(j,t)S_k(j,t)S_1(j,t)$$
(38)
shows the change of entropy with respect to the effective time $`k`$. The entropy of a dynamical system may be defined, by analogy with statistical systems, as the logarithm of an elementary phase volume , the latter, in our case, being $`|\delta f_k^{}(j,t)|`$. Thus the dynamical entropy is

$$S_k(j,t)\mathrm{ln}|\delta f_k^{}(j,t)|.$$
(39)
Then the entropy variation (38) becomes

$$\mathrm{\Delta }S_k(j,t)=\mathrm{ln}\left|\frac{\delta f_k^{}(j,t)}{\delta f_1^{}(j,t)}\right|.$$
(40)
With the notation for the mapping multiplier

$$m_k(j,t)\frac{\delta f_k^{}(j,t)}{\delta f_1^{}(j,t)}=\frac{f_k^{}(j,t)/t}{f_1^{}(j,t)/t},$$
(41)
the entropy variation (40) reduces to

$$\mathrm{\Delta }S_k(j,t)=\mathrm{ln}|m_k(j,t)|.$$
(42)
It is convenient to introduce the average multiplier $`\overline{m}_k(t)`$ by the relation

$$\frac{1}{|\overline{m}_k(t)|}\underset{j}{}\frac{1}{|m_k(j,t)|}.$$
(43)
Using equations (42) and (43), for the scenario probability (37), we have

$$p_k(j,t)=\left|\frac{\overline{m}_k(t)}{m_k(j,t)}\right|.$$
(44)
Note that the scenario probability (44) is normalized with respect to the summation over the pattern indices $`j`$ corresponding to different data-base scales. In a particular case of just one fixed scale, one could accomplish the normalization with respect to the summation over $`k`$, which would define the probability weights for a restricted data-base ensemble, as was postulated in reference . The derived scenario probability (44) concerns the general case of an arbitrary ensemble $`\{𝐃_k(j)\}`$ of the data bases (36). Being general, the approach of the present paper makes it possible to answer several principal questions.

One important question concerns the choice of the data-base order. More in detail, the problem is as follows. For a fixed time scale, labelled by $`j`$, we may analyze different data bases $`𝐃_k(j)`$, with varying $`k=1,2,3,\mathrm{}`$. Then, how many terms $`f_k^{(j)}`$ should we take? That is, when should we stop increasing $`k`$? The answer is straightforward: The data-base order $`k`$ has to be increased till we reach numerical convergence. More precisely, this means the following. Let us be satisfied by the results of an error $`\epsilon `$. Then we need to increase $`k`$ up to the saturation number $`N_j=N_j(\epsilon )`$ such that

$$|f_{k+n}^{}(j,t)f_k^{}(j,t)|<\epsilon $$
(45)
for $`kN_j`$, all $`n0`$, and $`t[0,1]`$. For this $`k=N_j`$, we have the saturated data base $`𝐃(j)𝐃_{N_j}(j)`$ and the related saturated forecast

$$f^{}(j,t)f_{N_j}^{}(j,t)$$
(46)
characterized by the saturated scenario probability

$$p(j,t)p_{N_j}(j,t).$$
(47)
Thus, varying the data-base scales, labelled by $`j`$, we get the saturated ensemble $`\{f^{}(j,t)\}`$ of scenarios (46), with the probability measure (47).

The most probable scenario from the ensemble $`\{f^{}(j,t)\}`$ is the forecast $`f^{}(j_0,t)`$ having the largest probability, such that

$$\underset{j}{\mathrm{max}}p(j,t)=p(j_0,t).$$
(48)
The latter, because of the form (44), is equivalent to the condition of the minimum for the absolute value of the multiplier $`m(j,t)m_{N_j}(j,t)`$, so that

$$\underset{j}{\mathrm{min}}|m(j,t)|=|m(j_0,t)|.$$
(49)
Having the probability measure (47), it is possible to define the average forecast, or the expected forecast

$$<f(t)>=\underset{j}{}p(j,t)f^{}(j,t).$$
(50)
The dispersion

$$\sigma ^2(t)<f^2(t)><f(t)>^2,$$

in the case of a market, describes the market volatility. The latter can also be characterized by the variance coefficient $`\kappa (t)\sigma (t)/<f(t)>`$.

The last question to be answered in order to have a completely self-consistent theory is how to choose the data-base time scale. In general, the moments of time $`t_n^{(j)}`$ should not be compulsory equidistant. However, for practical purpose, it looks more convenient to take them as such, defining the time step as $`\mathrm{\Delta }_jt_n^{(j)}t_{n+1}^{(j)}`$. It seems natural to start with the time scale $`\mathrm{\Delta }_0=1`$, which equals the prediction horizon. Then, one may decrease as well as increase the time step, for instance, according to the rule $`\mathrm{\Delta }_{2j}=2^j,\mathrm{\Delta }_{2j+1}=2^j`$, with $`j=0,1,2,\mathrm{}`$. One has to stop decreasing and increasing the data-base time scale at such $`j=j_{max}`$, when numerical convergence is reached. This implies that, for a given error $`\epsilon `$, one gets the inequality

$$|f^{}(j+m,t)f^{}(j,t)|<\epsilon $$

for $`jj_{max}`$, all $`m0`$, and $`t[0,1]`$.

5. Summary of main ideas

In this paper, a novel approach to analyzing time series has been presented. There are several principal points distinguishing this approach from the standard one. The aim of this paper has been to clearly describe these principal points forming the body of a self-consistent theory. Not yet all parts of this theory have been exploited in full for practical applications; the calculational work is in process. But some simplified versions of the approach have been illustrated by a number of examples for market time series . In this concluding section, I would like to emphasize again the main ideas the approach is based on.

To better stress the principal difference of the present approach from the standard way of analyzing time series, let us recall the basic idea of the latter: For a given set of data $`\{f_n|n=0,1,2,\mathrm{},k\}`$ one tries to invent a relation $`f_t=f(f_0,f_1,\mathrm{},f_k,t,\xi )`$ connecting the value $`f_t`$ at the moment of time $`t`$ with the past data. This relation can be in the form of an explicit function or in the form of a difference or differential equation, including a stochastic term $`\xi `$ modelling noise \[1–5\]. There are, to my mind, two principal deficiencies of such an approach. First, I think that no explicit equations, indifferently to how elaborated they are, can grasp all peculiarities of a realistic complex system, such as a market. Second, any given relations reflect only the past history, providing an interpolation for the learning historical period, while for predicting future one needs an extrapolation.

The self-similar approach is based not on attempts to invent a relation between the historical points but it tries to discover dynamic trends resulting in these points. For this purpose, instead of studying relations between points, it is necessary to analyze relations between sequences. Interpolative formulas are used here only as a starting step. The comparison of different interpolative expressions makes it possible to produce an extrapolation, that is, forecasting.

The basic philosophical idea of the self-similar analysis of complex systems is the Law of Self-Similar Evolution, formulated in Sec. 1. And the mathematical foundation is provided by the self-similar extrapolation of asymptotic series, described in Sec. 2. This extrapolation can be reformulated as forecasting for time series, as is done in Sec. 3. Since any forecasting can only be probabilistic, the way of defining the related probability measure is explained in Sec. 4. Some details of the self-similar analysis can be changed. For example, the number of points in a given data base can be reduced by replacing several neighbouring values $`f_n`$ by either arithmetic averages or by fitting the values on a large time interval with the help of simple splines . However, in the process of these kinds of averaging, some information on the considered time series will be inavoidably lost.

Another change could concern the definition of the effective control time $`\tau _n`$. The latter is to be defined from the minimization of a cost functional. As an example, the time-distance cost functional (21) was considered. But, in general, one could opt for another cost functional, depending on the available information and imposed constraints.

Also, it would be possible to deal not with the initially given data base but with some transforms of it. For instance, one could consider the set $`\{\mathrm{ln}f_n\}`$ instead of $`\{f_n\}`$. Or one could keep in mind a more elaborated transform, like the wavelet transforms, often employed for analyzing time series .

It is worth stressing once more that the notion of self-similarity exploited throughout the paper is understood here as the group self-similarity, which is a more general notion than the trivial geometric self-similarity one usually talks about in connection with fractals. In the latter case one assumes the existence of the scaling relation $`f(\lambda x)=\lambda ^\alpha f(x)`$ for the considered function $`f(x)`$. Such a relation, with a given boundary, or initial, condition $`f(x_0)=f_0`$, immediately results in the power-law function $`f(x)=f_0(x/x_0)^\alpha `$. One could consider a slightly more complicated scaling relation as $`f(\lambda x)=u(\lambda )f(x)`$, with a known function $`u(\lambda )`$. However again, with the given boundary condition, this immediately gives the answer $`f(x)=f_0/u(x_0/x)`$. All such scaling relations produce the considered function in an explicit form. Whereas the group self-similarity (13) provides an equation that is yet to be solved.

An important point is that the group self-similarity employed here, and which is the basis of the Self-Similar Approximation Theory \[13–23\], has to do not with a scaling of a variable but with the motion with respect to the effective time whose role is played by the approximation number. Here it is the motion on the manifold of approximants.

Finally, if we would like, digressing from mathematical foundations, to briefly conclude why the idea of the group self-similarity does work for extrapolating asymptotic series and forecasting time series, then we should return back to the Law of Self-Similar Evolution telling that all complex systems develop self-similarly, preserving their generic features in the course of their evolution. Such features may be not noticeable from the first glance but rather hidden somewhere in genes. The group self-similarity is a kind of genetic self-similarity. Remember also that the Lord created the man in a self-similar way .

Acknowledgenemts

I am grateful to M. Ausloos and E. Yukalova for the interest to my work and useful discussions.


format: text

text:
# Vacancy diffusion in the Cu(001) surface II: Random walk theory

## 1 Introduction

Diffusion is one of the most commonly observed stochastic processes. Random walks, the paths of diffusing particles, are among the first typical applications in probability theory textbooks. However, the diffusing particles, or their paths, are often invisible, and their effect can be observed only indirectly. This is the case for example in vacancy mediated tracer diffusion, when on a lattice of atoms a diffusing vacancy displaces atoms along its path. The process is observed by following a labeled (or tracer) atom, whose steps are slaved to the motion of the vacancy.

In this paper we treat theoretically the physical system of diffusing In atoms in a Cu(001) surface layer, described in detail in paper I . In section 2 we review the literature about tracer diffusion. We describe our discrete model in section 3, the classical tracer diffusion problem adapted to better suit our experimental system. Finally section 4 gives a simple continuum formulation whose solution allows one to fit the experiments without adjustable parameters.

We start with summarizing the experimental observations . On a Cu(001) surface with up to several hundred atomic spacing wide terraces, a few surface Cu atoms were substituted with In atoms. Then the area was periodically imaged with a scanning tunneling microscope (STM), and the position of the In atoms was followed in time.

To our initial surprise, the In atoms do not move in a typical diffusive way, but instead are stationary for some time, and then suddenly (unresolvable by STM) they jump to a nearby lattice site, often up to 5 or more nearest neighbor spacings away from their original position. These long jumps tend to happen at the same instant for different In atoms in the imaged area. To explain this, we found that the only mechanism permissible on physical grounds  is vacancy mediated tracer diffusion.

Based on Embedded Atom Model calculations  we expect  that surface vacancies at room temperature have a low concentration, of the order $`10^9`$ on Cu(001). They both are created and recombine at steps, and stay in the top layer of a terrace, as it is energetically unfavorable to dive deeper. From the typical jump rate of $`10^6`$ Hz and terrace widths of a few hundred atomic spacings, their lifetime is estimated to be at most of the order of milliseconds. (Vacancies reaching the middle of the terrace have the longest expected lifetime.) As the long jumps of In atoms are the effect of a single vacancy, this short lifetime explains why the dynamics of the jumps can not be resolved with STM (with imaging rate up to 10 Hz), while the long waiting time between jumps, of the order of 10 s, enables one to distinguish independent vacancies.

The measured waiting time between the jump events has an exponential distribution, supporting the explanation of the jumps as independent events. The distribution of the jump vectors, which has been measured, will be compared to the numerical model in this paper.

## 2 Tracer diffusion

The problem of vacancy mediated tracer diffusion was considered for a long time . It has been solved first for the simplest case , when the diffusion of the vacancy is unbiased (all diffusion barriers are equal; the tracer atom behaves in the same way as the other atoms), the lattice is two-dimensional and periodic or infinite. There is a single vacancy present, it takes a nearest neighbor step in a random direction at regular time intervals, and has an infinite lifetime, as there are no traps. The solution is constructed by separating the motion of the tracer and the vacancy. The correlation of the steps of the tracer atom is calculated from the probability that the vacancy returns to the tracer from a direction which is equal, perpendicular or opposite to its previous departure. The distribution of the tracer atom spreads with time, and on an infinite lattice for large times it approaches the following functional form:

$$P_t(𝐫)=\frac{2(\pi 1)}{\mathrm{log}t}\mathrm{K}_0\left(\frac{r}{[\mathrm{log}t/(4\pi (\pi 1))]^{1/2}}\right),$$
(1)
where $`\mathrm{K}_0`$ is the modified Bessel function, time $`t`$ measures the number of steps of the vacancy, and at $`t=0`$ the vacancy is near the tracer atom. The non-Gaussian shape of the spatial distribution function is typical for vacancy mediated diffusion. The average displacement of the tracer particle diverges with time. This is a direct consequence of the fact that two is the marginal dimension for the return probabilities in the random walk problem. For higher dimensional problems, the probability that the vacancy returns to the tracer particle is less than unity and the average displacement of the tracer particle remains finite. However, for dimensions equal to or smaller than two, the vacancy always returns to the tracer particle and as a consequence its displacement diverges.

The same problem has been solved in an alternative way for all dimensions . From this solution one can calculate the number of tracer-vacancy exchanges up to time $`t`$: in two dimensions its distribution is geometric, with mean $`(\mathrm{log}t)/\pi `$. The continuum version of this problem has been considered as well in the form of an infinite-order perturbation theory ; the solution matches the asymptotic form of the lattice model.

In a very recent study the lattice calculations were generalized to biased diffusion . The difference between the tracer atom and the substrate atoms was taken into account by having different vacancy-tracer and vacancy-substrate exchange probabilities, while the rate of vacancy steps was kept constant. Repulsive interaction reduces, moderately attractive interaction increases the spreading of the tracer distribution.

Although these exact solutions are closely related to our In/Cu(001) system, the differences, e.g. in boundary conditions and vacancy lifetime make a direct comparison with experiments impossible. For this purpose we develop a model of tracer diffusion which includes the essential properties of the experimental system.

## 3 Discrete model

In this section we describe a discrete model for In/Cu(001) diffusion, solve it numerically, and present the results.

Our model is defined on a two-dimensional simple square lattice of size $`l\times l`$, centered around the origin. This corresponds to the top layer of a terrace of the Cu surface, with borders representing steps. All sites but two are occupied by substrate atoms. At zero time the two remaining sites are the impurity (or tracer) indium atom, which we release at the origin, and a vacancy at position $`(1,0)`$. This corresponds to the situation immediately after the indium atom has changed places with the vacancy, e.g. for the first time.

The only allowed motion is the exchange of the vacancy with one of the neighboring atoms. The exchange rate depends on the local environment, i.e. on the relative position of the vacancy and the impurity atom. This takes into account the effect of the lattice stress induced by the tracer atom on the energy landscape observed by the vacancy. Each rate was assumed to be simply proporional to the Boltzmann factor $`e^{\frac{\mathrm{\Delta }E}{k_BT}}`$, where $`\mathrm{\Delta }`$E is the activation energy for the considered diffusion step and $`k_BT`$ is the thermal energy at temperature T. The diffusion barriers of the vacancy were calculated with the EAM method . We used the interatomic potentials proposed by Finnis and Sinclair , where the metallic cohesion is taken into account via the second moment approximation to the tight-binding model. The calculation was done on a $`15\times 15\times 15`$ slab of copper. The system with different relative positions of an In atom and a surface vacancy was fully relaxed until a threshold force of $`1.87510^{12}`$ N on the atoms. The barriers were taken at the point halfway between the stable sites. See figure 1 for the values for the barriers. Since the barrier differences are large compared to $`k_\mathrm{B}T`$, the difference in jump probabilities are extremely large (see below for typical values).

When the vacancy reaches the perimeter of the lattice, its random walk is terminated, corresponding to the physical process of its recombination at surface steps. During its lifetime, the vacancy displaces atoms along its path, many of them multiple times. Thus also the tracer atom can end up displaced from its original position at the time of recombination of the vacancy. Averaged over the random walks of many independent vacancies, this yields a probability distribution of the different displacement vectors that the tracer atom can make as a result of its encounter with a single vacancy. Due to the boundary conditions, introduced by the finite size of the lattice and due to the distribution of exchange rates, an analytic solution to this problem is no longer possible.

When treating the above model numerically, we separate the motion of the vacancy and the tracer atom, as has been performed also in some of the analytical treatments referred to in section 2. In our case of a finite lattice, this separation introduces an approximation, which is valid only if the tracer atom is relatively close to the middle of the lattice. First, we calculate the probabilities that the vacancy, released at one atomic spacing from the tracer, returns the first time to the tracer from equal ($`p_{\mathrm{eq}}`$), perpendicular ($`p_{\mathrm{perp}}`$) or opposite ($`p_{\mathrm{opp}}`$) directions; we also calculate the probability of its recombination ($`p_{\mathrm{rec}}`$) at the perimeter instead of returning to the tracer. Knowing these return and recombination probabilities, we turn to the motion of the tracer atom, which performs a biased random walk of finite length. The direction of each step with respect to the previous one, and the probability that this was the last step, are obtained from the return and the recombination probabilities. With this method we lose all temporal information about the random walk of the tracer, but the individual steps of the tracer atom are orders of magnitude too fast to be resolved with the STM, and therefore the STM observations are insensitive to this information anyway.

In practice, both the return probabilities and the motion of the tracer atom are obtained with direct evaluation of probabilities, which has better convergence properties than Monte-Carlo-type methods.

As an illustration of this enumeration method, let us consider the computation of the vacancy return probabilities. We assign a variable to each lattice site, which measures the probability that the vacancy after $`s`$ atomic steps is at that site, while it has not exchanged with the tracer yet. Initially all probabilities are zero except at (1,0) where it is unity: we release the vacancy from here. The boundary acts as a trap for the vacancy, as well as the site (representing the tracer) at the origin. For each atomic step of the vacancy we update the site variables parallel by distributing their probability to the four neighbors according to the respective exchange probabilities. As the probability flows into the trap at the origin, we record the cumulative flow in each of the four directions leading to that site, which gives the return probabilities at the end. This iteration converges fast, and the convergence can be measured by the sum of the probabilities still on the lattice. The other computation, the motion of the tracer atom, is similar but slightly more complex. In that computation we assign a variable to each incoming edge of each site, which measure the probability that after $`s`$ steps — each corresponding to a vacancy return — the tracer is at the given site and that it arrived from the given direction. In addition, each site has a variable which accumulates the probability that the tracer become immobile at that site. These probabilities for the tracer arrival and immobilization are updated iteratively according to the previously obtained vacancy return probabilities.

We first tested our model on the case of unbiased vacancy diffusion, which would correspond to infinite temperature. The shape of the tracer distribution was similar to the experimentally observed one, but to achieve quantitative agreement the only remaining parameter of the model — the lattice size $`l`$ — had to be tuned to astronomical sizes. This clearly shows that taking equal barriers is an unrealistic oversimplification for the In/Cu(001) system.

Using the EAM barriers, for typical parameters $`T=320`$ K and $`l=401`$ the values for the return probabilities were $`p_{\mathrm{eq}}=12.4\times 10^7`$, $`p_{\mathrm{perp}}=1.1\times 10^7`$, $`p_{\mathrm{opp}}=4.2\times 10^9`$, and the recombination probability $`p_{\mathrm{rec}}=1.1\times 10^8`$. These values depend weakly on $`l`$ \[e.g. the dependence of the mean square displacement, calculated later in this paper from the return probabilities, is logarithmic: $`r^2\mathrm{log}(l/l_0)`$\]. This is a consequence of the fact that two is the marginal dimension for the return problem of the random walker. In higher-dimensional space the vacancy does not necessarily return, the return probabilities are asymptotically independent of the lattice size, and the final distribution of the tracer — also independent of lattice size for large lattices — takes the following form, e.g. in three dimensions : $`p(𝐫)=C_1\delta (𝐫)+C_2r^1\mathrm{exp}(\mathrm{\Lambda }r)`$.

For the case of an indium impurity in a Cu(001) lattice, the diffusion barrier for a vacancy exchange with the indium atom is considerably lower than all other barriers. Therefore, in most cases the vacancy returns from the direction of its previous departure, and the individual moves of the tracer atom are strongly anti-correlated. Both the numerical and the theoretical treatment are simplified significantly if we do not have to follow the large number of ineffective “back and forth” exchanges of the vacancy and the tracer atom. For this purpose, consider the small probability $`ϵ`$ that the vacancy does not return from the same, equal direction. Thus $`p_{\mathrm{eq}}=1ϵ`$. In the case of indium in Cu(001) the EAM calculations yield $`ϵ=2.410^7`$ (see above). We now define

$`\widehat{p}_{\mathrm{perp}}`$ $`=`$ $`p_{\mathrm{perp}}/ϵ`$
$`\widehat{p}_{\mathrm{opp}}`$ $`=`$ $`p_{\mathrm{opp}}/ϵ`$ (2)
$`\widehat{p}_{\mathrm{rec}}`$ $`=`$ $`p_{\mathrm{rec}}/ϵ`$
and have

$$2\widehat{p}_{\mathrm{perp}}+\widehat{p}_{\mathrm{opp}}+\widehat{p}_{\mathrm{rec}}=1.$$
(3)
If we represent the quasi-bound state of the rapidly exchanging (on average $`1/ϵ`$ times) vacancy and impurity atom with the position of a bond of the original lattice, then the vacancy-tracer pair walks on the bonds of the original lattice. The pair steps to each of the four perpendicular bonds with probability $`\widehat{p}_{\mathrm{perp}}/2`$, and to each of the two parallel ones with $`\widehat{p}_{\mathrm{opp}}/2`$, see figure 2. (The factors $`1/2`$ reflect the probability for the vacancy to escape either at the right or the left side of the quasi-bound position, which is $`\frac{1}{2}`$ for both sides in the limit of vanishing $`ϵ`$. The advantage of this approach is twofold: the path of the tracer is made of fewer steps (beneficial for numerics), and the bond-to-bond steps are now independent (beneficial for theoretical treatment).

Using this, the tracer atom described as if it forms a pair with the vacancy on one of the bonds adjacent to its original site, walks on the bonds lattice, and at the end (which happens with probability $`\widehat{p}_{\mathrm{rec}}`$ after each step) it is released at either end of the last visited bond. Results for the probabilities of the different jump lengths (beginning-to-end vectors of these paths) are shown on figure 3. Note, that the model calculations in figure 3 contained no adjustable parameters.

A general advantage of numerical modeling is that we have access to quantities which are difficult or impossible to measure experimentally. One example in our simulation is the probability that a tracer atom had an encounter with a vacancy, but its net displacement was zero. The temperature dependence of this quantity is plotted on figure 4. Since the dependence is weak, the assumption in the experimental measurements to associate In-vacancy encounters with visible (non-zero) jumps of the In atom is justified.

For a given set of diffusion barriers, our model has two parameters: the temperature and the lattice size. When we compare results with experiments, both can be independently obtained, and in principle there are no adjustable parameters. For example the case of $`T=320`$ K (figure 3) — using the distance to the nearest step in the experiment as lattice size in the numerical model — gave a good match with the experiment, but for the measurements at other temperatures we had to adjust the lattice size to obtain good agreement. Although the best fit lattice size in some cases was a factor of 2-3 smaller than the measured distance to the nearest step, the change in the return probabilities was much smaller, as they depend logarithmically on the lattice size. Undetected defects — acting as trap for the vacancies — closer to the In atom than the nearest step could also be accounted for this difference.

## 4 Continuum model

In the previous section we solved numerically the version of the tracer diffusion problem which was relevant for the In/Cu(001) experiments. Although this already enables full comparison, the numerical solution has the disadvantage that it cannot be described with a few parameters. In this section we develop a simple continuum description, where the overall shape of the jump length distribution is described with a single parameter.

We use our previous results for the return and recombination probabilities of the vacancy, and consider the random walk of the tracer-vacancy pair on the bond lattice. Let $`\varrho (𝐫,n)`$ denote the probability that the tracer-vacancy pair is at position $`𝐫`$ and at instance $`n`$, where $`n`$ counts the number of steps the tracer-vacancy pair takes. Since the subsequent steps of the pair are independent, we can write an effective diffusion equation for the evolution of $`\varrho (𝐫,n)`$:

$$\frac{\varrho (𝐫,n)}{n}=D_{\mathrm{eff}}^2\varrho c\varrho .$$
(4)
The first term on the right hand side corresponds to the steps the pair takes on the bond lattice, here $`D_{\mathrm{eff}}`$ denotes the mean square displacement per step of the pair. The second term corresponds to the recombination of the vacancy111Note that since $`n`$ counts the number of steps of the vacancy-tracer pair on the bond lattice, the term $`c\varrho `$ does not imply that the vacancy can recombine at any lattice site — in fact it recombines at terrace steps, between subsequent returns to the In atom.; at this point the pair breaks up. In the continuum approximation for space and $`n`$, the solution for a Dirac-delta initial condition at the origin is

$$\varrho (𝐫,n)=\frac{1}{4\pi D_{\mathrm{eff}}n}\mathrm{exp}\left(\frac{r^2}{4D_{\mathrm{eff}}n}cn\right).$$
(5)
The final distribution of the tracer atom after the vacancy recombined is obtained by the integration of the loss term in Eq. (4):

$$p(𝐫)=_0^{\mathrm{}}c\varrho (𝐫,n)n=\frac{1}{2\pi }\frac{c}{D_{\mathrm{eff}}}\mathrm{K}_0\left(\frac{r}{\sqrt{D_{\mathrm{eff}}/c}}\right),$$
(6)
where $`\mathrm{K}_0`$ is the modified Bessel function of order 0. The functional form of this solution is similar to that of the infinite lattice system in section 2, in spite of the differences introduced by the finite vacancy lifetime and the different boundary conditions.

The mean square displacement $`r^2`$ is directly proportional to the square of the width of the Bessel solution, apart from lattice corrections. We can determine the width of the Bessel solution from the parameters $`D_{\mathrm{eff}}`$ and $`c`$, which are obtained from the return probabilities in the Appendix:

$$r^2\frac{D_{\mathrm{eff}}}{c}=\frac{\widehat{p}_{\mathrm{perp}}+\widehat{p}_{\mathrm{opp}}}{4\widehat{p}_{\mathrm{rec}}}.$$
(7)
This continuum solution is shown in figure 3. It closely follows the numerical solution of the model, even for relatively small distances from the origin, where one would expect stronger lattice effects.

## 5 Summary

In this paper we described a model for vacancy mediated tracer diffusion on a finite lattice with absorbing boundaries for the vacancy. These boundary conditions were appropriate to model the vacancy mediated diffusion of In atoms embedded in the top layer of Cu(001) terraces. In addition to the numerical solution of the discrete model, we set up a simple continuum formulation of the model. The spatial distribution of the tracer atom in the continuum solution has a modified Bessel function profile. This form of non-Gaussian distribution is typical for tracer diffusion assisted by other diffusing particles. In order to enable a quantitative comparison with the STM measurement of the In/Cu(001) system, we introduced modified vacancy diffusion rates near the In atom, calculated with the EAM method. The modified rates affect the width of the Bessel function, without changing the functional form of this characteristic distribution.

We thank L. Niesen and M. Rosu from the University of Groningen for providing us with the computer code for the EAM calculations. This work is part of the research program of the “Stichting voor Fundamenteel Onderzoek der Materie (FOM),” which is financially supported by the “Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO).”

## Appendix A Obtaining $`D_{\mathrm{eff}}`$ and $`c`$ from the return probabilities

In this Appendix we calculate the parameters of the continuum model, $`D_{\mathrm{eff}}`$ and $`c`$, from the return and recombination probabilities in Eq. (2). The probability $`P_0^{(n+1)}`$ that the vacancy-tracer pair is at site “0” (not necessarily the origin) after $`n+1`$ steps of the pair is contributed from sites $`1,2,\mathrm{},\mathrm{\hspace{0.17em}6}`$ (see figure 2a; no upper index refers to probabilities after $`n`$ steps):

$$P_0^{(n+1)}=\frac{\widehat{p}_{\mathrm{perp}}}{2}(P_1+P_2+P_3+P_4)+\frac{\widehat{p}_{\mathrm{opp}}}{2}(P_5+P_6).$$
(8)
During one step in $`n`$, the change in $`P_0`$ is (using Eq. 3)

$`P_0^{(n+1)}P_0=`$ $`{\displaystyle \frac{\widehat{p}_{\mathrm{perp}}}{2}}(P_1+P_2+P_3+P_44P_0)`$ (9)
$`+{\displaystyle \frac{\widehat{p}_{\mathrm{opp}}}{2}}(P_5+P_62P_0)\widehat{p}_{\mathrm{rec}}P_0.`$
The finite difference on the left hand side approximates the derivative, and the terms on the right hand side can be collected to form a lattice approximation to the Laplacian:

$$^2P=a_{\mathrm{perp}}^2(P_1+P_2+P_3+P_44P_0)$$
(10)
and

$$_x^2P=a_{\mathrm{opp}}^2(P_5+P_62P_0)\frac{1}{2}^2P,$$
(11)
where the distances $`a_{\mathrm{perp}}=1/\sqrt{2}`$ and $`a_{\mathrm{opp}}=1`$ are in units of nearest neighbor spacing of the atomic lattice. In the approximation of Eq. (11) we assumed on average equal second derivatives in all directions. After these substitutions

$$\frac{P_0}{n}=\underset{D_{\mathrm{eff}}}{\underset{}{\frac{\widehat{p}_{\mathrm{perp}}+\widehat{p}_{\mathrm{opp}}}{4}}}^2P\underset{c}{\underset{}{\widehat{p}_{\mathrm{rec}}}}P_0,$$
(12)
and the coefficients $`D_{\mathrm{eff}}`$ and $`c`$ can be read off.


text:
# 𝑆𝑈(2) Running Coupling Constant and Confinement in Minimal Coulomb and Landau GaugesTalk presented by A. Cucchieri.

## 1 INTRODUCTION

An essential step for understanding and extracting physical information from gauge theories is the elimination of redundant gauge degrees of freedom. This is usually done by choosing a representative on each orbit of gauge-related fields (gauge fixing). In Ref.  Gribov showed that the Coulomb and Landau gauge-fixing conditions do not fix the gauge fields uniquely, i.e. there are many gauge-equivalent configurations satisfying the Coulomb or Landau transversality condition. These Gribov copies do not affect perturbative calculations, but their elimination could play a crucial role for non-perturbative features of gauge theories, such as color confinement and hadronization.

In order to get rid of the problem of spurious gauge copies, Gribov restricted  the physical configuration space to the region $`\mathrm{\Omega }`$ of transverse configurations, for which the Faddeev-Popov operator is non-negative.111In the Coulomb case the transversality condition and the positiveness of the Faddeev-Popov operator are satisfied on each time slice.  This region is delimited by the first Gribov horizon, defined as the set of configurations for which the smallest, non-trivial eigenvalue of the Faddeev-Popov operator is zero. We now know that $`\mathrm{\Omega }`$ is not free of Gribov copies and that the physical configuration space has to be identified with the fundamental modular region . Nevertheless, the region $`\mathrm{\Omega }`$ is of interest in numerical simulations, since it is the space of configurations satisfying the usual lattice Coulomb or Landau gauge condition.

The restriction of the path integral, which defines the partition function, to the region $`\mathrm{\Omega }`$ implies a rigorous inequality  for the Fourier components of the gluon field. From this inequality, which is a consequence only of the positiveness of the Faddeev-Popov operator, it follows that the region $`\mathrm{\Omega }`$ is bounded by a certain ellipsoid $`E`$. This bound causes a strong suppression of the (unrenormalized) transverse gluon propagator $`D^{\mathrm{tr}}`$ in the infrared limit . More precisely, it was proven  that, in the infinite-volume limit, $`D^{\mathrm{tr}}`$ vanishes at zero momentum, although the rate of approach to $`0`$, as a function of the momentum or of $`L`$, was not established. This is in marked contrast to the divergence in the infrared limit of the free massless propagator.

Finally, because of entropy considerations , the Euclidean probability gets concentrated near the first Gribov horizon222This has been verified numerically in the minimal Landau gauge . where the inverse of the Faddeev-Popov matrix diverges. This causes an enhancement of the ghost propagator $`G(k)`$ in the infrared limit . This enhancement is a clear indication of a long-range effect in the theory that may result in color confinement.

The confinement scenario is particularly simple in the minimal Coulomb gauge where the ghost propagator determines directly the Coulomb interaction . In fact, in this case, confinement of color, i.e. the enhancement at long range of the color-Coulomb potential $`V(R)`$, is due to the enhancement of $`G(\stackrel{}{k})`$ at small momenta. At the same time, the disappearance of gluons from the physical spectrum is manifested by the suppression at $`\stackrel{}{k}=\stackrel{}{0}`$ of the propagator $`D_{ij}(\stackrel{}{k},k_4)`$ of 3-dimensionally transverse would-be physical gluons. Remarkably, $`V(R)`$ is the instantaneous part of the $`4`$-$`4`$ component of the gluon propagator $`D_{44}(\stackrel{}{x},t)`$, and is a renormalization-group-invariant quantity . Its Fourier transform $`\stackrel{~}{V}(\stackrel{}{k})`$ may serve to define the running coupling constant of QCD by considering $`x_0g_C^2(|\stackrel{}{k}|)=\stackrel{}{k}^2\stackrel{~}{V}(\stackrel{}{k})=g_0^2\stackrel{}{k}^2D_{44}(\stackrel{}{k})`$,where $`x_0=12N/(11N2N_f)`$ has been calculated in . Clearly, if the color-Coulomb potential $`V(R)`$ is governed by a string tension at large distances, i.e. $`\stackrel{~}{V}(\stackrel{}{k})`$ goes like $`1/\stackrel{}{k}^4`$ at small momenta, then $`g_C^2(\stackrel{}{k})1/\stackrel{}{k}^2`$ in the infrared limit.

Similarly, in Landau gauge, one can consider  the running coupling constant $`g_L^2(k)=g_0^2[k^2D(k)][k^2G(k)]^2`$. This is also a renormalization-group-invariant quantity since (in Landau gauge) $`Z_gZ_3^{1/2}\stackrel{~}{Z_3}=\stackrel{~}{Z_1}=\mathrm{\hspace{0.17em}1}`$. In this case we obtain $`g_L^2(k)k^2`$ if, for example, $`D^{tr}(k)\text{const}`$ and $`G(k)k^4`$ in the infrared limit. On the contrary, if the gluon progator goes to $`0`$ in the infrared limit and the ghost propagator blows up not faster than $`k^4`$ then $`g_L^2(k)`$ has an infrared fixed point .

We test these theoretical predictions with data from a numerical study of $`SU(2)`$ lattice gauge theory, without quarks, in the minimal Coulomb and Landau gauge at $`\beta =2.2`$. Simulations were done at different lattice volumes $`L^4`$, in order to check for finite-size effects and, if possible, to extrapolate to infinite lattice volume. Details of notation and numerical simulations are given, for the Coulomb case, in  and will be presented, for the Landau case, in . For these simulations we have used a cluster of ALPHA work-stations (Coulomb and Landau data) and a PC cluster (Landau data) at the Dept. of Mathematical Physics of the University of São Paulo (DFMA/USP).333We thank Jorge L. deLyra for kindly providing us with access to these clusters. In the Landau case we have also used a PC cluster44416 nodes and a server with 866MHz Pentium III and 256/512 MB RAM memory, operating with Debian Linux. at the Institute of Physics of the University of São Paulo, São Carlos (IFSC/USP).

## 2 RESULTS

Our data in the minimal Coulomb gauge  clearly show (see Fig. 1) that the equal-time transverse gluon propagator $`D^{\mathrm{tr}}(\stackrel{}{k},L)`$ passes through a maximum and decreases as the momentum $`\stackrel{}{k}`$ approaches $`\stackrel{}{0}`$ (for fixed $`L`$), and that $`\stackrel{~}{V}(\stackrel{}{k})`$ is more singular than $`1/\stackrel{}{k}^2`$ at low $`\stackrel{}{k}`$, which indeed corresponds to a long-range color-Coulomb potential. We have also obtained  an excellent 2-pole fit for $`D^{\mathrm{tr}}(\stackrel{}{k},L)`$. Our fit indicates that the poles occur at complex $`m^2=x(L)\pm iy(L)`$. In the infinite-volume limit we have $`m^2=0\pm iy`$, for $`y=0.375\pm 0.162`$ in lattice units, or $`y=0.330\pm 0.142`$ GeV2 for the location of the gluon poles in $`\stackrel{}{k}^2`$. (It follows from the Nielsen identities  that these poles are independent of the gauge parameters.)

Similarly, in Landau gauge (see Fig. 2), the transverse gluon propagator $`D^{\mathrm{tr}}(k,L)`$ is suppressed in the infrared limit, while the ghost propagator $`G(k)`$ is more singular than $`1/k^2`$. Also here a 2-pole fit can probably be used to fit the data for $`D^{\mathrm{tr}}(k,L)`$. In this case, however, we need to improve the statistics in order to have better control over the fits, and we need to simulate at larger lattice volumes (we went up to $`26^4`$) to probe the infinite-volume limit.

Finally, for the running coupling constant $`g^2`$ we consider  the fitting formula $`k^2=\mathrm{\Lambda }^2\mathrm{exp}[(bg^2)^1][(bg^2)^r+z(bg^2)^\alpha ]^1,`$ which implicitly defines $`g^2`$. Here $`r=102/121`$, and $`\mathrm{\Lambda },b,z`$ and $`\alpha `$ are fitting parameters. For small $`g^2`$, which corresponds to large momenta, this formula is dominated by the first term in the denominator, whereas for large $`g^2`$, i.e. small momenta, it is dominated by the second term in the denominator. In particular, $`\alpha `$ governs the strength of the singularity of $`g^2`$ in the infrared limit; the case $`g^21/k^2`$ corresponds to $`\alpha =1`$.

In the Coulomb case we obtain  $`\alpha =1.9\pm 0.3`$; this corresponds to $`g_C^21/|\stackrel{}{k}|`$, and $`V(\stackrel{}{k})1/|\stackrel{}{k}|^3`$ at low momentum. In the Landau case we have $`\alpha =2.5\pm 0.3`$, which implies $`g_L^21/k^{0.8}`$. In both cases we have fitted the data with a low momentum cut-off at $`\stackrel{}{k}^2=0.5`$. Of course, an extrapolation in $`\beta `$ will be necessary to determine the strength of these singularities in the continuum limit.


text:
# FLAT DIRECTIONS IN SUSY GUTS 11footnote 1Talk given by Borut Bajc at SUSY’01, Dubna, Russia

## 1 Introduction

One of the main reasons to consider supersymmetry as a phenomenological viable model is the hierarchy problem: why is the Higgs boson so light, since it is naively believed that it should get a mass of the order of the next new scale, for example the grand unified scale or the Planck scale. Supersymmetry solves technically this problem, assuring that the boson-fermion cancellations in loop diagrams do not change considerably the scalar mass. There are however two remnants of the hierarchy problem still in susy models.

The first one is to understand, why mass scales in nature are so different, for example, why is $`M_{GUT}M_{Planck}`$, or why is $`M_WM_{GUT}`$.

The second one is how to achieve the mass splitting between Higgs SU(2) doublets and Higgs SU(3) triplets, which inevitably appear in any GUT.

This short paper deals with the above problems, as well as with cosmological implications. For a longer version see .

## 2 Flat directions

It has been shown long ago  that flat directions in supersymmetric models can generate dynamically very different scales. The idea of such a mechanism, known also as dimensional transmutation or radiative symmetry breaking, is very simple: a) have a supersymmetric model with at least one flat direction; this flat direction is exact at any order of perturbation thery due to a supersymmetric non-renormalization theorem; b) break supersymmetry spontaneosusly or softly, in order to lift this flat direction; higher loop contributions to the effective potential can generate a nontrivial minimum; since these corrections are logarithmic in the flat direction field, this new minimum is exponentially far away from the original scale (for example the susy breaking scale or the soft mass term). In such a way one can reproduce with soft supersymmetry breaking the electroweak  or the GUT  scale from the Planck scale.

The simplest example of a GUT flat diection is given by

$$W=\lambda Tr(\mathrm{\Sigma }^3),$$
(1)
where $`\mathrm{\Sigma }`$ is an adjoint of the gauge SU(6). The flat direction is

$$\mathrm{\Sigma }=\sigma diag(1,1,1,1,1,1).$$
(2)
Such simple flat directions are present only in SU(2n) theories.

The solution (2) with arbitrary nonzero $`\sigma `$ breaks the original SU(6) gauge group to SU(3)$`\times `$SU(3)$`\times `$U(1). To break further one can use for example the Fayet-Iliopoulos D-term of an extra (anomalous) gauge U(1):

$$D_{U(1)}=q_H|H|^2+q_{\overline{H}}|\overline{H}|^2\zeta ,$$
(3)
where $`H`$ and $`\overline{H}`$ are the fundamental and antifundamental representations of SU(6) and $`q_H`$, $`q_{\overline{H}}`$ have the same sign as $`\zeta `$. The supersymmetric solution

$$H=\overline{H}=\sqrt{\zeta /(q_H+q_{\overline{H}})}(0,0,0,0,0,1)$$
(4)
together with (2) clearly breaks the GUT into the standard model.

## 3 The doublet-triplet splitting

The previous example has unfortunately light Higgs triplets and heavy doublets, which is the opposite of what we want. To obtain the right pattern we need  at least two adjoints ($`A`$, $`A^{}`$) and two singlets ($`S`$, $`S^{}`$):

$$W=\lambda Tr(A^2A^{})+\lambda _SSTr(AA^{})+\lambda _S^{}S^{}Tr(A^2).$$
(5)
The solution is given by $`A`$ in the SU(4)$`\times `$SU(2)$`\times `$U(1) direction and proportional to an undetermined (flat direction) $`S`$ and a zero $`A^{}`$. Using again the D-term (3), one gets exactly two light Higgs doublets and heavy triplets .

## 4 Cosmology

At high enough temperature the $`T^4`$ term dominates and the symmetry restoration takes place 222 An opposite behaviour is however possible due to different initial conditions ., since a large nonzero vev of the flat direction would diminish the number of light degrees of freedom . Due to the flatness of the potential, the critical temperature  is only of order $`(m_{3/2}M_{GUT})^{1/2}`$, which is safe for the production of monopoles . Also, the phase transition is first order, reducing further the production rate. Of course, all this is relevant only if the phase transition really takes place, which is the subject of further studies.

## Acknowledgments

We thank the organizers of SUSY’01 for great hospitality and extremely well organized conference. We are grateful to Lotfi Boubekeur for collaboration in the early stage of this work and to Gia Dvali for important comments. The work of B.B. is supported by the Ministry of Education, Science and Sport of the Republic of Slovenia; the work of I.G and G.S. is partially supported by EEC under the TMR contracts ERBFMRX-CT960090 and HPRN-CT-2000-00152. Both B.B. and R.G. thank ICTP for hospitality during the course of this work.

## References


text:
# Spiral waves in accretion discs

## 1. Introduction

In close binaries, it is the detailed process of angular momentum transport that determines the structure of the accretion disc and therefore the rate at which gas, supplied from the mass donor, is actually accreted by the compact object. Although so fundamental to the process of accretion through discs, our understanding of angular momentum dispersal is very limited. We can roughly divide the possible physical mechanisms that may provide the required angular momentum transport into two classes. Those that work on a local scale and exchange angular momentum among neighbouring parcels in the disc, and those that rely on global, large scale structures in the disc. The local processes are commonly referred to as ‘viscous processes’ though it is clear that the molecular viscosity of the accretion disc material itself is many orders of magnitudes too small. Viscous interaction in the sheared Keplerian disc allows some material to spiral inwards, losing angular momentum, while excess momentum is carried outwards by other parts of the flow. In the famous $`\alpha `$-parameterisation of Shakura & Sunyaev (1973), this viscosity was replaced by a single dimensionless constant, which allows one to solve the structure equations for thin, viscously heated accretion discs.

A very different way of transporting the angular momentum is via density waves in the disc. In self-gravitating discs, the ability of density waves to transport angular momentum is a direct result from purely gravitational interaction between the wave and the surrounding disc material. Waves of this type can still transport angular momentum in the absence of self gravity, provided some mechanism exists that exchanges momentum between the wave and the fluid. Sawada, Matsuda, & Hachisu (1986) conducted numerical simulations of mass transfer via Roche lobe overflow of inviscid, non self-gravitating discs, and witnessed the development of strong spiral shocks in the disc which were responsible for the bulk of the angular momentum transport throughout the flow. Such trailing spiral patterns are the natural result of a tidal deformation of the disc that is sheared into a spiral pattern by the (near) Keplerian rotation profile of the disc material (Savonije, Papaloizou, & Lin 1994).

Although predicted in the 80s, observational evidence for spirals relies on the ability to spatially resolve the accretion discs in interacting binaries. Indirect imaging methods are thus required to search for such global disc asymmetries. Since the wave pattern is a co-rotating structure close to the orbital plane, Doppler tomography of strong emission lines is the ideal tool at hand. Although providing an image of the line emission distribution in velocity space, and not spatial coordinates, spirals should be readily identified as they maintain their spiral shape in the velocity coordinate frame.

The first convincing evidence of spirals was found in 1997, in the accretion disc of the eclipsing dwarf nova IP Peg. Since then, such waves have been confirmed in several outburst events of IP Peg as well as in other cataclysmic variables. In the next section, we will review the current observational status. A more detailed review can be found in Steeghs (2001).

## 2. Observational overview

### 2.1. IP Peg: when the story begun

When a Doppler map was constructed of the hydrogen and helium emission from IP Pegasi, just after rise to one of its outbursts, a surprising emission pattern was found (Steeghs, Harlaftis, & Horne, 1997, 1998). The accretion disc was far from symmetric but instead disc emission was dominated by a two armed pattern in the lower left and upper right quadrants of the Doppler map. The spiral arm velocities ranged between 500 and 700 km/s, corresponding to the outer regions of the accretion disc. The emissivity contrast between the spirals and other parts of the disc, was about a factor of $``$3 for H$`\alpha `$, and $``$5 in the case of HeI 6678 emission. There was no evidence for line emission from the bright spot. The spiral arm in the lower quadrant, closest to the secondary, extended over an angle of $``$100, and was weaker than the arm in the opposite quadrant. Strong emission from the irradiated secondary star was also present, producing the prominent S-wave at low velocities.

A week later, during the same outburst, more spectroscopy was obtained before the decline of the outburst had started (Steeghs et al. 1996). The Doppler image showed that the spiral pattern persists throughout the outburst, and the secondary star made a considerably larger contribution (from 6% to 10%) to the line flux. The arm in the upper right quadrant was still stronger, and the location of the arms had not changed, although the upper right arm appeared shorter.

In a different outburst Harlaftis et al. (1999) secured a whole orbit of IP Pegasi, three days after outburst. The HeII emission line at 4686Å was observed to provide a comparison with the previously observed H$`\alpha `$ and HeI emission patterns. The trailed spectrogram of HeII, again showed the familiar behaviour of the double peaks from the disc, leading to a two armed spiral in the Doppler tomogram. A very similar emission pattern from the disc was also present in the nearby Bowen blend consisting of CIII/NIII emission and the Helium I line at 4471Å. The two spirals could be traced for almost 180 degrees, and the velocity of the arms varied from 495 km/s to 780 km/s. The widths of the arms were significant, and even changed as a function of azimuth. It is thus an intrinsically broad feature: the arms cover a substantial part of the disc, with the strongest emission from the outer regions. The dynamics of the majority of the accretion disc material are thus affected by the presence of these spirals.

Morales-Rueda, Marsh & Billington (2000) present spectroscopy of the dwarf nova IP Pegasi taken during two consecutive nights, 5 and 6 days after the start of the August 1994 outburst. Doppler maps show marked spiral structure in the accretion disc. The spiral shocks are present on both nights with no diminution in strength from one night to the next.

The detections by Steeghs et al. (1997) and Harlaftis et al. (1999) were made $`1.5`$ and 3 days after the start of the August 1993 and the November 1996 outbursts respectively, whereas the data presented by Morales-Rueda et al. (2000) were taken 5 and 6 days after the outburst had started. The spectroscopic data taken 8 days after the start of the August 1993 outburst by Steeghs et al. (1996) also hints at their presence. This indicates that the spiral shocks are long-lived structures, with time-scale of the order of days instead of hours.

During the 1999 outburst of IP Peg, spirals were again observed (Steeghs & Boffin, in preparation). One can then conclude that Doppler tomography of IP Pegasi during outburst invariably shows the presence of spiral shaped disc asymmetries. The two armed spiral dominates the disc emission from the start of the outburst maximum and persists for at least 8 days, corresponding to about 50 binary orbits. The spirals are present in a range of emission lines from neutral hydrogen to ionised helium, the latter indicating that the gas concerned has to be hot, although it is not yet clear if we are looking at direct emission from the shock, or recombination emission from the spiral arms. The asymmetry between the two arms that was observed in the discovery data, is also present at other epochs and in other lines. Moreover, as we will now see, IP Peg is no longer the only system with evidence for distortions of a spiral nature.

### 2.2. Disc asymmetries in other systems

The eclipsing dwarf nova EX Dra, with a 5-hour orbit, was the second object to show similar disc behaviour. Joergens, Spruit, & Rutten (2000) present Doppler images of EX Dra during its July 1996 outburst, and the similarity of the Doppler maps with those of IP Pegasi is obvious. The disc seen in the HeI line reconstructed by Doppler tomography shows a clear two-armed spiral pattern pointing to spiral shocks in the disc. The Balmer and HeII maps also give evidence for the presence of spirals. The asymmetry of the spirals in the EX Dra disc is smaller than in the IP Peg observations perhaps indicating that the shocks are weaker. The spiral arm in the upper right is stronger in intensity as well as in asymmetry than that one in the lower left. This pattern is also seen in the Doppler maps of IP Peg.

Time resolved spectroscopic observations of U Gem during its March 2000 outburst show strong spiral shocks in the accretion disc (Groot 2001). The HeII 4686 most clearly shows the spiral arms during outburst. During the plateau at maximum brightness, 9 days into the outburst, the spiral shocks contribute $``$14% to the total HeII flux. The two arms of the spiral show a distinctly different evolution during the outburst and decline, which indicates an asymmetric evolution in the disc. Before outburst, the trailed spectrum shows a clear single S-wave which is identified as arising from the hot spot, located at a position of 0.55 R$`_{\mathrm{L}_1}`$. The upper arm is stronger than the lower one at the end of the maximum brightness plateau but not during the decline of the outburst. Nine days after maximum brightness, the slope in the radial velocity - orbital phase plane of the lower arm of the spiral is –2600$`\pm `$700 km s-1 orbit-1 and for the upper arm this is –1300$`\pm `$200 km s-1 orbit-1. One day later, these numbers are –1740$`\pm `$400 km s-1 orbit-1 and –1100$`\pm `$200 km s-1 orbit-1 for the lower and upper arms, respectively, while during the decline, they become, respectively, 250$`\pm `$140 km s-1 orbit-1 and 900$`\pm `$170 km s-1 orbit-1. Thus, the decrease of the slope of the lower arm is much more rapid and pronounced than for the upper one.

As stated by Groot (2001), a possible evolutionary scenario for the spiral shocks one can derive from the comparison of these three systems is that the shocks appear immediately when the outburst starts (IP Peg; Steeghs et al. 1997), but not just prior to outburst (U Gem in quiescence; Groot 2001), grow in strength when the outburst reaches its maximum magnitude (Harlaftis et al. 1999 and Joergens et al. 2000), and continue to gain in strength, or at least remain constant, during the plateau phase characteristic of many DN outbursts (Steeghs et al. 1996, Groot 2001), and then fade during outburst decline (Groot 2001).

In another outburst of U Gem, Steeghs, Morales-Rueda, & Marsh (in preparation) also observed clear spirals in the HeII at the start of the outburst, which persisted during the maximum brightness plateau, 4 to 5 days after outburst (Fig. 1). On the other hand, during decline, 9 days after outburst, the Balmer emission reveals a relatively circular disc.

Indication for the presence of spiral patterns in the stars SS Cyg and V347 Pup are also discussed in Steeghs (2001). Very recently, prominent spiral arms have also been observed in the early phases of the 2001 WZ Sge outburst (Kuulkers et al., this volume).

We shall now turn to the theory and see if such spirals, apparently ubiquitous in the discs of cataclysmic variables, can be explained.

## 3. Theoretical analysis

Spirals are an ubiquitous phenomenon in astrophysical sites. They may appear, for example, in galaxies, protostellar discs, associated with Jovian planets in proto-planetary discs as well as in accretion discs in binary systems. Such waves may play a very crucial role in the transport of angular momentum and hence in the process of accretion.

Tidally generated spiral waves are the result of tidal torques of the companion star on the orbiting disc material. Initially triggered at the outer edge of the disc, where the tidal interaction between the disc and the companion star is strongest, they take the form of trailing spirals because the azimuthal velocity of the disc material monotonically increases inwards. Lin & Papaloizou (1979) were among the first to observe the spiral pattern in their study of tidal torques on accretion discs in binary systems with extreme mass ratios. In their case, the spirals were the result of the interaction with the inner Lindblad 2:1 resonance. It was however Sawada, Matsuda, & Hachisu (1986a,b) who showed, in their 2D inviscid numerical simulations of accretion discs in a binary of unit mass ratio, that spiral shocks could form which propagate to very small radii. Spruit (1987) and later Larson (1990) made semi-analytical calculations which were followed by numerous - mostly 2D - numerical simulations. A more detailed account can be found in Boffin (2001).

Savonije, Papaloizou, & Lin (1994) presented both linear and non-linear calculations of the tidal interaction of an accretion disc in close binary systems. They showed, in agreement with Lin & Papaloizou (1979), that the resonance is significant over a region proportional to the temperature. Thus, even in CVs with larger mass ratios, they concluded that the centre of the 2:1 resonance can still be thought of as lying in the vicinity of the boundaries of the disc and that the resonance can still generate a substantial wave-like spiral response in the disc, but only if the disc is large and the Mach number is smaller than about 10. For larger Mach number, Savonije et al. consider that wave excitation and propagation becomes ineffective and unable to reach small radii at significant amplitude.

### 3.1. Opening angle

The opening angle of the spirals is directly related to the temperature of the disc as, when the shock is only of moderate strength, it roughly propagates at sound speed. For illustration purposes, we present in Fig. 2 two numerical simulations of an accretion disc using the Smoothed Particle Hydrodynamics (SPH) method. An adiabatic equation of state (eos) is used with $`\gamma =1.2`$ in one case, and a pseudo-isentropic eos with $`\gamma =1.01`$ in the other. As is clear from the figure, the $`\gamma =1.2`$ case produces a very open spiral, while the cooler one shows only a hint of disturbance at the outer edge, with very tightly wound spirals. This is even more clear in the high resolution simulation shown in Fig. 7 of Boffin (2001). Figure 2 also shows the reason for the difference: in the former case, the Mach number has values as small as 5, corresponding to a very hot disc, while in the latter case, the Mach number is above 20. This effect is the principal reason why it was believed that spirals may not be present in the accretion discs of cataclysmic variables, which are not hot enough. This is in particular the case in quiescent discs. As a byproduct, Fig. 2 also shows that the origin of the spirals lies in the tidal force of the companion, as in these simulations, no mass transfer is allowed from the companion.

### 3.2. Structure of spiral shocks

In this section, we would like to have a more detailed look at the way spirals affect the flow and hence lead to angular momentum transfer and mass accretion. To this aim, we will use the result of a simulation with $`\gamma =1.2`$ using a modified version of the SPH method (Boffin, Yukawa, & Matsuda, in preparation). Figure 3 shows the flow in a close binary with mass ratio $`q=0.15`$ when we allow mass transfer from the inner Lagrangian point. Again, the now very typical open spiral can be seen. We then make three slices in this flow, at angles 0, $`\pi /4`$ and $`\pi /2`$, considering only the upper right quadrant, and show the properties of the particles as a function of distance from a primary star (Fig. 5). As can be seen, the spiral shocks have a clear effect on the density field but also on the velocity field. In particular, the radial velocity exhibits a clear jump at the position of the shock. A more detailed analysis of the role of the spirals in the transport of angular momentum, in terms of an effective $`\alpha `$ can be found in Boffin (2001). In Fig. 4 we show the corresponding Doppler map, that is the structure of the flow in the velocity-velocity plane. This allows a direct comparison with observations.

### 3.3. Comparison with observations

Looking at Fig. 4 and 1, one sees that a qualitative agreement can easily be reached. In fact, in Fig. 4, we have also plotted a scaled version of the position and intensity of the spirals observed during the August 1999 outburst of IP Peg (Steeghs & Boffin, in preparation; see Steeghs 2001). The agreement between theory and observations seems rather convincing. This very fact is at the base of a lot of controversy in the literature, because the discs obtained in the numerical simulations are considered too hot to be realistic, even for systems in outburst (see e.g. Savonije et al. 1994, Godon, Livio, & Lubow 1998).

But what is the Mach number in the discs of cataclysmic variables ? A rough estimate is generally taken assuming Keplerian rotation and a central temperature of a few $`10^4`$ K. This gives then Mach numbers in the range 15 to 50. One has to note however that, because of the tidal field, the flow in the accretion disc is far from being Keplerian, and large departures exist that can affect the geometry of the spirals, especially in the outer parts. Moreover, the tidal force can increase the temperature of the these outer parts (Steeghs & Stehle 1999). These two effects join together to bring the Mach number to values smaller than about 10, at least in the outer parts of the disc. This could explain the apparent agreement between observations and theory. This explains also why spirals can only be seen during the outburst, when the disc is large and hot. However, it has still to be explored if other parameters than the Mach number do not play a role in the geometry of the spiral pattern. A definitive conclusion will only be possible once we can make detailed three-dimensional numerical simulations including a realistic equation of state. One should also caution that the spirals observed in the different cataclysmic variables are not always similar and that even for a given system, their shape may vary in the outburst cycle. As Groot (2001) noticed, its tomogram of U Gem nine days after maximum compares favorably with the low Mach number simulations in Steeghs & Stehle (1999), while the IP Peg, EX Dra and U Gem in decline tomograms, all appear to have their maximum emissivity rotated anti-clock wise with respect to the simulations of Steeghs & Stehle (1999). On the other hand, as shown by Fig. 4, it is also possible to have an agreement between numerical models and one outburst of IP Peg. Here again, one needs detailed self-consistent numerical simulations and a complete observational coverage of the same outburst of one cataclysmic variable to reach a more definitive conclusion. The time of pioneering is over, now comes the time for a detailed study.

## References

Boffin, H. M. J. 2001, in LNP Proc. 573, Astrotomography, ed. H. M. J. Boffin, D. Steeghs, & J. Cuypers (Springer), 69

Godon, P. 1997, ApJ, 480, 329

Godon, P., Livio, M., Lubow, S. 1998, MNRAS, 295, L11

Groot, P. 2001, ApJL, 551, L89

Harlaftis, E. T., Steeghs, D., Horne, K., Martin, E., T. M., 1999, MNRAS

Joergens, V., Spruit, H. C., Rutten, R. G. M. 2000, A&A, 356, 33

Larson, R. B. 1990, MNRAS, 243, 358

Lin, D.N.C., Papaloizou, J. 1979, MNRAS, 186, 799

Morales-Rueda, L., Marsh, T., Billington, I. 2000, MNRAS, 313, 454

Savonije, G.J., Papaloizou, J.C.B., Lin, D.N.C. 1994, MNRAS, 268, 13

Sawada, K., Matsuda, T., Hachisu, I. 1986, MNRAS, 219, 75

Sawada, K., Matsuda, T., Hachisu, I. 1986, MNRAS, 221, 679

Shakura, N. I., Sunyaev, R. A. 1973, A&A, 24, 337

Spruit, H. C. 1987, A&A, 184, 173

Steeghs, D. 2001, in LNP Proc. 573, Astrotomography, ed. H. M. J. Boffin, D. Steeghs, & J. Cuypers (Springer), 44

Steeghs, D., Stehle, R. 1999, MNRAS, 307, 99

Steeghs, D., Harlaftis, E. T., Horne, K. 1997, MNRAS, 290, L28

Steeghs, D., Harlaftis, E. T., Horne, K. 1998, MNRAS, 296, 463

Steeghs, D., Horne, K., Marsh, T. R., Donati, J. F., 1996 MNRAS, 281, 626