text:
# References

Fractal Diffusion in Smooth Dynamical Systems
with Virtual Invariant Curves
B.V. Chirikov111Email: chirikov@inp.nsk.su and V.V. Vecheslavov222Email: vecheslavov@inp.nsk.su
Budker Institute of Nuclear Physics
630090 Novosibirsk, Russia

## Abstract

Preliminary results of extensive numerical experiments with a family of simple models specified by the smooth canonical strongly chaotic 2D-map with global virtual invariant curves (VICs) are presented. We focus on the statistics of the diffusion rate $`D`$ of individual trajectories for various fixed values of the model perturbation parameters $`K`$ and $`d`$. Our previous conjecture on the fractal statistics determined by the critical structure of both the phase space and the motion is confirmed and studied in some detail. Particularly, we have found additional characteristics of what we termed earlier the VIC diffusion suppression which is related to a new very specific type of the critical structure. A surprising example of ergodic motion with a ”hidden” critical structure strongly affecting the diffusion rate was also encountered. At a weak perturbation ($`K1`$) we discovered a very peculiar diffusion regime with the diffusion rate $`D=K^2/3`$ as in the opposite limit of strong ($`K1`$) uncorrelated perturbation but, to the contrary, with strong correlations and for a very short time only. As yet, we have no definite explanation of such a controversial behavior.

1. Introduction: virtual invariant curves

In a 2D-map (2.1) we are going to study here the diffusion crucially depends on the global invariant curves (GICs) which cut the 2D phase space (a cylinder, see next Section) of the motion. Even a single such curve is sufficient to completely block the global diffusion (GD) in the action variable along the cylinder. As is well known by now the existence of the GICs does depend not only on the perturbation strength but also on its smoothness. It is convenient to characterized the latter by the temporal Fourier spectrum of the perturbation. For analytical perturbation the Fourier amplitudes decay exponentially fast. In this case GD sets up if the perturbation $`ϵ\stackrel{>}{}ϵ_{cr}`$ exceeds some critical value. Otherwise, the chaos remains localized within relatively narrow chaotic layers of nonlinear resonances. As a result, the GD is either completely blocked by GICs or the rate of the diffusion as well as the measure of its domain decay exponentially in parameter $`1/ϵ`$ as $`ϵ0`$ (the so-called Arnold diffusion, for a general review see, e.g., ).

By definition, the Hamiltonian of a smooth system has the power–law Fourier spectrum with a certain exponent $`\beta +1`$ (see, e.g.,  and references therein). In this case the GD is always blocked for some sufficiently small perturbation strength $`ϵ<ϵ_{cr}(\beta )`$ provided the smoothness parameter $`\beta >\beta _{cr}`$ exceeds the critical value. This is similar to the analytical Hamiltonian except that the critical perturbation depends now on the Hamiltonian smoothness ($`ϵ_{cr}(\beta )0`$ as $`\beta \beta _{cr}`$).

To the best of our knowledge, the most strong rigorous result reads: $`\beta _{cr}<4`$ for a 2D-map like in this paper (see ). However, a simple physical consideration  leads to even smaller value $`\beta _{cr}=3`$ which is still to be confirmed somehow, theoretically or numerically. In any event, the smoothness of our model here $`\beta =2`$ is even less.

Until recently, the behavior of dynamical systems in the opposite case $`\beta <\beta _{cr}`$ of a poor smoothness remained rather vague. Even though the most numerical data seemed to confirm the simplest behavior of some universal GD (see, e.g,, ) a few counterexamples were observed too (see, e.g., ). In the latter some trajectories remained within a certain restricted part of the phase space for a fairly long computation time. No clear explanation of these strange events has been given as yet.

Meanwhile, about 20 years ago (!) a number of mathematical studies revealed various possibilities for the existence of GICs in smooth systems with $`\beta <\beta _{cr}`$ (see, e.g., ). To us, the most comprehensive analysis of this problem was given by Bullett  who rigorously proved a strange survival of infinitely many GICs amid a strong local chaos. Surprisingly, all these interesting results remain essentially unknown, at least to physicists. Apparently, this is because the abovementioned mathematical papers were restricted (perforce !) to what could be done rigorously that is to the invariant curves only without any attempt to analyse the very interesting and important transport processes such as diffusion. This is still in reach to the physical analysis and numerical (or laboratory) experiments only. As a result, only after recent accidental rediscovery of GICs in chaos by Ovsyannikov  (which are still unpublished (!), see  for the full text of Ovsyannikov’s theorem) the intensive physical studies of this interesting phenomenon have begun .

Interestingly, both authors  made use of exactly the same model for which a strange locked-in trajectory was observed still much earlier . Apparently, this is because such a model (a particular case of our model with parameter $`d=1/2`$, see next Section) is the simplest one possessing those curious GICs (for discussion see ). Perhaps the main surprise was in that the GICs include the separatrices of nonlinear resonances which were always considered before as ones destroyed first by almost any perturbation. The principal difference is that now the invariant curves, separatrices including, do exist for the special values of the system parameters only (say $`K=K_m`$).

Even though there are infinitely many such special values of the parameter, and the infinitely many GICs for each of the value of which just a single GIC does completely block the GD the probability of the latter that is the measure of such $`K`$-values is apparently zero. So, a principal question to be answered is: what would be the behavior of that system for an arbitrary value of $`K`$? In  we conjectured that even though the set of $`K_m`$ is not generally everywhere dense  the density of this set is apparently rather high, so that one may expect some change (presumably suppression) of the diffusion for every $`K`$ value as compared to the ”usual” (familiar) dynamical system. In other words, we guessed that the structure of the phase space and of the motion therein may be changed by the formation of GIC at a close $`K`$ value even if at almost any $`K`$ there is actually no GIC. This is why we call now such a neighbour-$`K`$ invariant curve the virtual one (VIC) with respect to any $`K`$ .

Preliminary numerical experiments presented in  did confirm our conjecture. This was done by the prompt computation of the average diffusion rate $`D(K)`$ as a function of parameter $`K`$ in the domain with GICs, real or virtual ones. The experiments revealed a very strong suppression of the diffusion, up to many orders of magnitude, restricted only by the computation time. But what turned out to be even more interesting was a very complicated (apparently fractal) structure of the dependence $`D(K)`$. This seems to be a result of a very complicated structure of the model phase space itself. Preliminary, the latter looks like a so-called critical structure (see, e.g., ) but a rather specific one due to a forest of VICs.

In the present paper we begin studying this seemingly new type of such a structure. Specifically, we start with investigation of the statistical properties of the diffuson as one of the characteristic processes in the chaotic motion.

2. Model: the same again

For reader’s convenience we repeat below the discription of the model which is the same as in . The model is specified by a map in canonical variables action (momentum) $`p`$ \- phase $`x`$:

$$\overline{p}=p+Kf(x),\overline{x}=x+\overline{p}mod1.$$
$`(2.1)`$
where $`K=\epsilon >0`$ is perturbation strength (not necessarily weak), and ”force” $`f(x)`$ has a form of antisymmetric ($`f(y)=f(y),y=x1/2`$) piecewise linear ”saw” of period 1. The phase space of the model is a cylinder: $`0<x<1,\mathrm{}<p<+\mathrm{}`$.

As in  we are going to actually consider a family of maps with another parameter $`d`$ (see Fig.1 in ) and the force

$$f(x)=\{\begin{array}{cc}2x/(1d),\hfill & \text{if }|x|(1d)/2,\hfill \\ 2y/d,\hfill & \text{if }|y|d/2,\hfill \end{array}$$
$`(2.2)`$
where $`y=x1/2`$, and the second parameter $`d`$ ($`0d1`$) is the distance between the two ”teeth” of the saw $`|f(x)|=1`$ at points $`y=y_\pm =\pm d/2`$. The most studied particular case of the family corresponds to $`d=1/2`$ when the two-teeth saw $`f(x)`$ is symmetric. In the limit $`d=0`$ the two teeth merge in one, and all the invariant curves are destroyed. This was observed and explained in  for $`K>0`$. In the opposite case $`K<0`$ (which is equivalent to $`K>0,d=1`$) the dynamics of the model is completely different, and we will not consider it in this paper (for a brief discussion see ). In our 2D-map (2.1) the GIC supports rotation of phase $`x`$ around cylinder which bars any motion in $`p`$ over GICs. Unlike this the local invariant curve (LIC), surrounding, for example, a domain of regular motion (see, e.g.,  and Section 5 below), corresponds to oscillation in phase $`x`$ which allows other trajectories to bypass that obstacle.

The GICs, separatrices including, do exist in the whole interval $`0<d<1`$ but for special $`K`$ values only . Particularly, the invariant curves are completely absent  for sufficiently large parameter

$$K>K_B(d)=\frac{2d^2}{1+d},0<d<1.$$
$`(2.3)`$
If $`KK_B`$ (see below) the physical quantity of the main interest for us, the diffusion rate $`D`$, can be approximately calculated from the Fourier expansion of force (2.2) (for detailes see )

$$f(x)=\underset{n1}{}\frac{f_n}{n^\beta }\mathrm{sin}(2\pi nx),$$
$`(2.4)`$
where

$$f_n=\frac{2}{\pi ^2}\frac{\mathrm{cos}(n\pi )\mathrm{sin}(n\pi d)}{d(1d)},\beta =2.$$
$`(2.5)`$
Particularly, in the limit $`d=0`$

$$f_n=\frac{2}{\pi }\mathrm{cos}(n\pi ),\beta =1.$$
$`(2.6)`$
the smoothness parameter $`\beta `$ becomes less by one but the both values are less than critical $`\beta _{cr}=3`$.

The calculation of the diffusion rate and other quantities are done using the standard analysis of the nonlinear resonances and their interaction (overlap) (see, e.g., ). The calculation is especially simple if one can neglect the variation of coefficients $`|f_n|const`$ in (2.4). This simplification is exact for $`d=0`$ (2.6), and remains reasonably accurate  for

$$K\stackrel{>}{}\mathrm{\hspace{0.17em}3}K_B=\frac{6d^2}{1+d}$$
$`(2.7)`$
Then the diffusion rate is approximately given by a very simple standard relation

$$D(K)=\frac{\overline{(\mathrm{\Delta }p)_t^2}}{t}\frac{256}{\pi ^5}K^{5/2}0.57K^{5/2}$$
$`(2.8)`$
where $`t`$ is the motion time in map’s iterations, and parameter $`K1`$ is assumed to be sufficienly small. The latter expression in (2.8), which we will use below, is the result of extensive numerical experiments in  confirmed also in  for $`K\stackrel{<}{}0.1`$ (see  and next Section).

Notice that the dependence $`D(K)K^{5/2}`$ is different from the usual, or better to say, the simplest one $`D(K)K^2`$. This is explained by the dynamical correlation of motion which is determined by the frequency of the phase oscillation on nonlinear resonances

$$\mathrm{\Omega }_n=\sqrt{\frac{2\pi Kf_n}{n^{\beta 1}}}\mathrm{\hspace{0.17em}2}\sqrt{K}\mathrm{\Lambda }_n(K)1$$
$`(2.9)`$
Here $`\mathrm{\Lambda }_n`$ stands for the Lyapunov exponent charaterizing the local exponential instability of the motion which is the main criterion for dynamical chaos. Notice that for $`\beta =1`$ both $`\mathrm{\Omega }_n`$ and $`\mathrm{\Lambda }_n`$ do not depend on the Fourier harmonic number $`n`$. The exact value of the Lyapunov exponent in the limit $`d=0`$ is given by

$$\mathrm{\Lambda }=\mathrm{ln}(1+K+\sqrt{2K+K^2})\sqrt{2K}1$$
$`(2.10)`$
The latter expression is the approximation for small $`K`$ (cf. Eq.(2.9)) which is fairly well within the region of application of Eq.(2.8) ($`K\stackrel{<}{}0.1`$) with accuracy $`1\%`$. As the time in our model is discrete (the number of map’s iterations) both correlation characteristics, (2.9) and (2.10), must be small, hence the above restriction on parameter $`K`$ too.

In the opposite limit $`K1`$ the correlation between successive $`x`$ values are negligible, and one arrives at the ”usual” relation for the diffusion rate:

$$D(K)=K^2_0^1f^2(x)𝑑x=\frac{K^2}{3}$$
$`(2.11)`$
independent of the parameter $`d`$. In intermediate region ($`K1`$) the correlation causes the decaying oscillation (see ) which is beyond the scope of the present paper.

3. Diffusion without any invariant curves: averages and moments

As was already mentioned above there are no invariant curves for $`d=0`$. Moreover, the motion is ergodic that is of the simplest structure of the phase space (cf. Section 4 below). Therefore, this particular case is not of the main interest to us by itself. Nevertheless, it is a good introduction to our central problem considered in Section 6 below. A similar approach was taken in our previous paper .

First, we consider the time dependence of the diffusion rate $`D(K;t)`$. The semicolon here instead of usual comma is intended to emphasize that this time dependence is not a real physical contribution to the diffusion but rather a combination of two different processes: the proper diffusion via accumulation of random perturbation effects and a stationary regular oscillation of the diffusing variable ($`p`$ in our case) which is a sort of background for the diffusion. Such a phenomenon can be roughly represented by a simple relation:

$$D(K;t)D_{\mathrm{}}(K)+\frac{B(K)}{t}$$
$`(3.1)`$
where $`B(K)`$ is some function of the perturbation (see, e.g.,  and Eq.(3.5) below). In other words, in many cases, the present studies including, the nondiffusing stationary part can be separated from the diffusing part which much simplifies the analysis of this complicated process. All this can be describe, of course, via the standard method of the correlation of perturbation. However, this would lead to a much more intricate theoretical relations and, besides, to much less information on the diffusion dynamics (see, e.g., ).

An example of the diffusion kinetics is presented in Fig.1. The computation was done as follows. A number of trajectories $`M1`$ with random initial conditions homogeneously distributed within the unit area of the phase cylinder ($`0x_0<1,0p_0<1`$) were run for a sufficiently long time with successive outputs at certain intermediate moments of time $`t`$ as shown in Fig.1. Remember that $`t`$ is measured in the number of map’s iterations. Each output includes the diffusion rate $`<D>`$ averaged over all $`M`$ trajectories, and the dimensionless variance of that

$$V_M=\frac{<D^2><D>^2}{2<D>^2}$$
$`(3.2)`$
For Gaussian distribution of action $`p`$ this variance would be unity. This is indeed the case for a sufficiently long motion time when the measured diffusion rate reaches its asymptotic value $`D_{\mathrm{}}`$ (3.1). A quite different dependence $`V_M(t)`$ for previous smaller time is of no surprise (nor is it very interesting) as over there $`D(t)`$ depends on a completely different physical process one needs to pass over.

A real surprise was the very beginning of the diffusion, the plateau in Fig.1. This looks as a real diffusion unlike the following part of a stationary oscillation. Moreover, the diffusion rate on the plateau $`D_0=K^2/3`$ is the maximal one (2.11) as for large $`K1`$. Another interesting observation is the duration of this strange diffusion

$$t_0\frac{1}{\mathrm{\Lambda }}\frac{1}{\sqrt{2K}}$$
$`(3.3)`$
which is close to the inversed Lyapunov exponent, the rise time of the local exponential instability of the underlaying chaotic motion. The last but not least curious property is the fast increase of variance (3.2):

$$V_M(t)\frac{t}{3},2t\stackrel{<}{}t_0$$
$`(3.4)`$
as shown in Fig.1. This is qualitatively different from the behavior of the same diffusion rate for large $`K1`$ when the variance $`V_M1`$ is usual. The dynamical mechanism of this strange transitional diffusion is not completely clear and requires further studies. Apparently, it is related somehow to the main correlation (2.9) on the dynamical scale (3.3). Even though the initial ”diffusion” is relatively fast it goes on a short time only so that the relative change of the initial distribution of trajectories $`|\mathrm{\Delta }p|/|\mathrm{\Delta }p|_0\sqrt{D_0/\mathrm{\Lambda }}K^{3/4}1`$ is negligible for small $`K1`$ unless the initial distribution $`|\mathrm{\Delta }p|_0\stackrel{<}{}K^{3/4}`$ is very narrow. However, in the latter case the dependence $`D(t)`$ is very sensitive to the form of the initial distribution in $`p`$ as a few our preliminary numerical experiments reveal. The variance of $`D(t)`$ is especially strong for small $`tt_0`$ in the region of that mysterious plateau but eventually decays as $`t\mathrm{}`$ when the diffusion approaches its limit value $`D_{\mathrm{}}`$. Apparently, this is related to a complicated fine structure of the phase space and/or of the motion correlations. This interesting question certainly deserves further studies but in the present paper we consider the simplest, homogeneous, distribution of the trajectory initial conditions on the phase cylinder.

In this particular case, a very simple and surprisingly accurate empirical relation for the diffusion time dependence has been found starting from a qualitative picture (3.1). It has the form:

$$D(t)\frac{D_0+\tau D_{\mathrm{}}}{\left(1+\tau ^\gamma \right)^{1/\gamma }},\tau =c\mathrm{\Lambda }t$$
$`(3.5)`$
Here $`\tau `$ is the dimensionless time with an empirical fitting parameter $`c1`$ which is very close to one. The second empirical parameter $`\gamma 4`$ is less definite but it affects the turn of the dependence $`D(t)`$ at $`\tau 1`$ only. This relaxation of the diffusion rate has two time scales: (i) the plateau $`\tau _{pl}=1`$ or $`t_{pl}=1/c\mathrm{\Lambda }1/\sqrt{2K}1`$, and (ii) the relaxation $`\tau _R=D_0/D_{\mathrm{}}1/\sqrt{K}1`$ or $`t_R1/K`$ which is still much longer. Interestingly, the usual diffusion spreading of a very narrow initial $`p`$-distribution on the relaxation time scale

$$|\mathrm{\Delta }p|_R^2=D_{\mathrm{}}t_R=D_{\mathrm{}}(D_0/D_{\mathrm{}})/c\mathrm{\Lambda }=D_0/c\mathrm{\Lambda }=|\mathrm{\Delta }p|_{pl}^2$$

is exactly equal to the spreading on the plateau. Hence, the full relaxation spreading is twice as large which is also directly seen from the empirical relation (3.5):

$$|\mathrm{\Delta }p|_R^2=D(\tau _R)\frac{\tau _R}{c\mathrm{\Lambda }}\frac{D_0+D_0}{\left(1+\tau _R^\gamma \right)^{1/\gamma }}\frac{\tau _R}{c\mathrm{\Lambda }}K^{3/2}1$$

and which is still much less than the unit $`p`$-period.

In Fig.1 the empirical relation (3.5) is presented and compared to the numerical data in the dimensionless variables $`\tau `$ and $`D^{}=D/D_{\mathrm{}}`$ where $`D_{\mathrm{}}`$ is the asymptotic (”true”) diffusion rate (2.8). In these variables the curves with various $`K`$ values are similar, and converge in the limit $`\tau \mathrm{}`$.

Another interesting scaling can be done as follows. Let us calculate the diffusion rate $`D_{\mathrm{}}(D(\tau ))=D_{th}`$ from Eq.(3.5), and plot the ratio of that to the true rate (2.8):

$$\frac{D_{th}}{D_{\mathrm{}}}\frac{D(\tau )\left(1+\tau ^\gamma \right)^{1/\gamma }D_0}{\tau D_{\mathrm{}}}\mathrm{\hspace{0.17em}1}$$
$`(3.6)`$
Then, within the accuracy of scaling (3.5) and of fluctuations, this ratio must be always close to unity. This is indeed the case except the plateau ($`t\stackrel{<}{}t_0`$) where the rate $`D(\tau )`$ is almost independent of $`\tau `$ (see Fig.1).

The next important statistical property are the fluctuations of the diffusion rate. One characteristic of those is the dispersion of tracjectories which is characterized by the variance (3.2). If all the trajectories would be statistically independent the dispersion of the mean diffusion rate were

$$\left(\frac{\mathrm{\Delta }<D>}{<D>}\right)^2=\frac{2V_M}{M\mathrm{\hspace{0.17em}1}}$$
$`(3.7)`$
By construction, the trajectories are independent indeed with respect to their initial conditions but not necessarily to the corresponding diffusion rate. To check this we repeated the computation of diffusion $`N`$ times with new and independent initial conditions, and then calculated the second (new) dimensionless variance for the average diffusion rate:

$$V_N=\left(\frac{<<D>^2>_N}{<<D>>_N^2}\mathrm{\hspace{0.17em}1}\right)\frac{M\mathrm{\hspace{0.17em}1}}{2V_M}\mathrm{\hspace{0.17em}1}$$
$`(3.8)`$
Again, if Eq.(3.7) holds true the variance $`V_N`$ should be close to one.

The time dependence of both variances, $`V_M(t)`$ and $`V_N(t)`$, is shown in Fig.1. Remarkably, their behavior is qualitatively different. The first variance $`V_M(t)`$ depends on the distribution function of $`p`$ in the ensemble of trajectories while the second variance $`V_N(t)`$ is affected by the statistical dependence (or independence) among trajectories whatever their distribution function. The results of our numerical experiments presented in Fig.1 clearly demonstrate that the distribution in $`p`$ quickly deviates from the Gaussian one during the diffusion on plateau, and come back only in the limit $`t\mathrm{}`$ when the diffusion rate $`DD_{\mathrm{}}`$ approaches the asymptotic value without any nondiffusing part. Unlike this, the trajectories remain statistically independent during the whole process of the diffusion relaxation. We will come back to discussion of this interesting point in the conclusion to this paper (Section 7).

Now we turn to the most informative statistical characteristic, the distribution function $`f(D)`$ of the diffusion rate.

4. Diffusion without any invariant curves: the distribution function

In the main part of our paper (Section 6) we shall be primarily interested in the distribution tail $`D0`$ of very low diffusion rate. The shape of this tail is known to be an important characteristic of the critical structure of the motion (see, e.g., ). First indications of such a structure in the presence of the virtual invariant curves have been observed in . Here we continue these studies.

Since the statistics of the far tail is always rather poor we make use, as in , of a special version of the integral distribution

$$F(D)=_0^Df(D^{})𝑑D^{}\frac{j}{J}$$
$`(4.1)`$
the so-called ”rank-ordering statistics of extreme events” (see, e.g., ). To this end the following simple ordering of the $`D(j)`$ values (events) of the diffusion rate is sufficient: $`D(j+1)>D(j),j=1,2,\mathrm{},J`$. Then the integral probability is approximately given by the ratio $`j/J`$ as shown in Eq.(4.1).

In computation we typically ran $`M`$ trajectories by $`N`$ times (see Section 3), so that the maximal number of events reached $`J=M\times N=10^4\times 10=10^5`$. To obtain the lowest possible $`D`$ values and minimize, at the same time, a rather big output we ordered all the computed events but printed out much less of those $`J_0J`$ in such a way to get some all $`J_1<J_0`$ first (the smallest) $`D_j`$ while the rest were printed out in a logarithmic scale. An example of such distribution is presented in Fig.2 for $`K=0.001`$ in variables $`D^{}=D/<D>`$ and $`F(D^{})=j/J`$ where $`<D>`$ is some average diffusion rate (see below). The upper distribution corresponds to a rather long motion time $`t=10^41/K`$ when the mean diffusion rate is already very close to the limit $`D_{\mathrm{}}`$. For the lower distribution $`t=10`$ is very short and corresponds to the plateau.

At least in the former case when the $`p`$ distribution is Gaussian (see Section 3) the distribution

$$f(D)=\frac{\alpha ^\lambda }{\mathrm{\Gamma }(\lambda )}D^{\lambda 1}\mathrm{e}^{\alpha D}$$
$`(4.2)`$
is the so-called Pearson $`\mathrm{\Gamma }`$-distribution with the two moments

$$<D>=\frac{\lambda }{\alpha },(\mathrm{\Delta }D)^2=<D^2><D>^2=\frac{\lambda }{\alpha ^2}$$
$`(4.3)`$
which are the mean and variance, respectively. For Gaussian $`p`$-distribution the reduced variance (3.2) $`V_M=1`$ whence

$$\left(\frac{\mathrm{\Delta }D}{<D>}\right)^2=\frac{1}{\lambda }=\mathrm{\hspace{0.17em}2}$$
$`(4.4)`$
and $`\lambda =1/2`$ independent of $`\alpha `$. If, moreover, we introduce the dimensionless diffusion rate

$$DD^{}=\frac{D}{D_{\mathrm{}}}$$
$`(4.5)`$
with average $`<D^{}>=1`$, we obtain from Eq.(4.3) $`\alpha =\lambda =1/2`$, too. Then, the new distribution becomes

$$f(D^{})=\frac{(D^{})^{1/2}\mathrm{e}^{D^{}/2}}{\sqrt{2\pi }}$$

and

$$F(D^{})=_0^D^{}f(D^{})𝑑D^{}\sqrt{\frac{2}{\pi }D^{}}$$
$`(4.6)`$
where the latter expression gives the asymptotics $`D^{}0`$ we need. This asymptotics is in a very good agreement with empirical data in Fig.2 even at $`D^{}0.1`$ (!). For very small $`D^{}`$ the accuracy of the agreement is limited by the fluctuations due to a few remaining points. The smallest $`D^{}=8.3\times 10^{11}`$ corresponds to the estimate $`D_{min}^{}1/J^2=10^{10}`$.

Since the distribution $`f(D^{})`$ in (4.6) is also Gaussian of $`\sqrt{D^{}}`$ the integral $`F(D^{})`$ admits a very simple approximation found in :

$$F(D^{})\{\begin{array}{cc}1\frac{\mathrm{exp}(D^{}/2)}{\sqrt{D^{}}+\mathrm{\hspace{0.17em}1}},\hfill & D^{}>\mathrm{\hspace{0.17em}1}/2\hfill \\ \sqrt{\frac{2D^{}}{\pi }},\hfill & D^{}<\mathrm{\hspace{0.17em}1}/2\hfill \end{array}$$
$`(4.6a)`$
The relative accuracy $`|\mathrm{\Delta }F/F|<0.05`$ of that approximation is better than 5% in the whole range of $`F`$. Actually, the accuracy is even much better except a narrow interval at $`D^{}1/2`$.

Thus, the upper distribution in Fig.2, which describes the real diffusion at sufficiently long motion time, is well in agreement with the available theory. This is no longer the case for the lower distribution on the plateau. In itself, this is not a surprise as in the latter case, unlike the former one, the measured diffusion rate is mainly determined by nondiffusive processes. However, a very interesting feature of this nondiffusive distribution is in that the exponent of the power-law tail remains exactly the same as if the $`p`$-distribution were again a Gaussian one. The simplest explanation, quite plausible to us, is in that the far tail still represents such a distribution which is a part of the whole distribution according to our original picture expressed by estimate (3.1). One immediate inference would be decrease of the tail probability if we still use the same variable $`D^{}=D/D_{\mathrm{}}`$. This is indeed the case according to the data in Fig.2 !

A more difficult problem is the quantitative estimate of the distribution shift for the motion time $`t\stackrel{<}{}1/K`$ when the ratio $`<D^{}>=<D(t)>/D_{\mathrm{}}>1`$. This shift can be characterized either via the probability decrease by $`R_F`$ times for a fixed $`D^{}`$ or via the increase of $`D^{}`$ itself by $`R_D`$ times for a fixed probability. Notice that on the tail $`R_D=R_F^2`$ due to the square-root dependence (4.6). Characteristic $`R_D`$ seems to us more preferable since it describes the shift not only of the tail but also (qualitatively) of the whole distribution $`F(D^{})`$.

After some playing with the data we have found the following empirical relation for the tail shift:

$$R_D(D_{})D_{}^a$$
$`(4.7)`$
where the new diffusion ratio

$$D_{}(\tau )\frac{D_0}{\tau D_{\mathrm{}}}+\mathrm{\hspace{0.17em}1}$$
$`(4.8)`$
and the fitted exponent $`a=0.45`$.

The philosophy behind this relation is following. We start with our original picture of a combined diffusive/nondiffusive process (3.1) which is almost our final choice (4.8). However, at the beginning we seemed to improve the original relation by inclusion of our surprising discovery, the plateau. Specifically, we tried to make use of Eq.(3.5) which is in a good agreement with the empirical data for the dependence $`D(t)`$ (see Fig.1). Also, we have found that it partly described the distribution $`F(D)`$ too, except on that mysterious plateau ! Then, our final, so far, step was the change from (3.5) back to a version of (3.1) in the form (4.8).

How strange it may seem this did work with a reasonable accuracy as the insert in Fig.2 demonstrates. The remaining question ”why?” is still to be answered in farther studies. Actually, this is a general serious problem of the dynamical mechanism underlaying the plateau formation and statistics.

Our empirical relation (4.7) can be represented in a different way. Namely, instead describing the actual distribution tail shifted with respect to the limiting asymptotics (4.6) we may introduce the scaled diffusion rate

$$D\frac{D}{R_D}$$

whence

$$D^{}\frac{D^{}}{R_D}$$
$`(4.9)`$
The result is shown in Fig.2 as a beam of 10 scaled distributions which are scattered now around asymptotics (4.6).

5. Diffusion amid virtual invariant curves: the Lyapunov exponents

Above we considered a very particular and most simple limiting case of our model (2.2) with parameter $`d=0`$. In this case the motion is ergodic  which greatly simplifies the problem under consideration. Nevertheless, we obtained a number of new results which form the firm foundations for further studies.

The most important new feature of the motion for $`d>0`$ is the so–called divided phase space of the system that is a mixture of both chaotic as well as regular components of the motion. This is a typical structure of a few–freedom dynamical system (see, e.g., ).

First of all, we need to exclude the regular trajectories from further analysis of the diffusion statistics. The standard well known method to do this is simultaneous computation for each trajectory of the so–called Lyapunov exponent $`\mathrm{\Lambda }`$ that is the rate of the local exponential instability of the motion (see, e.g.,  and references therein). In a 2D canonical (hamiltonian) map like our model (2.2) there are two Lyapunov exponents whose sum is always zero: $`\mathrm{\Lambda }_1+\mathrm{\Lambda }_2=0`$. For a chaotic trajectory one exponent, say, $`\mathrm{\Lambda }_1=\mathrm{\Lambda }_+>0`$ is positive while another one $`\mathrm{\Lambda }_2=\mathrm{\Lambda }_{}<0`$ is negative. As a result, according to the standard definition of the Lyapunov exponent in the limit $`t\mathrm{}`$, any tangent vector $`(dx,dp)`$ of the linearized motion approaches the eigenvector corresponding to $`\mathrm{\Lambda }_+>0`$.

A simple well known procedure for computing $`\mathrm{\Lambda }_+`$, we made use in the present work too, is the following. For each of $`M`$ trajectories with random initial conditions $`x_0,p_0`$ we chose the tangent vector $`(dx,dp)`$ of random direction and unit modulus: $`d\rho ^2=dx^2+dp^2=1`$. Then both maps, the main one and the second one linearized with respect to the main reference trajectory $`x(t,x_0,p_0),p(t,x_0,p_0)`$ were run simultaneously during some time $`t`$. Finally, the current $`\mathrm{\Lambda }(t)`$ was calculated from the standard relation:

$$\mathrm{\Lambda }(t)=\frac{<\mathrm{ln}\rho (t)>}{t}$$
$`(5.1)`$
where the brackets denoted the averaging over $`M`$ trajectories. Unlike the formal mathematical definition of $`\mathrm{\Lambda }`$ in the limit $`t\mathrm{}`$, in numerical experiments the Lyapunov exponent $`\mathrm{\Lambda }(t)`$ is always time dependent, perforce.

In Fig.3 a few typical examples of the $`\mathrm{\Lambda }`$ distribution are depicted for the number of events in (4.1) $`J=M`$ equal to that of trajectories while the number of printed–out points $`J_0=M^{}M`$ is less except the case $`d=0`$. The simplest one is for ergodic motion ($`d=0`$). It has a form of almost vertical step which derivative $`dF/d\mathrm{\Lambda }10^4`$ is very narrow $`\delta `$-function. Notice that regular chain of points along $`F`$ axis has no special physical meaning but simply reflects a particular accepted type of the distribution $`F(\mathrm{\Lambda }_j)=j/J`$ (4.1) with integer $`j`$. The mean $`\mathrm{\Lambda }`$ depends only on $`K`$ (see Eq.(2.10)) but not on the initial conditions. This example in Fig.3 shows the empirical/theoretical ratio which is very close to unity as expected.

Two other examples correspond to the same $`K=0.45`$ and $`M=10^4`$ but different motion time $`t=10^4`$ and $`10^5`$ iterations. Both distributions have the same step at the largest $`\mathrm{\Lambda }`$ which corresponds to diffusive components (not necessarily a single one) of the motion similar to the ergodic case. However, the most interesting part of the former is the rest of the distribution which represents a rich motion structure contrary to a dull one in the ergodic motion.

The largest (but again not the most interesting) part of this structure is related to the distribution steep cut-off at small $`\mathrm{\Lambda }`$. Comparison of the two distributions for different motion time $`t=10^4`$ and $`10^5`$ shows that the $`\mathrm{\Lambda }`$ values of the trajectories in this region decrease with increasing time, approximately as $`\mathrm{\Lambda }1/t`$. This would mean that all these trajectories are regular (see Eq.(5.1)) because the tangent vector $`\rho `$ does not grow. The relative number of such trajectories gives the total area of regular motion on system’s phase cylinder. In a particular example under consideration it amounts to $`A_{reg}=3177/100000.318(t=10^5)`$. Generally, that value depends on a particular cut–off border from above chosen (see arrow in Fig.3). This delicate experimental problem is considerably mitigated by a lucky feature of $`\mathrm{\Lambda }`$ distribution in our model, namely, a relatively wide plateau of $`F(\mathrm{\Lambda })`$ immediately above the cut–off with only a few trajectories on it. However, the statistical accuracy

$$\frac{\mathrm{\Delta }A_{reg}}{A_{reg}}(MA_{reg})^{1/2}$$
$`(5.2)`$
is typically much worse, and can be improved by increasing the number of trajectories (and the computation time) only.

Another interesting feature of $`\mathrm{\Lambda }`$ distribution in our model is a characteristic ”fork” shape of the cut–off. This is a result of negative $`\mathrm{\Lambda }`$ for many regular trajectories. Such a peculiar representation is obtained by ordering $`\mathrm{\Lambda }(t)`$ values with their signs but plotting out the moduli $`|\mathrm{\Lambda }(t)|`$ only. Thus, the lower prong of the fork corresponds to negative $`\mathrm{\Lambda }(t)<0`$ while on the upper one $`\mathrm{\Lambda }(t)>0`$ are positive. This is due to the complex–conjugate Lyapunov exponents which results in a strictly bounded oscillation of the tangent vector $`(dx,dp)`$ in this case. However, the area (5.4) $`A_\pm 0.20<A_{reg}0.318`$ is noticeably less than the total regular domains $`A_{reg}`$. The rest is filled with the trajectories which are also regular but linearly unstable. This means the linear in time growth of the tangent vector $`\rho (t)t`$ so that $`\mathrm{\Lambda }(t)0`$ remains positive but is vanishing in the limit $`t\mathrm{}`$. This is the so–called marginal local instability with both $`\mathrm{\Lambda }_\pm =0`$ equal zero (for discussion see ). A curious point is that this seemingly exceptional case becomes the typical one in a nonlinear oscillator system due to the dependence of oscillation frequencies on the trajectory initial conditions. In fact, the bounded $`\rho `$ oscillation producing negative $`\mathrm{\Lambda }(t)`$ is the exceptional case instead. The origin of this peculiarity is in a piecewise linear force in our model (2.2). As a result, the motion in the main (and the biggest for large $`K`$) regular domain around fixed point $`x=1/2,p=0`$ is plainly a harmonic oscillation with the frequency ($`K<d`$)

$$\mathrm{\Omega }=\mathrm{arccos}\left(1\frac{K}{d}\right)1.47$$
$`(5.3)`$
which remains the same in the whole regular domain of area

$$A_\pm =\frac{2\pi K}{d}y_\pm ^2\left(1\frac{K}{2d}\right)0.20$$
$`(5.4)`$
Here $`y_\pm =x_\pm 0.5=\pm d/2`$ is the position of two singularities of the force (see Eq.(2.2) and below) which restrict the size of the regular domain surrounded by the limiting ellipse to which both lines of singularity $`y_\pm =d/2=0.25`$ are tangent. This ellipse is determined by the initial conditions

$$p_0=\mathrm{\hspace{0.17em}0},x_0=\mathrm{\hspace{0.17em}0.5}+y_\pm \left(1\frac{K}{2d}\right)\mathrm{\hspace{0.17em}0.5}\pm \mathrm{\hspace{0.17em}0.185}$$
$`(5.5)`$
All the numerical values above correspond to $`K=0.45`$ and $`d=1/2`$. Within the ellipse the motion of tangent vector obeys the same equation as the main motion, the only difference being an arbitrary length $`\rho `$ of the tangent vector (for details see  and references therein).

Coming back to Fig.3 notice the decrease of the measured area $`A_\pm `$ with increasing motion time. This is explained by penetration of trajectories into a very complicated critical structure at the chaos border surrounding each regular domain (for details see, e.g., ). For the same reason, the direct measurement of the whole regular region $`A_{reg}0.40`$ by a single chaotic trajectory as long as $`10^9`$ iterations gives a noticeably larger value as compared to $`A_{reg}0.318`$ from $`10^4`$ trajectories for $`10^5`$ iterations each.

In all the curiosity of the $`\mathrm{\Lambda }(t)`$ distribution in regular components of the motion our main interest in the present study is the intermediate region between the regular cut–off at most small $`\mathrm{\Lambda }(t)0`$ and the chaotic step at maximal $`\mathrm{\Lambda }`$ independent of $`t`$. The distribution in this region also does not depend on motion time and characterizes the proper critical structure of the chaotic motion. In an example in Fig.3 this structure is presented by a relatively small probability step $`\mathrm{\Delta }F0.06`$ at $`\mathrm{\Lambda }0.03`$. In the next Section a few other examples will be given too.

6. Diffusion amid virtual invariant curves: the Critical statistics

In Fig.4 we present three characteristic examples of the effect of the critical structure on the diffusion statistics. The dashed curve shows the ”unperturbed” distribution $`F(D^{})`$ (4.1) of the normalized diffusion rate $`D^{}=D/D_{norm}`$ where the normalizing rate $`D_{norm}`$ to be chosen in each particular case (see below). The term unperturbed means here the ergodic case $`d=0`$ without any invariant curves and critical structure (Section 4, the problem of critical structure in this case is not as simple as it may seem, see below and Section 7). In this case the normalizing rate $`D_{norm}=D_{\mathrm{}}`$ is the true asymptotic diffusion rate (4.5).

Now we are interested in the effect of the critical structure which typically arises in a nonergodic motion with its barriers for the chaos, or chaos borders. The latter are particular, and very important, case of an invariant curve which is transformed into itself under system’s dynamics. As was discussed already above (Section 1) there are several different types of invariant curves (ICs).

One is well studied and rather familiar chaos border surrounding any domain with regular motion. In this paper we call it the local invariant curve (LIC) which does not block the global diffusion (GD) around such a domain. An important property of LIC is the robustness which means that a small change of the system, say, of a parameter $`K`$ or $`d`$ cannot distroy the LIC but may only deform it slightly. This implies that LICs are always present in any divided phase space.

Here we are mainly interested in a different IC type, the global invariant curves (GICs). Each GIC cuts the whole phase–space cylinder ($`xmod1`$) of our model, and thus completely prevents GD in $`p`$. Such ICs are less known, especially the most surprising of them, the separatrix of a nonlinear resonance. However, those GICss are not robust in the model under consideration (see ), being destroyed by almost any arbitrarily small perturbation of the system, particularly by a change of even a single its parameter. In other words, such GICs do exist for the special values, say, $`K=K_m`$ only. Even though there are typically infinitely many such special values the probability to find a GIC in a randomly chosen system is zero. This is why we are interested in a more generic situation when our model has no GICs at all. Yet, the effect of those still persists in a certain domain around each $`K_m`$ ! For this reason we call such GICs the virtual invariant curves (VICs) in analogy with other virtual quantities in physics like, for example, virtual energy levels in qantum mechanics. Notice that unlike a GIC the VIC is robust and, hence, generic.

Both LICs and GICs produce the so–called critical structure of the motion (see, e.g., ) which is typically characterized by a power–law distribution of principal quantities. The corresponding exponents $`c_n`$ are called the critical exponents. Their values are shown in Fig.4 at the related distributions. Notice that the opposite is generally not true that is a particular power law does not necessarily indicates any critical structure. In our model this is just the case for the ergodic motion where the diffusion rate distribution is also characterized by an asymptotic ($`D0`$) power law with exponent $`c_0=0.5`$ (see above and Section 7). However, an important difference between ergodic and nonergodic dynamics is in that all critical exponents in the latter case $`c_n<c_0`$ are less than (generally noncritical) ergodic exponent $`c_0`$. This is the main physical result of our preliminary numerical experiments we can present and discuss already right now (see Fig.4).

Let us start with the distribution for $`K=0.45`$ (upper solid line) which is far in the region without VICs (the border of this region is at $`K_B(d=1/2)=1/3`$, see Eq.(2.3) above and ). However, the regular trajectories ($`A_{reg}0.318`$) together with LICs and the related critical structure are present. As a result the distribution (with $`D_{norm}=D_{\mathrm{}}`$) is well deviates from the unperturbed one for ergodic motion with $`d=0`$. This type of the critical structure in a relatively narrow layer around a LIC is well studied by now (see, e.g., ) including some deviation of a typical distribution from a pure power law. The latter would mean the exact scale invariance of the underlaying critical structure in both the system phase space as well as its motion time.

The critical structure is described by the so–called renormalization group, or renormgroup for brevity. On the other hand, the motion equations for any dynamical system also form a certain (dynamical) group. Such a fundamental similarity allows one to interpret the critical structure as a certain dynamics which was called the renormdynamics . In this picture the exact scale invariance with a pure power–law distribution corresponds to the simplest, periodic, renormdynamics even though the original dynamics may be the most complicated chaotic motion. The resolution of this apparent paradox is in that the latter complexity of the original dynamics is ”transferred” to the dynamical infinitely dimensional space of the renormdynamics leaving behind the most simple renormdynamics itself (sometimes !).

The latter limit is most studied simply because it is the simplest one. However, the generic case is just opposite that is a typical renormchaos is also chaotic . Particularly, this implies a certain chaotic oscillation of the characteristic distribution around some average power law. This is just the case for the upper distribution in Fig.4 under discussion. It is characterized by the average critical exponent $`c_1=0.3`$ with fluctuations $`(c_1^{}c_1)=0.1`$. Such an interpretation of the critical structure in question is known to be typical but not necessarily unique (see below). The really unique property of this critical structure is the infinite power law, with or without the fluctuations. The term ”infinite” corresponds here to the range of a renormdynamical variable $`\mathrm{ln}D\mathrm{}`$ with unrestricted variation even though the diffusion rate itself $`D>0`$ is strictly bounded from below.

This is no longer the case for a new type of critical structure which we have encountered in our problem and which is produced by VICs (=robust GICs) rather than by the robust LICs. As was already explained above the principal difference between the two is in that the VIC is not an invariant curve at all. In terms of renormdynamics it means that VIC can mimic a GIC for relatively large $`\mathrm{ln}D`$ only. This is clearly seen in Fig.4 in the upper part of the distribution with local critical exponent $`c_2=0.09`$ and parameter $`K=0.335`$ (points). Here we have taken $`D_{norm}=10^6<D_{\mathrm{}}2\times 10^5`$ much smaller than the true diffusion rate $`D_{\mathrm{}}`$. This shifts the whole distribution to the right to avoid the overlapping with other distributions. This value is slightly above the border $`K_B(1/2)=1/3`$ (see Eq.(2.3)) where there is a lot of VICs without any GIC. As a result the range of characteristic critical exponent $`c_2`$ is very short: $`\mathrm{\Delta }\mathrm{ln}D^{}5`$ as compared with the total available range $`25`$. The rest of the distribution remains fairly close to the unperturbed one. This would mean the absence of the critical structute over there or, at least, its sharp change at $`\mathrm{ln}D\stackrel{<}{}2`$. In case of the former interpretation the renorm-motion stops in the latter region.
In turn, this would imply a ”dissipative” rather than ”Hamiltonian” renormdynamics. Notice that the main part of the distribution is close but not identical to the unperturbed one including a slight difference in the characteristic exponent. Does it means a certain very slow renorm-motion remains a very interesting open question yet to be studied. Interestingly, the larger critical exponent $`c_2^{}=0.45`$ is also close to the local critical exponent $`c_1^{}=0.4`$ in the region without VICs or GICs and was interpreted above as a random fluctuation in the renormchaos. Is it really true still remains unclear.

Finally, the third distribution in Fig.4 (lower solid line) is actually coincides with the unperturbed distribution ($`D_{norm}D_{\mathrm{}}`$) even though it corresponds to the region with many VICs and strong suppression of the diffusion ($`K=0.3294`$, see Fig.3 in ). A deviation for very small $`D^{}`$ is due to a poor statistics at this end. Notice that the coincidence of both distributions is not only asympotic ($`F0`$) but complete, including the opposite limit $`F1`$. And this is in spite of a rather large regular region $`A_{reg}0.581`$. The origin of this peculiarity for a particular $`K`$ value remains unclear. One possibility is that, for some reasons, the area of the critical structure at the chaos border around this regular domain is unusually small. Examples of such a peculiarity in different models are known (see ). Actually, in the latter work the critical structure was found to be unusually large but hidden. In other words, the motion was ergodic but with strong correlations (cf. the unusual diffusion rate (2.8) for $`K1`$ in ergodic system at $`d=0`$). Coming back to this case in Fig.4 we can conclude that our ”unperturbed” power–law distribution with exponent $`c_0=1/2`$ (dashed line) may well represent a peculiar critical structure related to the strong hidden temporal correlations rather than to a purely spatial geometry of the phase space. If this is true, indeed, the correlation decay may well be not a power–law one at all as is the case in model  where such a hidden decay is purely exponential (see Fig.6 over there).

At last, let us mention another peculiarity of the critical structure in question: all the critical exponents found so far are less, if only a little, than ”unperturbed” or ”hidden” one $`c_0=1/2`$. The physical meaning of this universal inequality is in that the critical structure under consideration does always increase the probability of very low diffusion rate $`D0`$. The general mechanism of this is known (see, e.g., ), and is explained by the ”sticking” of a trajectory within a complicated critical structure which slows down the diffusion. Interestingly, that the sign of the sticking effect may be opposite when the sticking accelerates the diffusion up to the absolute maximum $`D(t)t`$ of homogeneous diffusion rate .

To summarize, we see that our ”simple” model in the present paper really reveals a great variety of critical structure still to be farther studied and understood.

7. Conclusion: A hidden critical structure ?

In this paper we present some preliminary results of the numerical experiments with a family of simple models specified by the smooth canonical 2D-map (2.1) with global virtual invariant curves (VICs). As in our previous paper  we make use here of the same strongly chaotic model, and focus again on the statistics of the diffusion rate $`D`$ which proves to be of a very complicated (apparently fractal) type determined by the so–called critical structure of both the phase space and the motion (see, e.g., ). In previous paper  we studied the statistics of the mean diffusion rate $`<D(K)>`$ averaged over ensemble of trajectories with random initial conditions. Our main result there was observation of very big and irregular fluctuations of the dependence $`<D(K)>`$, and a long and very slowly decaying tail of $`<D>`$ distribution for $`<D>0`$. We termed the latter effect the VIC diffusion suppression.

In the present paper we continue studying this interesting phenomenon in more details. To this end we turn from the statistics of averages $`<D(K)>`$ as a function of model parameter $`K`$ to that of individual trajectories for a given $`K`$. In principle, such an approach provides the most deep insight into a statistical problem. As the main statistical characteristic we have chosen the integral distribution $`F(D)`$ in the form (4.1) for a poor statistics at $`D0`$. Preliminary results of our extensive numerical experiments presented in Fig.4 confirm, indeed, our earlier conjecture on a critical structure underlaying the fractal dependence $`<D(K)>`$ in , the true sign of such a structure being various power–law distributions found. Moreover, besides the familiar well known critical structure exemplified in Fig.4 by the case with parameter $`K=0.45`$ we did observe many cases of a rather different structure as one with $`K=0.335`$. The principal difference of the latter is its finite size in the structure variable $`\mathrm{\Delta }\mathrm{ln}D\stackrel{<}{}5`$. The natural explanation of this difference is the following. First of all, the VIC is not a true invariant curve like GIC. The latter completely blocks the GD while the former can, at most, inhibit the diffusion only. The mechanism of inhibition is known to be the sticking of a trajectory inside a very complicated critical structure. In turn, the sticking is the stronger (longer) the smaller is the spatial and/or the longer is the temporal scale of the critical structure. But for the VIC structure both are strictly restricted. On the other hand, such a restriction is the weaker the higher is the VIC density. In the system under consideration the latter is rather large, and hence the restriction leaves enough freedom for a strong suppression of GD for almost any $`K`$. Moreover, since the critical exponent of the VIC structure is typically very small (for example $`c_2=0.09`$ in Fig.4) the probability of large suppression is high even for a short critical structure (cf.  for a different characteristic of this phenomenon). This slow decaying suppression probability is well ascertained in our numerics but we have not, as yet, any theoretical explanation of such a behavior.

Now we come, perhaps, to the most interesting result of our current studies. How strange it may seem, this brings us to the apparently simplest case of our model with $`d=0`$ when the motion is ergodic. Can it still reveal any structure on the grounds that the distribution $`F(D)`$ is also a power law (Fig.4)? That is the question! Certainly, it is not the case if, in addition, parameter $`K1`$ is big, and the diffusion rate has the standard form (2.11): $`DK^2`$. But what if $`K1`$ is small? At least, the diffusion rate becomes qualitatively different: $`DK^{5/2}`$. What does it mean? Generally, nothing! But in a particular case under consideration such dependence $`D(K)`$ can be, and actually was, derived  from the resonance structure of the motion. If the system were not ergodic (with divided phase space) this structure would be clearly seen in the phase space. The question is what happens for ergodic motion with the same dependence $`D(K)`$? In  we conjectured that some structure would still persist in the form of the correlations which determine the diffusion rate that is in some ”hidden” form unseen directly in the picture of the motion in phase space. An example of such hidden critical structure was found in  (see Section 6 above). However, in that case a particular distribution function was exponential rather than a power–law one (?). Hence, the question is if this qualitative difference could depend on a particular characteristic of the critical structure? Still another question arises from a very strange temporal behavior of the diffusion rate in the same ”simple” case of ergodic motion for $`d=0`$. We mean a ”mysterious” plateau at the very beginning of diffusion under a weak perturbation ($`K1`$, see Fig.1). In this case the dependence $`D(K)=K^2/3`$ is the same as in the opposite limit of strong ($`K1`$) uncorrelated perturbation (?) but for a very short time only, the shorter the stronger is the perturbation (?!). Moreover, the correlations on the plateau are not only very large as in the weak–perturbation limit $`K0`$ but even increasing during the whole plateau regime (see Fig.1, dashed lines for variances $`V_M(\tau )`$ (3.2)). As yet, we have no definite explanation for such a controversial behavior. A discreet current conjecture is the following. The duration of plateau is $`\tau _{pl}1`$, or $`t_{pl}1/\mathrm{\Lambda }1/\mathrm{\Omega }`$ (see Eq.(2.9)). But the latter expression gives the period of phase oscillation on the critical nonlinear resonance which determines the diffusion rate . Then, one can image that this period does characterize not only the correlation decay, as usual, but also the correlation uprise. Still, the invariable diffusion rate over the whole plateau region is to be explained yet.

In any event, we are very curious to continue these most interesting studies!

Acknowledgements. This work was partly supported by the Russia Foundation for Fundamental Research, grant 01-02-16836. We are grateful to Ms. L.F. Hailo for her permanent and very important assistance in computer experiments.

Figure captions

* Diffusion relaxation $`D^{}(\tau )=D(\tau )/D_{\mathrm{}}1`$ in model (2.2) with parameter $`d=0`$ (without invariant curves) is presented as a function of dimensionless time $`\tau `$ (3.5) for two values of $`K=0.01`$ (circles) and $`3\times 10^5`$ (crosses). Two smooth solid lines show the empirical relation (3.5) with two fitting parameters $`c=1`$ and $`\gamma =4`$. Dashed lines are variances $`V_M(\tau )`$ (3.2), and dotted lines show variances $`V_N(\tau )`$ (3.8) (see text). In the lower part the scaling (3.6) is presented reduced by factor 10 to avoid overlapping with other data. The full volume of empirical data was $`J=M\times N=10^4\times 10=10^5`$ (see text).
* Distribution function $`F(D^{})`$ (4.1) of the reduced diffusion rate $`D^{}`$ (4.5) in model (2.2) without invariant curves ($`d=0`$). Thick dashed straight line represents asymptotics (4.6) of the integrated $`D`$-distribution (4.2) for the Gaussian $`p`$-statistics. Two lower wiggly lines correspond to large deviations from the latter: $`D_{}=42(K=10^3)`$ and $`461(K=3\times 10^5)`$ (see Insert). A group of 10 $`D`$-distributions in a large interval ($`10D_{}461`$) are brought together using empirical relation (4.9). Insert: shift factor $`R_D`$ (see text) vs. deviation $`D_{}`$ (4.8) for $`K=10^3`$ (circles) and $`3\times 10^5`$ (crosses); the straight line is empirical relation (4.7).
* Examples of distribution function $`F(\mathrm{\Lambda })`$ of type (4.1) but for the Lyapunov exponent in model (2.2) for $`d=0,M=M^{}=80,t=10^4`$ (the right-most step $`F(\mathrm{\Lambda })`$, ergodic motion), and for $`d=1/2,M=10^4,M^{}=1000,t=10^4,10^5`$ (nonergodic motion, see text); in all cases $`K=0.45`$. Horizontal line indicates the total share $`A_{reg}0.318`$ of the motion regular components. Arrow at $`\mathrm{\Lambda }=10^4`$ shows the lower border of chaotic trajectories chosen for further analysis (for $`t=10^5`$, see text).
* Three characteristic examples of the diffusion statistics in the critical structure including virtual invariant curves ($`d=1/2`$). Shown are integral distributions $`F(D^{})`$ (4.1) of the normalized diffusion rate $`D^{}=D/D_{norm}`$ (see text). The numbers at curves are critical diffusion exponents $`c_m`$. The biggest one $`c_0=1/2`$ corresponds to the ergodic motion ($`d=0`$) without any critical structure (dashed curve). Two straight lines show the average ($`c_1=0.3`$) and local ($`c_1^{}=0.4`$) critical exponents for $`K=0.45`$ (solid line connecting 500 values of $`F(D^{})`$). The distribution for $`K=0.335`$ with two local critical exponents ($`c_2=0.09`$ and $`c_2^{}=0.45`$) is presented by 300 points shifted to the right to avoid overlapping with two other distributions. The third distribution (a solid line through 1000 points, $`K=0.3294`$) is surprisingly close to that in ergodic case (dashed line). In all three examples $`M=10^4,t=10^5`$.


format: text

text:
# Global Optimization by Adiabatic Switching

## I Introduction

Determination of the lowest energy configuration for a cluster of $`N`$ atoms is a nontrivial task . The complexity arises in part from the exponentially (in $`N`$) large number of minima in the potential energy surface . Furthermore, the geometry of the potential energy landscape itself can make it computationally hard.

The problem is, however, simple to describe: Find the lowest energy minimum of a $`N`$ body potential energy surface, $`𝒱(\stackrel{}{𝐫})`$,

$$𝒱(\stackrel{}{𝐫})=\underset{i,j}{}V(r_ij)$$
(1)
where $`\stackrel{}{𝐫}`$ are the atomic coordinates, and $`V(r_{ij})`$ is the interatomic potential of interaction between atoms $`i`$ and $`j`$. For small $`N`$ one can hope to enumerate all possible minima and decide the lowest of these, but even for moderate $`N`$ and for the simplest $`V`$ such as the Lennard–Jones (LJ) potential, this becomes difficult. A case in point is the 38 atom LJ cluster which has the so–called “double funnel” structure; the global minimum, which has octahedral symmetry, is marginally lower than the first excited state which has icosahedral symmetry. These were respectively found by the basin hopping technique  and a genetic algorithm method .

A number of techniques of global optimization have been applied to this problem  and by now there are extensive compilations of global minima for a number of different clusters  , notably those described by two–body or many–body potentials which are commonly applied in atomistic simulations. A major difficulty is in ensuring that the algorithms reach the global minimum without being trapped in local minima. One method of overcoming such trapping  is by transforming the PES, broadening the thermodynamic transitions so as to increase the probability of finding the global minimum at temperatures where the free energy barriers are almost unsurmountable. For example, addition of a linear term to the PES provides a compressing effect which has been shown to be successful in locating the true minima in multiple funneled global structures . Locatelli and Schoen  used such transformations to locate the global minimum for non-icosahedral clusters.

We have recently proposed a new method of global optimization wherein time–dependence is introduced in the potential energy landscape . The evolving landscape is designed in a manner such that asymptotically, the potential energy surface develops into the surface of interest. A number of other techniques can be used to follow the evolving minima. Applications have been made to determining the ground state configurations of simple cluster systems .

In the present paper we apply this method to determine the ground state configurations and energies of atomic clusters described by the many–body Sutton–Chen potential  by switching from a known ground state. Initially, the interaction is chosen to be the Lennard–Jones potential,

$$𝒱_0(\stackrel{}{𝐫})=\underset{ij}{}V_{LJ}(r_{ij})=\underset{ij}{}4ϵ[(\sigma /r_{ij})^{12}(\sigma /r_{ij})^6]$$
(2)
while the surface of interest is the potential

$$𝒱_f(\stackrel{}{𝐫})=ϵ\underset{i}{}[\frac{1}{2}\underset{ji}{}(\frac{a}{r_{ij}})^nc\sqrt{\rho }_i],\rho _i=\underset{ji}{}(\frac{a}{r_{ij}})^m$$
(3)
One choice for the time dependent potential energy surface is

$$𝒱(t)=𝒱_0(\stackrel{}{𝐫})g(t)+𝒱_f(\stackrel{}{𝐫})h(t)$$
(4)
with $`g(t)`$ an adiabatically varying switching function that interpolates between 1 and 0, and $`h(t)`$ doing the reverse.

In the next section we describe the method as applied to the problem of ground state energy determination for Sutton–Chen clusters. Detailed results are presented for one cluster size, while the more general application and results are indicated in brief. This is followed by a discussion and summary.

It is a pleasure to dedicate this article to Steve Berry who has directly and indirectly influenced much of the development in the area of cluster studies over the past few decades. We have learned a lot from him, both in conversation as well as through his many articles and reviews .

## II Adiabatic optimization

The adiabatic optimization method , is a heuristic technique for locating minima. The essential idea is as follows.

Time dependence is introduced into the potential energy landscape directly by the incorporation of slowly varying terms as discussed in Eq. (4). A given choice is made for the switching functions $`g(t)`$ and $`h(t)`$, though in practice the choice does not affect the results greatly. A similar application of the adiabatic principle to determine semiclassical ground states of multidimensional systems  has noted the insensitivity of the technique to the precise form of the switching function, so long as the induced variation of the potential energy surface is slow enough. We note, parenthetically, that the switching principle has wide applicability, and in recent work has been used in the computation of the free–energy of finite clusters .

Location of the evolving minima can be done by any of a number of techniques. The simplest procedure is to introduce damping into the equations of motion and allow the system to evolve to a position of rest in a potential minimum; by starting with an ensemble of initial configurations and varying the available parameters, a number of minima can be located, and the putative global minimum can be recognized. Elsewhere  we have suggested the conjugate gradient  or simulated annealing (SA)  as other possible methods for locating the minima. It is likely that of these, the conjugate gradient technique will be more efficient as compared to SA though some SA variants  may also provide a suitable method for following the evolving minima.

The overall procedure can be summarized as follows:

1. Take the initial configuration of the $`N`$ atom cluster to be the ground state for the LJN cluster .
2. Choose some switching function, say $`g(t)`$. Similarly choose $`h(t)`$, and the simplest choice, which we make here, is $`h(t)=1g(t)`$. We have explored a large variety of switching functions and in the present application we use $`g(t)=\mathrm{cos}^2(3\pi \zeta t)\mathrm{exp}(\zeta t)`$, where $`\zeta `$ is the adiabaticity parameter.
3. Perform molecular dynamics simulations for this $`N`$–particle cluster with forces deriving from Eq. (4), with an additional damping term, namely the equations of motion

$$m\ddot{\stackrel{}{𝐫}_k}+\gamma \dot{\stackrel{}{𝐫}_k}+\frac{𝒱}{\stackrel{}{𝐫}_k}=0,k=1,2\mathrm{},N,$$
(5)
where $`\stackrel{}{𝐫}_k`$ is the position vector for the $`k`$th particle, $`m`$ is its mass and $`\gamma `$ is the damping constant.
4. Vary $`\zeta `$ and $`\gamma `$, keeping in mind the natural timescales of the problem. Evolve to a minimum energy configuration, namely when the particle velocities become zero; the lowest energy found in an ensemble of simulations is the ground state energy predicted by the present method.

### A Results for Sutton–Chen global minima

Here we attempt to switch from the minimum of the LJN system to the minimum of the SCN system. Both sets of minima have been extensively studied earlier and are tabulated in the Cambridge Cluster Database . A point of interest is that for the Sutton-Chen 9-6 family of potentials, the symmetries of the global minimum configurations are frequently different from the symmetries of the Lennard–Jones minima, so that in the adiabatic switching process, the cluster atoms must also move so as to adopt a different symmetry.

We present detailed results for the cluster size $`N=15`$, though we have applied this technique to larger clusters and obtained results in agreement with the current standards . For the 15 atom LJ cluster, the ground state has the point group symmetry C2v while for the SC cluster the symmetry is D6d. Shown in Fig. 1 is a plot of the potential energy versus time for a particular choice of $`g(t),\zeta `$ and $`\gamma `$. Also shown is the effect of instantaneous switching, namely taking the limit $`g(t)=0`$, where it can be seen that the system finds the nearest available local minimum from which it does not move. The time–dependence in the potential effectively permits the system to explore the multidimensional potential energy landscape of the SC cluster in an efficient manner. Finding a local minimum does not trap the system since there is always kinetic energy until the adiabatic switch is essentially over. Inset in the figure is a schematic of the cluster configuration at different times during the process, showing how the cluster both contracts as well as rearranges to eventually reach the minimum of the SC surface.

It should be added that we have performed simulations for a variety of cluster sizes and in all cases we find that the procedure successfully finds the tabulated minima of SC clusters; these are not presented here since the details are repetitive. As we have emphasized elsewhere , the present method is heuristic, and thus some exploration of different switching functions, variation in the adiabaticity and damping parameters, and indeed the choice of initial potential, $`𝒱_0(\stackrel{}{𝐫})`$ is necessary.

## III Summary and Discussion

In this paper we have presented the outline of a general procedure for global optimization with specific application to the problem of cluster ground state geometry determination. The application here, to the determination of the minimum of model metallic (Sutton-Chen) clusters by adiabatically deforming the potential energy surface relevant to model rare–gas (Lennard Jones) clusters is meant to be illustrative rather than exhaustive: the method introduced here is one of a class of techniques that employs time–dependence in the potential energy surface to enhance the exploration of phase space in contrast to other means of achieving the same objective .

A multiplicity of techniques is needed to approach hard problems such as global optimization. Few rigorous results are available, and application of most techniques is not guaranteed, with few possible exceptions, to give reliable (or certifiable) results. The present adiabatic switching method locally solves the optimization for an evolving surface, and thus mimics other methods of making large scale excursions in configuration or phase space.

We are presently studying this technique in detail with respect to the variation of parameters as well as to functional variations. One of the main issues of concern, and one that we are addressing in current work, is the relative efficiency of this method in comparison to other global optimization techniques. In a number of applications, we find that this method gives very encouraging results, and permits the determination of fairly reliable minimum energy configurations for a wide variety of cluster systems. The flexibility of choice of a number of starting potentials including the free particle case  as well as the flexible choice of switching functions and parameters, and finally the flexibility in the dynamical evolution all combine to suggest that while the method is heuristic, it holds promise.

## Acknowledgment:

This work is supported by a grant from the Department of Science and Technology. We thank Subir Sarkar for discussions.


text:
# Soft X-ray emission lines from photoionized accretion discs: constraints on their strength and width

## 1 Introduction

Features observed in the hard X-ray spectra of active galactic nuclei (AGN) are one of the most direct probes of the structure of accretion discs close to the central black hole. The prominent iron K$`\alpha `$ emission line at 6.4 keV and Compton hump at $``$20–30 keV in the spectra of many Seyfert 1 galaxies indicates that optically thick material is being illuminated by the hard X-ray source (Pounds et al., 1990). As pointed out by Fabian et al. (1989), if this reflected emission arose from the accretion disc, then relativistic and Doppler effects would broaden the emission lines in a characteristic way. Observations of such a line allows a determination of not only the accretion disc structure close to the black hole, but also the space-time geometry in the strong gravity regime. The ASCA observatory (1993–2001) finally provided the sensitivity to search for relativistically broadened Fe K$`\alpha `$ lines, and detected a clear example in the spectrum of the bright Seyfert 1 galaxy MCG–6-30-15 (Tanaka et al., 1995). Since that initial discovery, broad Fe K$`\alpha `$ lines have been found in many other Seyferts (Nandra et al., 1997; Yaqoob et al., 2002), and are now an integral part of AGN phenomenology (Fabian et al., 2000).

Calculations of the reprocessed emission from an irradiated slab of gas have shown that the Fe K$`\alpha `$ line, by virtue of its large fluorescent yield and relatively high cosmic abundance, is the strongest line in the reflection spectrum (George & Fabian, 1991; Matt, Perola & Piro, 1991). Yet, strong lines from abundant elements such as carbon, nitrogen and oxygen are also present at lower energies, and, if detected, would provide further information on the ionization state of the accretion disc. Models of reflection from photoionized discs have shown that the soft X-ray features are very sensitive to the ionization parameter of the disc material (Ross & Fabian, 1993; Ross, Fabian & Young, 1999). Determining the ionization state of the disc would constrain the magnitude of the illuminating radiation, and would therefore place limits on the amount of accretion energy that is released in the corona.

Although the soft X-ray spectra of Seyfert galaxies are complicated by features such as warm absorbers and soft excesses, the sensitivity and resolution provided by the instruments on board XMM-Newton and Chandra could allow relativistic soft X-ray emission lines to be found amongst the clutter and noise. Indeed, claims of such lines from hydrogenic carbon, nitrogen and oxygen have already been made based on XMM-Newton Reflection Grating Spectrometer (RGS) spectra of MCG–6-30-15 and Mrk 766 (Branduardi-Raymont et al., 2001; Sako et al., 2002). Clearly, the time is right for a closer look at the properties of the soft X-ray accretion disc lines.

Here, we employ detailed models of photoionized accretion discs to investigate the properties (specifically, the strength and width) of the C, N & O lines found in the reflection spectra. A range of ionization parameters are considered, and both constant density slabs and hydrostatic atmospheres are used for the disc structure. The calculations and assumptions are outlined in the following section, with the results presented in Section 3. Finally, we discuss the implications of our results in light of the recent XMM-Newton observations and draw conclusions in Section 4.

## 2 Computations

Models of illuminated accretion discs have greatly increased in sophistication over the last couple of years. Initially, the calculations had assumed that the surface of the disc was a constant density slab (Ross & Fabian, 1993; Życki et al., 1994; Magdziarz & Zdziarski, 1995; Ross et al., 1999), and so the resulting gas structure and reflection spectrum could be characterized by the ionization parameter

$$\xi =\frac{4\pi F_\mathrm{X}}{n_\mathrm{H}},$$
(1)
where $`F_\mathrm{X}`$ is the X-ray flux incident on the slab, and $`n_\mathrm{H}`$ is the hydrogen number density of the material. While such a simple structure has some theoretical justification (a standard optically-thick disc has a roughly constant vertical density profile in the radiation pressure dominated limit; Shakura & Sunyaev 1973), the heated surface of the slab was out of pressure equilibrium with the remainder of the material. Recently, more realistic calculations were performed with the surface of the disc required to remain in hydrostatic balance (Nayakshin, Kazanas & Kallman, 2000; Ballantyne, Ross & Fabian, 2001). Under this condition, the irradiated gas takes on a more ‘two-phase’ structure with an outer hot ($``$ Compton temperature), ionized layer in pressure balance with a denser, colder and more neutral inner zone. The spectral features in the reflection spectra depend on the sharpness of the transition between the two layers (Nayakshin et al., 2000; Ballantyne & Ross, 2002), but qualitatively appear to be a diluted version of the constant density models.

In this section, we detail ionized disc calculations made with both the constant density code of Ross & Fabian (1993) and the hydrostatic code of Ballantyne et al. (2001). Since these latter models are dependent on the unknown accretion disc parameters, this paper concentrates on the results from the constant density models. However, an example hydrostatic calculation will be shown to illustrate the similarities and differences between the two approaches.

Important updates were made to both codes in order to better determine the soft X-ray emission spectrum. First, more elements were treated, and more recombined ions of each element were included. This was necessary so that nitrogen lines would be included in the output spectrum and so that the opacity would be treated accurately throughout the soft X-ray portion of the spectrum. In addition to fully-ionized species, the following ions are included in the calculations: C iii–vi, N iii–vii, O iii–viii, Ne iii–x, Mg iii–xii, Si iv–xiv, S iv–xvi and Fe vi–xxvi. Details of the atomic physics used will be presented in a subsequent paper. Second, the appropriate fraction of emission-line photons that escape (in narrow profiles) directly from the gas without undergoing any further continuum absorption or Compton scattering is kept track of separately from the rest of the diffuse radiation originating within the gas (Ross, 1979). Finally, the L shell lines of Fe xviixxii were changed from 3d$``$2p transitions to 3s$``$2p transitions. Like the other prominent lines, the Fe L-lines are produced by recombination. Liedahl et al. (1990) have pointed out that recombination cascades favour the production of 3s$``$2p lines over 3d$``$2p lines, and this effect is enhanced by resonance trapping of 3d$``$2p lines in optically thick gas. We assume that all recombinations to Fe xvii-xxii, except for radiative recombinations directly to the 2p level, ultimately result in 3s$``$2p emission.

The results of Ross & Fabian (1993) and Ross et al. (1999) have shown that, over a wide range of $`\xi `$, the Ly$`\alpha `$ line of O viii is typically the strongest soft X-ray emission line in the reflection spectrum. Therefore, this is the line that is most likely to be detected by current instruments and so we computed models where it is expected the be at its most prominent. Specifically, reflection spectra were calculated for $`\xi =`$250, 500 & 1000 erg cm s-1. The incident spectrum was assumed to be a power-law with photon-index $`\mathrm{\Gamma }=1.8`$, typical of a Seyfert 1 galaxy (e.g., Nandra et al., 1997). Furthermore, a 10 eV blackbody illuminated the base of the slab (at a Thomson depth $`\tau _\mathrm{T}=10`$) with six times the flux of the power-law. An example of one of the computed spectra is shown in Figure 1.

To investigate the effects of different metal abundances, a set of models was computed with five times the solar oxygen abundance, and another set with five times the solar nitrogen abundance. The initial set of abundances was taken from the list of Morrison & McCammon (1983).

The hydrostatic model requires many input parameters as the height of the disc at the irradiated point has to be explicitly calculated. For generality, we chose values typical of a Seyfert 1 galaxy: black hole mass $`M=10^7`$ M, accretion rate $`\dot{m}=0.01`$ (where $`\dot{m}\dot{M}/\dot{M}_{\mathrm{Edd}}`$ and $`\dot{M}_{\mathrm{Edd}}=48\pi Mm_\mathrm{p}G/\sigma _\mathrm{T}c`$, is the Eddington accretion rate), and the reflection occurs at a disc radius of 7 Schwarzschild radii. Radiation-pressure dominated boundary conditions were assumed. As with the constant density models, the surface of the disc was illuminated with a $`\mathrm{\Gamma }=1.8`$ power-law, but with an incident-flux-to-disc-flux ratio $`F_\mathrm{X}/F_{\mathrm{disc}}=1`$. These parameters resulted in a reflection spectrum similar to the $`\xi =500`$ erg cm s-1 constant density model (see Sect 3.2) and thus isolates the differences due to the different density structures111If $`F_\mathrm{X}/F_{\mathrm{disc}}=1/6`$ rather than 1 for the hydrostatic model, the added Compton cooling decreases the Compton temperature and hence the effective ionization parameter of the atmosphere. The actual differences between these two reflection spectra in the energy range of interest are fairly small, with the most significant change being narrower lines due to the O emission region lying within 0.5 Thomson depths of the surface.. Only a single model with solar abundances was computed with these parameters.

## 3 Results

### 3.1 Constant density models

It is expected that for most Seyfert 1s, the X-ray power-law will be observed directly along with the reflection spectrum, so our analysis concentrates on the sum of the incident and reflection spectra (i.e., a reflection fraction of one). In Figure 2 we plot the resulting nine total spectra between 0.3 and 1.2 keV.

The plots show many lines in this region of the spectrum which vary in shape and strength as the parameters are changed: C vi Ly$`\alpha `$ at 0.37 keV, N vi (resonance+intercombination) at 0.43 keV (usually very weak), N vii Ly$`\alpha `$ at 0.50 keV, O vii (resonance+intercombination) at 0.57 keV, O viii Ly$`\alpha `$ at 0.65 keV, and a cluster of Fe L lines (particularly from Fe xviixix) between $``$0.7 and 1.0 keV. In Table 1 the equivalent widths (EWs) of the most prominent lines are listed for each of the nine models, as well as the width of their blue wings (WBW).

The EWs were calculated by integrating the spectrum over the line between two hand-picked values of the continuum on either side. The WBW is simply the difference between the continuum on the high-energy side of the line and the value of the spectrum at the line peak, and is the minimum width of the line. Since the continuum can be difficult to determine at times (particularly in the cases with high O abundance), these measurements will not necessarily be very precise, but will accurately show trends in the data.

Concentrating on the models with solar abundances (the left-hand column in Fig. 2), we see that the O viii line is the most prominent soft X-ray emission feature at all three ionization parameters with a maximum EW$``$73 eV when $`\xi =250`$ erg cm s-1. As $`\xi `$ is increased to 1000 erg cm s-1, the N and Fe L lines have all but disappeared, and the O viii line is drastically weaker (EW$``$31 eV) with broad wings on either side of the line core due to Compton scattering (WBW$``$60 eV).

According to Table 1 all of the lines have a WBW on the order of tens of eV, even at lower ionization parameters. This is because the majority of the O viii line emission occurs at a non-negligible Thomson depth and Compton scattering is a major contributor to the shape of the line. As an illustration, in Figure 3 the oxygen ionization fractions in the $`\xi =500`$ erg cm s-1 model are plotted as functions of Thomson depth.

In this case, large fractions of O viii and O vii are produced only at Thomson depths $`1`$. Therefore, the line emission from these species must pass through an optically thick scattering layer before escaping. As $`\xi `$ increases the lines originate from a larger Thomson depth and the photons are scattered more often.

Increasing the O abundance by a factor of five primarily increases the strength of the O vii line. This is because the incident spectrum is held constant and the number of possible ionizations from O viii to O ix is fixed. With the increase in O abundance, it is no longer possible to fully-strip a large zone of O, and the curves seen in Fig. 3 are shifted to the left. Furthermore, since all of the ionizing photons have been used up, a larger number of oxygen ions can recombine to the O vii configuration. Therefore, the O vii line strongly increases with the oxygen abundance. In fact, at $`\xi =250`$ erg cm s-1 the EWs of both the O vii and O viii lines are $``$260 eV which would hopefully be strong enough to be detected in a high signal-to-noise spectrum of a Seyfert 1.

In contrast to O, when N is increased by a factor of five, only the N vii line grows significantly. This is due to the fact that, even with a five times overabundance, N is not abundant enough for the N vii edge to saturate in these conditions. The N vii line is strongest in the $`\xi =250`$ erg cm s-1 N-overabundant model, but still has an EW of only 21 eV. In all other cases, the EW of N vii is $`20`$ eV. We conclude that at the current time it is unlikely for N vii Ly$`\alpha `$ emission to be observed from an irradiated accretion disc.

The C vi line is usually the second strongest line in the reflection spectrum with EWs between 10 and 30 eV. In the solar abundance model the C vi EW only decreases by about a factor of two as the ionization parameter is increased from 250 erg cm s-1 to 1000 erg cm s-1. These results certainly suggest that the C vi, O vii and O viii lines are the best candidates for observable soft X-ray accretion disc lines.

### 3.2 Hydrostatic model

The summed reflection and incident spectrum from the hydrostatic model is shown in the left panel of Figure 4.

Qualitatively, the spectrum looks similar to the $`\xi =500`$ erg cm s-1 constant density model with O viii Ly$`\alpha `$ dominating the emission. However, many of the lines have broader bases than in the analogous constant density spectrum. The broad bases are particularly apparent for the C vi (WBW=43 eV), and Fe L lines. Furthermore, as seen in Table 2, the lines are actually quite weak with the EW of O viii only reaching 17 eV. These values are in good agreement with the results of similar computations done by Różańska et al. (2002) and Nayakshin (private communication).

The reason for these differences is that the more realistic density structure results in a hot scattering layer overlying the line emitting material. In the right-hand panel of Figure 4 the vertical temperature profile of the hydrostatic model is overlaid on a plot of the oxygen ionization structure. We see that the temperature reaches $`kT`$1 keV at the surface of the atmosphere, which is nearly an order of magnitude hotter than the surface of the $`\xi =500`$ erg cm s-1 constant density model. Oxygen is fully stripped until a Thomson depth of 0.9, where, because of the rapidly increasing gas density, it quickly recombines to O iii, the lowest ionization state that is treated. Most of the O viii emission therefore originates from a narrow zone which is not very dense. Thus the O viii line is intrinsically weaker in a hydrostatic disc model. The region in which large fractions of O viiviii are produced is so restricted that the O vii line is very much weaker in the hydrostatic model than in the constant density models. Therefore, the hydrostatic models predict that only O viii is detectable.

The existence of the hot, fully ionized skin on the surface of the disc dilutes the spectral features which are formed underneath by Compton scattering the emission lines. This results in the broad bases seen in Figure 4. Although most of the O viii emission originates from only a single Thomson depth into the atmosphere, the temperature of the ionized skin is such to result in significant dispersion. An ionized skin at the surface of the atmosphere is a common feature of hydrostatic disc models (Nayakshin et al., 2000; Ballantyne et al., 2001). If the illuminated hydrostatic models are a more realistic description of accretion discs than the constant density ones, then soft X-ray accretion disc lines are expected to be weak and broad features, and will therefore be challenging to detect.

## 4 Discussion and Implications

The ultimate goal of X-ray spectroscopy of AGN is to learn about accretion flows around supermassive black holes. Detecting soft X-ray emission lines from the accretion disc would complement the information inferred from the Fe K$`\alpha `$ line. In principle this should be relatively straightforward as the effective area of most X-ray detectors peak in the soft X-ray band, and gratings observations would allow an accurate determination of any warm absorber. Unfortunately, the calculations presented in the previous section suggest that most of the emission lines will be quite weak and broad (even before relativistic blurring), and therefore difficult to pull out of the noise with current instruments222If, as the recent numerical simulations of Turner, Stone & Sano (2002) suggest, the disc is clumpy or inhomogeneous, then it is much more difficult to predict the properties of the reflection spectrum.. The one exception may be the O viii Ly$`\alpha `$ line which has an EW that can approach 100 eV, depending on the density structure of the disc. If that line could be measured accurately enough then, as seen in Fig. 2, it would give a good constraint on the ionization parameter of the disc. If significant O vii emission is also seen then abundance information may also be derived. Although the unambiguous detection of individual lines will likely be a challenge, soft X-ray accretion disc lines have the potential to be important to the study of AGN.

### 4.1 The soft X-ray features of MCG–6-30-15

Recently, Branduardi-Raymont et al. (2001) and Sako et al. (2002) claimed to have detected relativistically blurred C vi Ly$`\alpha `$, N vii Ly$`\alpha `$ and O viii Ly$`\alpha `$ lines in XMM-Newton RGS observations of the Seyfert 1 galaxies MCG–6-30-15 and Mrk 766. This interpretation was challenged by Lee et al. (2001) who favour a dusty warm absorber model to explain the features seen in their Chandra gratings data of MCG–6-30-15. The predictions of ionized disc models should be able to provide a data-independent view of which interpretation is most probable. Ballantyne & Fabian (2001) fitted the ionized disc models of Ross & Fabian (1993) and Ballantyne et al. (2001) to the ASCA data of MCG–6-30-15, concentrating on the region around the Fe K$`\alpha `$ line. The resulting best fit models predicted O viii, O vii and Fe L emission with EWs $`<`$ 20 eV, incompatible with the results of Branduardi-Raymont et al. (2001).

The results presented in Sect. 3 enable us to take an orthogonal point of view. Ignoring the Fe K$`\alpha `$ line entirely, we can ask the question: is there a region in the explored parameter space which can account for the putative soft X-ray emission lines? The updated analysis of MCG–6-30-15 by Sako et al. (2002) is very specific: they only see emission lines of O viii, N vii and C vi with EWs of $`162\pm 8`$ eV, $`54.4\pm 3.2`$ eV and $`24.9\pm 2.5`$ eV, respectively. From Table 1 it seems that a constant density model with $`\xi =250`$ erg cm s-1 and an O overabundance could account for the C and O lines, but none of the models can account for the strength of the N line, unless there is also a large (i.e., greater than 5 times solar) overabundance of N. However, these models predict significant O vii and Fe L emission which is not seen by Sako et al. (2002), although detectable given the limits placed on the other lines. The only way to remove these features would be to increase the ionization parameter, but this will weaken and smear out the O viii line. Of course, it will be even more difficult to obtain the required line strengths with hydrostatic models because the emission features are intrinsically weaker.

Our calculations point to another serious problem with the relativistic emission line interpretation. Namely, the blue wing of the putative O viii line is very sharp, dropping down at 700 eV to the continuum in $`\mathrm{\Delta }E10`$ eV, as judged from the high resolution Chandra data where the drop is $`<`$3 eV (Fig. 1 in Lee et al. 2001). The line must therefore have a relative width of much less than one per cent. Even a single scattering from an electron at 106 K will impart a rms width of 2 per cent to the line, rising to 6 per cent at 107 K. A sharp observed feature from a line must therefore be due to the unscattered line core333Broadening by electron scattering is relevant to all narrow emission lines produced by reflection, and could strongly modify predictions, e.g. those of Li, Gu & Kahn (2001)., which usually carries little flux.

The ionized disc models predict a minimum width of the O viii line of over 40 eV (Table 1). Relativistic blurring makes the width even larger. Figure 5 shows the results of blurring the $`\xi =250`$ erg cm s-1solar abundance model (the top-left model from Fig. 2) with the parameters given by Sako et al. (2002).

With the relativistic blurring effects, the O vii and O viii lines are combined together and peaks at 0.6 keV. At 0.7 keV the blue wing of this blended line falls by only 20 per cent over a range of 50 eV. This is significantly weaker than the factor of 2–3 drop over $``$3 eV seen in the Chandra and XMM-Newton data. The 0.7 keV drop is weaker in this reflection model because the Fe L lines have been blurred together and redshifted against the edge of the oxygen lines. The oxygen lines can be made sharper by moving to lower $`\xi `$, but only at the expense of increasing the strength of the O vii and Fe L lines. As before, hydrostatic models make the situation worse because they provide additional Compton scattering due to the ionized skin. Of course, these considerations also apply to Mrk 766.

Although we have obviously not covered the entire available parameter space, we conclude that it is very difficult, if not impossible, for reflection from ionized discs to produce the lines required by the interpretation of Branduardi-Raymont et al. (2001) and Sako et al. (2002). First, there does not seem to be a situation where one can obtain strong O viii, N vii and C vi Ly$`\alpha `$ emission without simultaneously having significant O vii and Fe L lines. Second, in any situation where O viii does dominate the emission spectrum, the disc has to be fairly highly ionized and so the lines are Compton broadened with widths of 10s of eV. This width will further increase upon relativistic blurring, and therefore cannot possibly have the sharp, $`10`$ eV drop at the high-energy side that the MCG–6-30-15 data requires. More realistic hydrostatic ionized disc models exacerbate these problems because of the hot ionized skin on the surface of the disc. In general, any soft X-ray relativistic line will tend to create broad weak features in the observed spectrum. Any large sharp drops in this spectral range must be due to absorption edges.

It is likely there are weak broad features due to ionized oxygen in the spectrum of MCG–6-30-15. They are a plausible identification for the soft excess required in spectral modelling of Chandra data by Lee et al. (2001) and for the residuals from a dusty warm absorber fit to XMM-Newton data by Sako et al. (2002, see their Fig. 4). On the basis of our calculations, we conclude that the sharp drop seen in all soft X-ray spectra of MCG–6-30-15 at 0.7 keV cannot plausibly be the blue wing of an O viii emission line produced by reflection.

## Acknowledgements

The authors thank the referee, Masao Sako, for his helpful comments and suggestions. DRB acknowledges financial support from the Commonwealth Scholarship and Fellowship Plan and the Natural Sciences and Engineering Research Council of Canada. ACF and RRR acknowledge support from the Royal Society and the College of the Holy Cross, respectively.