text:
# Testing and forecasting the time series of the solar activity by singular spectrum analysis

## 1 Introduction

It has long been observed that solar activity depends on a number of spots visible on its disk. During about 11 years which is called a solar cycle, this number varies over a wide range. Accompanying this process changes in the structure of magnetic fields of the Sun affect the Earth climate and have a probable connection with the natural catastrophes. In addition, intensity of the solar radiation (the frequency of the Sun flares, etc.) apparently exerts an influence on all areas of human life including socially-historical activity. Thus, in view of the quite large significance of the magnetic activity of the Sun its prediction is a subject of much current interest.

At the present time, to describe the dynamics of the solar activity many approaches are used. Among them the Wolf number associated with the effective amount of the sunspots is much more convenience method. The dynamics of the Wolf numbers has more or less quasiperiodic nature but, owing to the facts that models of this process do not take into account many essential factors of the solar magnetic activity, its prediction is difficult. It is known that during last 250 years the cycle period of the sunspots changes its value not more than $`20\%`$. At the same time, the amplitude (i.e. the averaged number of sunspots) has varied over an order or even more. Detailed models of the solar processes do not describe such variations.

In the last few years sufficiently many methods devoted to the prediction and reconstruction (to the past) the dynamics of the Wolf numbers have been proposed (see, e.g.,  and references cited therein). It should be noted however, that these methods have certain limitations and demerits. That is the reason why prediction of the sunspot dynamics based on the data of observations only (i.e. without a modelling of the process) is a quite perspective approach. In this way time-series analysis (see ) can give an essential assistance. But in this case there are many problems related to the fact that the Wolf series is not strictly deterministic system, it does not have clearly defined dimension , and its length is not so large.

In the present paper, for investigations of the Wolf time series a sufficiently efficient method of singular spectrum analysis is proposed. It is shown that this method provides a quite high reliability in the description of amplitudes of the 11-year solar cycle. Moreover, it is appropriate for revealing more long cycles and forecasting the further solar activity during one and a half of 11-year cycle, i.e. the nearest 16 years.

## 2 Singular spectrum <br>analysis

The method of singular spectrum analysis (SSA) , which is used in the present article allows us the following.

* Recognise certain components in the time series which have been obtained from the observable at regular intervals;
* Find periodicities that are not known in advance;
* On the basis of the chosen components to smooth out the initial data;
* Extract components with the known period;
* Predict the further evolution of the observed dependence.

The method of SSA is a reasonably new one, but now it is quite clear that it is sufficiently competitive with numerous smoothing methods . Moreover, in certain cases forecasting the system evolution on its basis gives much more reliable results in comparison with the other known algorithms (see  and references therein).

The SSA method is based on the passage from the investigation of an initial linear series $`(x_i)_{i=1}^N`$ to the analysis of a many-dimensional series consisted of components of some length which contain except for the value $`x_i`$, certain quantities $`x_{ij}`$, $`j=1,\mathrm{},\tau `$, at the previous instants of time.

Let us describe the central steps of the application of SSA to a series $`(x_i)_{i=1}^N`$.

(1) On the first step, a one-dimensional series is transformed into a many-dimensional one. For such a transformation it is necessary to take a certain number of delays $`\tau [N/2+1]`$, where $`[]`$ is an integer part, and represent the initial $`\tau `$ values of the sequence as the first column of some matrix $`X`$. For the second column, the values of the sequence from $`x_2`$ to $`x_{\tau +1}`$ are chosen. Thus, the last $`\tau `$ elements of the sequence $`x_n,\mathrm{},x_N`$ correspond to the last column with the number $`n=N\tau +1`$. Therefore the transformed series has the following matrix form:

$$X=\left(\begin{array}{ccccccc}x_1& x_2& \mathrm{}& x_\tau & \mathrm{}& x_n& \\ x_2& x_3& \mathrm{}& x_{\tau +1}& \mathrm{}& x_{n+1}& \\ x_3& x_4& \mathrm{}& x_{\tau +2}& \mathrm{}& x_{n+2}& \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}& \\ x_\tau & x_{\tau +1}& \mathrm{}& x_{2\tau 1}& \mathrm{}& x_N& \end{array}\right).$$

It is obvious that for the constructed matrix the expression $`x_{ij}=x_{i+j1}`$ holds. In general, the matrix $`X`$ is a rectangular one. But, in a limit case, i.e. at $`\tau =N/2`$ and an even $`N`$, it degenerates into a square matrix.

(2) After this transformation, for the matrix $`X`$ the corresponding covariance matrix

$$C=\frac{1}{n}XX^T$$

should be obtained.

(3) Now it is necessary to find eigenvalues and eigenvectors of the matrix $`C`$. For this, the matrix $`C`$ should be factored as follows: $`C=V\mathrm{\Lambda }V^T`$, where $`\mathrm{\Lambda }=\left(\begin{array}{cccc}\lambda _1& 0& \mathrm{}& 0\\ 0& \lambda _2& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ 0& 0& \mathrm{}& \lambda _\tau \end{array}\right)`$ is the diagonal matrix of eigenvalues, $`V=(V^1,V^2,\mathrm{},V^\tau )=\left(\begin{array}{cccc}v_1^1& v_1^2& \mathrm{}& v_1^\tau \\ v_2^1& v_2^2& \mathrm{}& v_2^\tau \\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ v_\tau ^1& v_\tau ^2& \mathrm{}& v_\tau ^\tau \end{array}\right)`$ is the orthogonal matrix of eigenvectors of the matrix $`C`$. It is clear that $`\mathrm{\Lambda }=VCV^T`$, $`_{i=1}^\tau \lambda _i=\tau `$ and $`detC=_{i=1}^\tau \lambda _i`$.

(4) For the next step, the matrix $`V`$ of eigenvectors should be presented as a matrix of the conversion to the principal components, $`Y=V^TX=(Y_1,Y_2,\mathrm{},Y_\tau )`$, of the initial series. Here $`Y_i`$, $`i=1,2,\mathrm{},\tau `$, are rows of the length $`n`$. Therewith, the eigenvalues $`\lambda _1,\lambda _2,\mathrm{},\lambda _\tau `$ can be considered as a certain contribution of the principal components to a general information content of the series $`(x_i)_{i=1}^N`$. Then, by means of these principal components it is possible to reconstruct the initial matrix $`X`$:

$$X=(V^1,V^2,\mathrm{},V^\tau )\left(\begin{array}{c}Y_1^T\\ Y_2^T\\ \mathrm{}\\ Y_\tau ^T\end{array}\right)=\underset{i=1}{\overset{\tau }{}}V^iY_i^T.$$

In turn, by the matrix $`X`$ one can reconstruct the time series $`(x_i)_{i=1}^N`$. It should be noted that for the reconstruction not all principal components $`Y_1,Y_2,\mathrm{},Y_\tau `$ are usually applied. Only a part of them can be involved. This depends on the goal which we pursue and the informative content of the used components (see ). This means that each vector-row $`Y_i`$ can be considered as a result of some projection of a $`\tau `$-dimensional totality on a direction corresponding to the eigenvector $`V^i`$. Thus, the series is presented via a set of $`\tau `$ components $`Y_i`$. Therewith, the weight of the component $`Y_i`$ in the initial sequence $`(x_i)_{i=1}^N`$ can be defined by the corresponding eigenvalue $`\lambda _i`$ which is, in turn, the eigenvalue of the eigenvector $`V^i`$.

Each $`i`$-th eigenvector includes $`\tau `$ components, $`V^i=\left(\begin{array}{c}v_1^i\\ v_2^i\\ \mathrm{}\\ v_\tau ^i\end{array}\right)`$. Let us construct a dependence of the component values $`v_k^i`$, $`k=1,2,\mathrm{},\tau `$, as a function of their number: $`v^i=v^i(k)`$. Then using the orthogonality property of eigenvectors, the further analysis of the sequence $`(x_i)_{i=1}^N`$ can be performed by means of diagrams plotted by the analogy with the Lissajous figure. Namely, along the axes the components $`v_k^i`$, $`v_k^j`$ are plotted in pairs. If the constructed diagrams are similar to a circle, then the functions $`v^i=v^i(k)`$, $`v^j=v^j(k)`$ are approximated by certain periodic functions with almost coincided amplitudes and a phase lag about a quarter of the period.

Thus, for some pairs of the eigenvectors $`V^i`$, $`V^j`$ one can find a value which has a meaning of the period. Therefore, such a graphical representation provides a certain estimate of the component frequencies inherent in the initial time series $`(x_i)_{i=1}^N`$.

For the given parameter $`\tau `$ the number of all possible pairs of principal components is $`\tau ^2`$. It is obvious that even for a sufficiently small $`\tau `$ the analysis of all such pairs is a quite laborious problem. Moreover, at a large values of $`\tau `$ only a small part of diagrams has a helical form. Thus, before a graphical analysis it would be reasonable to restrict our search. This can be done if we arrange $`V^i`$ and $`Y_i`$ in order of decreasing their eigenvalues and consider only such pairs of eigenvectors which have close enough eigenvalues $`\lambda `$. In the diagram $`\lambda =\lambda (i)`$ at a quite large $`\lambda `$ these pairs form a decreasing (with the growth of $`i`$) step function.

(5) Suppose now that for the further reconstruction we have only $`r`$ leading components. Thus, for the reconstruction of the initial matrix $`X`$ one should use $`r`$ leading eigenvectors $`V^i`$. In this case

$$\stackrel{~}{X}=(V^1,V^2,\mathrm{},V^r)\left(\begin{array}{c}Y_1\\ Y_2\\ \mathrm{}\\ Y_r\end{array}\right)=\underset{i=1}{\overset{r}{}}V^iY_i,$$

where $`\stackrel{~}{X}`$ is a reconstructed matrix with $`n`$ columns and $`\tau `$ rows. Then the initial time series obtained from this matrix is defined as follows:

$$\stackrel{~}{x}_s=\{\begin{array}{cc}\frac{1}{s}\underset{i=1}{\overset{s}{}}\stackrel{~}{x}_{i,si+1},1s\tau ,\hfill & \\ \frac{1}{\tau }\underset{i=1}{\overset{\tau }{}}\stackrel{~}{x}_{i,si+1},\tau sn,\hfill & \\ \frac{1}{Ns+1}\underset{i=1}{\overset{Ns+1}{}}\stackrel{~}{x}_{i+sn,ni+1},\hfill & \\ nsN.\hfill & \end{array}$$

Described way of the reconstruction is said to be a SSA-smoothing of the initial time series $`(x_i)_{i=1}^N`$ by the leading $`r`$ components.

(6) In the next stage of the SSA application, a prediction procedure of the initial time sequence can be considered (see ). This means that the series $`(x_i)_{i=1}^{N+p}`$ which is the extension of the known data $`(x_i)_{i=1}^N`$, is constructed. In turn, extrapolation to $`p`$ points forward is reduced to the application of $`p`$ times of the prediction procedure to the one point. The basic idea of the computation of the point $`x_{N+1}`$ is the following.

Consider the sequence $`x_1,x_2,\mathrm{},x_N`$ and construct a sample in the form of matrix $`X`$. As a basis of the surface containing this sample one can take the chosen before vectors $`V^1,V^2,\mathrm{},V^r`$ of the matrix $`C`$. Let us write the parametric equation of this surface as $`S(P)=\underset{i=1}{\overset{r}{}}p_iV^i`$, where a set of the parameters $`p_i`$ corresponds to the value $`S(P)`$ which is a column with $`\tau `$ elements. In this case the set of parameter values $`P^k=(p_1^k,p_2^k,\mathrm{},p_r^k)`$ corresponds to $`k`$-th, $`k=1,2,\mathrm{},n`$, column of the matrix $`X`$. Therefore, $`X^1=S\left(P^1\right)`$, $`X^2=S\left(P^2\right)`$, $`\mathrm{}`$ , $`X^n=S\left(P^n\right)`$. Now, to predict the value $`x_{N+1}`$ it is necessary to find the $`(n+1)`$-th column $`X^{n+1}`$ which, in turn, fits the parameters $`P^{n+1}=(p_1^{n+1},p_2^{n+1},\mathrm{},p_r^{n+1})`$. Using the data $`(x_i)_{i=1}^N`$ these parameter values can be obtained from the expression $`S(P)=\underset{i=1}{\overset{r}{}}p_iV^i`$. Thus, the predicting column is written as follows: $`X^{N+1}=S\left(P^{n+1}\right)`$.

Let us introduce the following designations:

$$V_{}=\left(\begin{array}{cccc}v_1^1& v_1^2& \mathrm{}& v_1^r\\ v_2^1& v_2^2& \mathrm{}& v_2^r\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ v_{\tau 1}^1& v_{\tau 1}^2& \mathrm{}& v_{\tau 1}^r\end{array}\right);$$

$$\stackrel{~}{p}=\left(\begin{array}{c}\stackrel{~}{p}_1^{n+1}\\ \stackrel{~}{p}_2^{n+1}\\ \mathrm{}\\ \stackrel{~}{p}_r^{n+1}\end{array}\right);Q=\left(\begin{array}{c}x_{N\tau +2}\\ x_{N\tau +3}\\ \mathrm{}\\ x_N\end{array}\right);$$

$$V_\tau =(v_\tau ^1,v_\tau ^2,\mathrm{},v_\tau ^r).$$

The set of the values $`(p_1^{n+1},p_2^{n+1},\mathrm{},p_r^{n+1})`$ is easily found as a solution of the system $`V_{}\stackrel{~}{P}=Q`$ with respect to $`\stackrel{~}{P}`$. Thus, the final expression of the forecasting value has the following form:

$$x_{N+1}=\frac{V_\tau V_{}^TQ}{1V_\tau V^T}.$$

In the simplest case, to predict the next values it is necessary to change the matrix $`Q`$ in a corresponding way and multiply it by the value $`V_\tau V_{}^T/\left(1V_\tau V^T\right)`$. In addition, however for each predicted point one can completely repeat the SSA algorithm. Then the matrixes $`V_\tau `$ and $`V_{}`$ will be changed.

(7) At the final stage of the SSA application one should dwell on the choice of the main parameter — the number of delays $`\tau `$ which are used for the many-dimensional sample $`X`$. As in the case of selection of the principal components this value essentially depends on the investigated problem.

Consider the smoothing procedure of a time series by the SSA method. In this case, as noted above, selection of a principal component is the filtration of the time series with the transition filtering function in the form of the eigenvector of this principal component. If the delay value $`\tau `$ is the greater, the greater the number of parallel filters, and a bandwidth of each of them is more narrow. Thus, for a large enough $`\tau `$ we have a sufficiently efficient smoothing of the time series.

If it is necessary to define unknown (hidden) periods in the observed sequence then one should take the value of $`\tau `$ as large as possible. Next, after omitting close to zero eigenvalues, the delay value should be shortened.

Suppose now that we should select only the one known periodic component. In this case it is necessary to choose the delay $`\tau `$ which is equal to the required period.

Finally, let us consider the problem of some extension of an observed sequence to a given value, i.e. the problem of forecasting the evolution of the process under investigation. Then one should take a maximum allowable value of the delay $`\tau `$ and thereafter select the number $`r`$.

## 3 Forecasting the solar activity by SSA

To estimate of the practicality of the SSA method it is necessary to use the time series of a natural origin. In the present article, as such a time series the sequence of the Wolf numbers characterised in a certain way the solar activity was chosen. Taken alone, the Wolf numbers defined by the visible sunspots cannot give a quantitative information related to the solar activity. However there is a large enough correlation between the Wolf numbers and the $`F10.3`$ emission. That is the reason why investigations of the relative variations in this sequence can give a specific information concerning the solar activity.

For the first time, in 1848 a Swiss astronomer R.Wolf proposed that a measure of the solar activity can be characterised by the number of sunspots. To this end he recommended to consider the union of the total number of spots visible on the Sun and tenfold numbers of regions in which these spots are placed. This last summand should coordinate the results of measurements performed under several conditions. Thus, since the year 1849 the results of daily measurements have become available. Using more earlier observations and various sources, R.Wolf has reconstructed the data of the solar activity (with an admissible accuracy and negligible gaps) up to the year 1818. Now the averaged number of the sunspots is called the Wolf number.

Later the average monthly values of the Wolf numbers up to the year 1749 (namely this series is used in the present paper) and their average annual values up to the year 1700 have been reconstructed. In the last case however, the error can reach more than 10 percent. The chosen data covers the wide time interval without gaps and a quite high time resolution. The investigated sequence is shown in Fig.1. The time with the interval of one month is plotted on the abscissa. The corresponding Wolf numbers are laid as ordinates.

On the first step of the application of SSA it is necessary to take the maximum permissible value of parameter $`\tau `$. For our investigations we used $`\tau =500`$. This value allows us to make up periodicities up to period of 42 years. The use of the larger values leads to essential computer problems. Moreover, increasing the parameter $`\tau `$ (up to 600) does not yield any significant change in the results of the first expansion in principal components.

Owing to a large value of $`\tau `$ the dependence of the roots of eigenvalues of the covariance matrix (ordered in decreasing) decays exponentially. In addition, the sum of the leading five eigenvalues exceeds $`99\%`$ of their sum total. Coupled with a quite large numbers of initial points, this leads to the fact that the first principal component yields a small enough smoothing of the initial series, and by means of the leading four-five components one can reconstruct this series. Moreover, at small values of the parameter $`\tau `$, say at $`\tau =5`$, the first component weakly changes its form. This is due to the fact that the SSA method is stable with respect to this parameter. Therefore, application of a sufficiently large $`\tau `$ is justified only for the prediction.

To check possibilities of the prediction let us cut off the initial sequence of the averaged monthly values for a length of 216 points (18 years) and reconstruct it by the following way. Evaluate optimal parameter values of the reconstruction algorithm by means of an additional truncation and resolve the obtained series into $`\tau `$ components. Therewith it is necessary to choose such a number of the leading components at which we get the best coincidence of the reconstructed values with the additionally truncated data. Then using the obtained parameters, we will find the initially shortened part of 216 points.

By means of the direct exhaustion one can define that the best results we get at $`r=150`$ (the number of the chosen components). Once again let us take the reduced sequence (for a length of 18 years) and apply the chosen $`r`$ for the prediction. One can further improve the quality of the prediction if we will decompose the predicted interval (216 points) into subintervals and after prediction of each part calculate the principal components. In the ideal case, after the prediction of each point it is necessary to realise this process. However, such an approach leads to an increase of the computer time.

Fig.2. illustrates the reconstructed series at which recalculations have been realised every 72 points (i.e. 3 times). However, as it turns out, this result is almost identical to the prediction of the whole interval of 216 points, without decomposition into subintervals.

Theoretically, one can analyse the series components with the aim of detection the existence of the known and any other periods. However a number of these periods and thus, the resemblance of the components $`Y_i`$ of the initial series makes this analysis a quite cumbersome procedure. In addition, selection of sufficiently short periodicities (about several months) is a very difficult computer task. For these reasons, the choosing a large enough subinterval of time is the more simply problem.

Let us consider now a series with the average annual Wolf numbers. This series contains only 248 points. Thus, the maximum allowable value of $`\tau `$ is 123. Choose this value as an initial one. The corresponding eigenvalues is shown in Fig.3. The first number presents the principal component related to a trend. The following step-like data form the pairs of components with the numbers 2–3, 4–5, 6–7, 8–9 and 11–12. Beginning with 14-th number this dependence gives way to an exponential tail. Eigenvectors for the pairs 2–3, 4–5, 8–9 and 11–12 (see Fig.4) fit periodicities with 11-year period. The corresponding helical dependence for the components 2 and 3 is shown in the left side of Fig.5.

It should be noted that except for the evident the eleven-year solar cycle it is possible to guess a supposed eighty-year Gleisberg cycle (see Fig.5). Here we keep in mind the pair of the eigenvectors 6 and 7. The corresponding eigenvalues is not exactly equal to each other, i.e. the step-like data is oblique (see Fig.3) and the phase lag is not $`\pi /2`$. That is the reason why the form of the function diagram is not helical. In spite of this fact and a small enough eigenvalues, it is quite possible to trace a periodicity.

For the best selection of the eighty-year dynamics one can identify parameter $`\tau `$. By the numerical analysis we have found that $`\tau =80`$ is the most suitable value. In this case, vectors 4 and 7 are associated with such a periodicity. The diagram for 4-th and 7-th components is shown in Fig.6. A possible eighty-year solar cycle obtained by the reconstruction only via these two components is shown in Fig.7.

Consider now a much more interesting problem concerning the possibility of the prediction of the average annual Wolf numbers by the SSA method. Shorten this series from the right hand side for a length of 18 points (that means years) and reconstruct it. To make this let us truncate additionally an interval of 11 points and try to restore it in the best way. Choosing the suited numbers of eigenvectors we will use such a procedure in various parts of the series.

For the series of 219 points ($`219=2481811`$, where 248 is the number of all points, see above) the maximal possible $`\tau =109`$. As follows from numerical analysis, prediction for 11 points strongly depends on the picking of the components. At the same time, in a sufficiently wide range a qualitative guide of the prediction is not bad. However, for the quantitative prediction the given $`\tau `$ is too large. For the lesser values of $`\tau `$ it is quite possible to find the necessary number of components. For example, for $`\tau =33`$ and choosing the leading 11-th components we get fairly good results of the prediction.

Thus, let us use these parameters for the prediction of the removed 18 points. In contrast to the case of the average monthly data and $`\tau =500`$ now, at $`\tau =33`$ we have a possibility to recalculate eigenvectors with consideration of the last predicted point. This recalculations can be made at the stage of the selection of vectors as well as the prediction. The result of the prediction with the correction at each step is shown in Fig.8.

In addition, we have studied a sequence of natural logarithms obtained from the initial series. Taking the logarithm is often used at the data processing (for example, in the case of the correlation analysis) and allows to get more interesting and complete results. However, eigenvalues and eigenvectors after such an operation have not principally changed. Thus, in this case taking the logarithm is not necessary because the basic information can be obtained from the analysis of the initial series.

In the closing stage let us consider application of SSA for a real prediction of the solar activity. For this purpose, the average annual sequence of the Wolf numbers from the year 1748 to the year 1996 has been chosen. The end of this series corresponds to a minimum of the solar activity. Therefore it is interesting to describe the further activity of the Sun and predict its two next maxima. To realise this idea it is necessary to extend the average annual sequence for a length of 18 year points.

We resolved the series into 33 components and choose for the prediction the leading eleven of them. The result of the prediction is shown in Fig.9. As follows from this figure, in the nearest future, in comparison with the two previous maxima, the Sun will be in a relatively quiescent state. In addition, the level of the forthcoming maximum in the year 2011 will not be so high. On the basis of our numerical investigations concerning the application of SSA, the prediction during more than one and a half of 11-th year solar cycle is not justified. Nevertheless, we believe that the obtained values of the Wolf numbers shown in Fig.9 are a quite reliable.

For comparison, a series of the average monthly Wolf numbers has been analysed. However, its investigation did not yield an essential modifications to the obtained results.

## 4 Conclusion

It seems to be true that in the nearest future investigations of time series by means of SSA will occupy a more important and deserving place among different ways related to the processing and forecasting many experimentally obtained sequences. Resolving the initial series into components the analytical form of which is not fixed, this method allows us on a quite high level to reveal periodic subsequences and forecast the dynamics of this series. Therewith, restrictions to the numbers of points and characteristic periods in the investigated data as a rule, is essentially less than at applications of the other methods (for example, at correlation or Fourier analyses).

In the present article a possibility of the SSA application to the sequence of the Wolf numbers characterising the solar activity, is considered. In spite of a relatively small length of this sequence, SSA allows to reveal the components corresponding the known solar cycles and, on the basis of only certain constituents, admits reconstruction of the data. Moreover, it is found that by means of SSA it is possible to extend short time series with an acceptable accuracy.

Like any other method, SSA is not devoid of drawbacks. First, there is a certain difficulty in the problem of revealing unknown frequencies in the investigated sequence. These frequencies can be more easily obtained by the Fourier analysis. Secondly, SSA does not contain clear rules concerning the choice of components, especially in the case of forecasting. Finally, SSA does not give the correct prediction of the cycle period. Therefore in the predicted sequence a systematic phase lag is accumulated. Nevertheless, as follows from the performed analysis it is a useful supplement to the existing methods of the experimental data processing.

In the forthcoming investigations devoted to the analysis of the solar activity, we will consider a possibility of a certain correction of the cycle phase by means of additional methods and via special feature of the sunspot dynamics. Besides, to resolve this problem one can use an empirical dependence of the amplitude and the phase of the solar cycle . Moreover, it seems to be reasonable that it is quite possible to improve control parameters which are used in SSA.

On the whole we can say that the described method of singular spectrum analysis is a sufficiently advantageous and promising way for the prediction of dynamics of the solar activity.


format: text

text:
# Localization on short-range potentials in dissipative quantum mechanics

\[

## Abstract

In this article the problem of the existence of a state localized on a weak short-range attractive potential in the presence of dissipation is considered. It is shown that, contrary to the pure quantum case, a localized state is produced in any number of dimensions, while in low dimensions dissipation leads to much stronger localization. The results have physical implications for the dissipative dynamics of objects such as heavy particles in Fermi liquids and for superconductivity in high-$`T_c`$ materials.

\]

Research in the theory of strongly interacting fermionic systems was rekindled soon after the discovery of high temperature cuprate superconductors when it was realized that Landau Fermi liquid theory can not describe the normal state of the underdoped materials. Since then, many ideas have been proposed to exhibit possible routes to non-Fermi liquid behavior.  While these theories differ in details, they all are intimately related to a peculiar property of fermionic liquids known as Anderson’s orthogonality catastrophe.  This property implies that the wave function of the liquid before and after an adiabatic insertion of an impurity has a vanishing overlap:

$$\psi |\psi ^{}N^A,$$
(1)
where $`N`$ is the number of fermions in the liquid and $`A`$ is the orthogonality exponent which can be expressed in terms of impurity scattering phase shifts. This result shows that interactions between fermions have singular effect and the fact that the non-interacting quasiparticles in Landau theory can often be defined is rather an exception than a principle.

Kondo  introduced a closely related overlap of the wave function of the liquid with an impurity at 0 and that with the impurity moved slowly to $`𝐱`$:

$$\psi (0)|\psi (𝐱)(\mathrm{\Lambda }t)^{K(𝐱)},$$
(2)
where $`t`$ is a characteristic time of the impurity motion and $`\mathrm{\Lambda }`$ is a high energy cut-off. The exponent $`K(𝐱)`$ is quite difficult to evaluate in general; however, it was shown by Schönhammer  that for small $`𝐱`$ ($`\mathrm{}=1`$ throughout the letter):

$$K(𝐱)\frac{\eta 𝐱^2}{2\pi }.$$
(3)
Here $`\eta `$ is the friction coefficient that determines the viscous force $`𝐅=\eta \dot{𝐱}`$ exerted on the moving impurity by the fermions in the liquid. That there is such a force can be easily seen by making a Galilean transformation from the reference frame of the liquid to that of the impurity. There, the fermionic liquid flows past the impurity and the net force produced by scattering on the impurity has the same origin as the well-known residual resistance of metals. Moreover, this transformation allows a simple calculation of $`\eta `$ which can be expressed in terms of the impurity’s transport scattering cross-section.

The significance of Eq. (3) is that it relates an adiabatic effect — orthogonality catastrophe — with an essentially non-adiabatic effect — energy dissipation caused by the viscous force. In problems where the particle of interest can be singled out (like, for example, diffusion of heavy particles in Fermi liquids, the dynamics of topological defects, solitons, collective coordinates etc.) this connection provides a powerful tool that reduces the original strongly interacting problem to what is known as dissipative quantum mechanics of a single particle.  As for the Fermi liquid itself, one can not, in general, deduce the breakdown of Landau theory using these simple arguments. However, once the breakdown occurs, its direct cause — Anderson’s orthogonality catastrophe — can be linked to the dissipative dynamics of fermions. In the case of cuprate superconductors this idea is supported by experimental evidence of strongly damped quasiparticles.

Motivated by these connections, I study the problem of localization of a particle on attractive short-range potentials in $`D`$ dimensions in the presence of Ohmic dissipation at zero temperature. This problem is relevant to the localization of heavy particles on impurities in fermion liquids, pinning of vortices in superconductors  etc. It is also relevant to the problem of interaction between dissipative particles and to superconductivity. The question that I address is whether a short range attractive potential leads to localization of the particle and, if so, what are the properties of the localized state. I find that localization occurs for arbitrarily weak attractive potentials in any number of dimensions and for any strength of dissipation. This should be compared to the corresponding problem in the non-dissipative case  where it is known that weak potentials localize the particle in $`D=1`$, lead to exponentially weak localization in $`D=2`$ and do not lead to localization in $`D3`$. I suggest that the much stronger localization in the dissipative case may be relevant to the physics of Anderson localization and the problem of high critical temperature in cuprate superconductors.

Consider a particle of mass $`m`$ that moves in $`D`$ dimensions in the presence of an attractive short-range potential $`V(𝐱)=\alpha \delta (𝐱)`$. The particle also experiences the action of a friction force that is linear in the particle’s velocity: $`𝐅=\eta \dot{𝐱}`$ (Ohmic case). The well-accepted microscopic approach to this problem is to think of the friction force as originating from coupling to an effective environment. The environmental degrees of freedom are eventually traced over to produce an effective description of the particle’s dynamics. I use the imaginary (single) time path integral representation of the partition function.  This way, tracing over the environmental degrees of freedom reduces to path-integration over environmental variables, yielding the following effective action description:

$`Z`$ $`=`$ $`{\displaystyle 𝒟𝐱e^S},`$ (4)
$`S`$ $`=`$ $`{\displaystyle \frac{T}{2}}{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}\left(m\nu _n^2+\eta |\nu _n|\right)|𝐱_n|^2`$ (5)
$``$ $`\alpha {\displaystyle _0^{1/T}}𝑑\tau \delta (𝐱(\tau )).`$ (6)
Here $`\nu _n=2\pi Tn`$ are Matsubara frequencies: $`𝐱(\tau )=e^{i\nu _n\tau }𝐱_n`$.

The $`\eta |\nu _n|`$ term in the effective action $`S`$ is responsible for dissipation. Note that it has a non-analytic $`|\nu _n|`$ frequency dependence. This fact strongly complicates the problem because it leads to a non-local effective action when written in time representation. However, there are two limits in which the problem is somewhat simplified: the pure quantum limit and the dissipative limit. Indeed, the interplay between different terms in $`S`$ produces a characteristic energy scale $`\omega _0`$. Then, in the pure quantum limit $`\omega _0\eta /m`$ the dissipative term in $`S`$ is much smaller than the kinetic energy term and can, therefore, be neglected. The resulting pure quantum problem can be solved exactly.  In the opposite dissipative limit $`\omega _0\eta /m`$ the kinetic energy term can be neglected; its only purpose then is to set a high-energy cut-off.

To analyze the problem I employ the variational method. This method exploits the fact that for an arbitrary trial action $`S_0`$ the following inequality puts an upper bound on the free energy:

$`FF_0+TSS_0_{S_0}.`$ (7)
Here $`F=T\mathrm{ln}Z=T\mathrm{ln}𝒟𝐱e^S`$, $`F_0=T\mathrm{ln}Z_0=T\mathrm{ln}𝒟𝐱e^{S_0}`$ and $`SS_0_{S_0}=Z_0^1𝒟𝐱(SS_0)e^{S_0}`$. We are interested in the ground state energy: $`E=lim_{T0}F`$. Taking the low temperature limit in Eq. (7) we get:

$`EE_0+\underset{T0}{lim}TSS_0_{S_0}.`$ (8)
Here $`E_0=lim_{T0}F_0`$. I take for $`S_0`$ the action of an Ohmic-damped harmonic oscillator:

$`S_0={\displaystyle \frac{T}{2}}{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}\left(m\nu _n^2+\eta |\nu _n|+m\omega _0^2\right)|𝐱_n|^2.`$ (9)
Here, the $`\omega _0`$ of the oscillator will serve as a variational parameter. As we shall see below, all quantities that are necessary for the evaluation of Eq. (8) can be expressed through the expectation values of the operators that are diagonal in the $`x`$-representation. This means that we only need to know the diagonal elements of the equilibrium density matrix $`\rho _0`$ that corresponds to $`S_0`$. However, the fact $`S_0`$ is quadratic implies that $`\rho _0`$ is Gaussian. Thus, in the $`T0`$ limit we can write:

$`\rho _0(𝐱,𝐱)={\displaystyle \frac{1}{(2\pi x^2_0)^{\frac{D}{2}}}}\mathrm{exp}\left\{{\displaystyle \frac{𝐱^2}{2x^2_0}}\right\}.`$ (10)
Here, $`x^2_0=\frac{1}{D}𝐱^2_0`$ is the ground state expectation value of $`x^2`$ for the damped harmonic oscillator that corresponds to $`S_0`$. This quantity can be easily evaluated:

$`x^2_0=\underset{T0}{lim}T{\displaystyle \underset{n=\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{1}{m\nu _n^2+\eta |\nu _n|+m\omega _0^2}}.`$ (11)
In the $`T0`$ limit the sum can be substituted by an integral. Defining $`r=2m\omega _0/\eta `$, we can write the resulting expression as:

$`x^2_0={\displaystyle \frac{1}{\pi \eta }}\times \{\begin{array}{cc}\frac{1}{\sqrt{1r^2}}\mathrm{ln}\frac{1+\sqrt{1r^2}}{1\sqrt{1r^2}},\hfill & r<1\text{;}\hfill \\ \frac{2}{\sqrt{r^21}}\mathrm{arctan}\sqrt{r^21},\hfill & r>1\text{.}\hfill \end{array}`$ (12)
In particular, in the dissipative regime $`r1`$ we get: $`x^2_0=(2/\pi \eta )\mathrm{ln}(1/r)`$, while in the pure quantum regime $`r1`$ we obtain the correct quantum expression: $`x^2_0=1/2m\omega _0`$; at $`r=1`$ $`x^2_0`$ is analytic and equals $`2/\pi \eta =1/\pi m\omega _0`$. Note that $`x^2_0`$ for $`r>1`$ can be obtained from the expression for $`r<1`$ by an analytic continuation.

Next, using Eq. (10) we can evaluate:

$`\underset{T0}{lim}TSS_0_0`$ $`=`$ $`\alpha \delta (𝐱){\displaystyle \frac{m\omega _0^2}{2}}𝐱^2_0`$ (13)
$`=`$ $`{\displaystyle \frac{\alpha }{\left(2\pi x^2_0\right)^{\frac{D}{2}}}}{\displaystyle \frac{Dm\omega _0^2}{2}}x^2_0.`$ (14)
To evaluate $`E_0`$ we apply Hellmann-Feynman theorem: $`E_0/\omega _0=m\omega _0𝐱^2_0`$; this leads to:

$`E_0=Dm{\displaystyle _0^{\omega _0}}𝑑\omega _0^{}\omega _0^{}x^2_0.`$ (15)
Substituting Eqs. (13,15) into Eq. (8) and integrating by parts, we get a variational bound on the ground state energy:

$`E{\displaystyle \frac{Dm}{2}}{\displaystyle _0^{\omega _0}}𝑑\omega _0^{}\omega _0^2{\displaystyle \frac{}{\omega _0^{}}}x^2_0{\displaystyle \frac{\alpha }{\left(2\pi x^2_0\right)^{\frac{D}{2}}}}.`$ (16)
An immediate consequence of this result is the following statement: For any combination of non-zero values of dimensionality $`D`$, mass of the particle $`m`$ and friction constant $`\eta `$ the ground state energy $`E<0`$. And since for any delocalized state $`E0`$ this implies that the particle is always localized! To see this, notice that for small enough values of $`\omega _0`$ it follows from Eq. (12) that the first term in Eq. (16) is positive and $`\omega _0^2`$ while the second term is negative and its absolute value is $`(\mathrm{ln}\frac{1}{\omega _0})^{D/2}`$. Therefore for sufficiently small $`\omega _0`$ the second term dominates and the energy of the ground state is bounded from above by a negative value.

Minimizing Eq. (16) with respect to $`\omega _0`$ we obtain an equation on $`\omega _0`$:

$`{\displaystyle \frac{\alpha }{(2\pi )^{\frac{D}{2}}}}=m\omega _0^2x^2_0^{\frac{D}{2}+1}.`$ (17)
The upper bound on the ground state energy can be found numerically by solving Eqs. (16,17); this solution is presented below. However, in two limiting cases, $`\omega _0\eta /m`$ and $`\omega _0\eta /m`$, approximate analytical expressions can be obtained. These two limits are realized for weak and strong pinning potentials respectively. To be more precise, let us define $`\alpha _c`$ as the strength of the potential at which the minimizing value $`\omega _0`$ satisfies $`r=2m\omega _0/\eta =1`$. Substituting this in Eqs. (12,17) we obtain:

$`\alpha _c={\displaystyle \frac{2^{D1}}{\pi m}}\eta ^{1\frac{D}{2}}.`$ (18)
For weak attractive potentials, such that $`\alpha \alpha _c`$, we have the dissipative regime: $`\omega _0\eta /m`$. In this case, using the limiting form of Eq. (12) for $`r1`$ and solving Eqs. (16,17) we obtain:

$`E_{\mathrm{diss}}{\displaystyle \frac{\eta }{2\pi m}}\left({\displaystyle \frac{\alpha }{\alpha _c}}\right)\left(\mathrm{ln}{\displaystyle \frac{\alpha _c}{\alpha }}\right)^{\frac{D}{2}}.`$ (19)
This result shows that the localization on sufficiently weak attractive potentials is dominated by dissipative dynamics and is much stronger compared to the non-dissipative case. Indeed, note the linear dependence of $`E_{\mathrm{diss}}`$ on $`\alpha `$ in Eq. (19) vs. exact results in the pure quantum case: $`E_{\mathrm{quant}}=m\alpha ^2/2`$ in $`D=1`$; $`E_{\mathrm{quant}}=(m\lambda ^2)^1\mathrm{exp}(2\pi /m\alpha )`$ in $`D=2`$ ($`\lambda `$ is the short-distance cut-off); no localization on weak potentials in $`D3`$.

In the opposite limit of sufficiently strong potentials, such that $`\alpha \alpha _c`$, we have the pure quantum regime: $`\omega _0\eta /m`$. Here the situation depends on the dimensionality. In $`D=1`$, using the limiting form of Eq. (12) for $`r1`$ and solving Eqs. (16,17), we obtain:

$`E_{\mathrm{quant}}^{(D=1)}{\displaystyle \frac{m\alpha ^2}{\pi }}.`$ (20)
This bound is reasonably close to the exact quantum-mechanical result $`E=m\alpha ^2/2`$. In dimensions $`D2`$ the pure quantum problem of a strong ($`\alpha \alpha _c`$) attractive potential acquires a short distance cut-off (bandwidth) dependence; this is a direct consequence of the fact that $`\delta `$-function is poorly defined in $`D>1`$. However, it should be emphasized that the particle is always localized in any $`D`$ (see the arguments above)! While the question about the ground state energy in the strong potential case lies outside the main line of this letter, it could be mentioned that, in the first approximation, $`E(\alpha \alpha _c)`$ is bounded from above by $`E_{\mathrm{diss}}(\alpha \alpha _c)\eta /m`$ for any $`D`$.

Finally, let us discuss the numerical minimization of Eq. (16). The results of the solution for $`D=1`$ are presented in Fig. 1. A crossover between the $`\alpha ^2`$ behavior in the pure quantum regime (large $`\alpha `$) and the $`\alpha `$ behavior in the dissipative regime (small $`\alpha `$) occurs at $`\alpha \alpha _c`$. Note the much stronger localization that occurs for weak potentials compared to the pure quantum case. One can see that the two limiting expressions, Eqs. (19,20) describe the numerical results well in their respective limits.

To understand the qualitatively different localization in the dissipative case let us first recall how localization occurs in the non-dissipative case. There, it is determined by the competition of two tendencies of the particle: to localize and thus lower its potential energy and to diffuse away thus lowering its kinetic energy. In the dissipative case, however, the diffusion is strongly suppressed by the Anderson’s orthogonality catastrophe (see Eq. (2)). Indeed, as was shown by Hakim and Ambegaokar , the mean-squared displacement of a free Ohmic-damped particle at $`T=0`$ (this is known as quantum Brownian motion problem) grows very slowly with time:

$$x^2=\frac{2}{\pi \eta }\mathrm{ln}\frac{\eta t}{m}.$$
(21)
(This can also be shown  to be a direct consequence of the fluctuation-dissipation theorem of Callen and Welton ). Thus, strictly speaking, the diffusion coefficient of the Ohmic-damped particle is zero. Therefore, the tendency to lower potential energy by localization dominates over diffusion and leads to localization for any $`D`$, in contrast to the pure quantum case. Meanwhile, in one dimension dissipation leads to much stronger localization on weak potentials compared to the non-dissipative case.

Physical implications of these results are quite interesting. For example, if one considers the problem of diffusion of a heavy particle in fermionic liquids with impurities one arrives at the conclusion that the theory of Anderson localization breaks down due to the presence of dissipation (orthogonality catastrophe). In particular, no metal-insulator transition is possible for heavy particles in $`D=3`$ since in any dimensionality all states are localized.

Perhaps, a more important implication is for the interaction between damped particles. The well-known effect of attractive short-range interactions is superconductivity. In conventional superconductors superconductivity is known to arise from Cooper pair formation by electrons near the Fermi surface. For such electrons the density of states can be taken constant. Therefore, the problem of Cooper pairing becomes analogous (after a transformation to the center-of-mass and relative coordinates) to single-particle localization on attractive short-range potentials in $`D=2`$ where a single particle has constant density of states. In the latter problem, the binding energy is exponentially small:  $`EW\mathrm{exp}(1/\alpha \nu )`$, where $`W`$ is the band-width, $`\alpha `$ is the strength of the potential and $`\nu =m/2\pi `$ is the single-particle density of states in $`D=2`$. This explains why the critical temperature $`T_c`$ is so low in conventional superconductors. While no single accepted theory of superconductivity in the high-$`T_c`$ cuprate superconductors exists yet, most attempts to explain high $`T_c`$ begin by considering strong interactions: $`\alpha \nu 1`$, while still following the BCS ideas (extended to the strong-coupling case by Eliashberg). Here, I suggest that in systems where the Landau Fermi liquid theory breaks down and the dynamics of fermions is dissipative (governed by the orthogonality catastrophe) the nature of Cooper pairing changes qualitatively. Most importantly, the binding energy becomes much bigger and is essentially linear in $`\alpha `$ for weak attractive interactions (see Eq. (19)). This means that any weak second-order (fluctuation-induced) process is, in principle, capable of producing sufficiently high $`T_c`$.

It should be emphasized, however, that this model of superconductivity is a phenomenological one since the arguments given in the beginning of the letter can not be used to deduce the breakdown of the Landau theory. However, if the breakdown does occur, this model is believed to capture the qualitative features of the physical picture.

I would like to thank F. D. M. Haldane for sharing his insight into the problem. I have also benefited from discussions with L. P. Gor’kov and O. Motrunich. This work was supported by NSF MRSEC DMR-9809483.


text:
# Magnetic-history-dependent nanostructural and resistivity changes in Pr0.5Ca0.5Mn0.98Cr0.02O3

## Abstract

We show that nanostructure and resistivity of Pr0.5Ca0.5Mn0.98Cr0.02O3 are sensitive to whether the sample is zero-field-cooled (ZFC) or field-cooled (FC) either in the ‘self-magnetic field (H $``$ 2 T)’ of the electron microscope or under the external magnetic field of 2 T. FC resistivity at H = 2 T is lower than ZFC values below 140 K. The average value of the charge-orbital modulation vector (q = 0.44) of the FC crystallites is lower than that of the ZFC crystallites (q = 0.48) and the FC crystallites exhibit numerous defects like discommensuration, dislocations and regions with loss of superstructures compared to the ZFC crystallites.

PACS NO: 61.14.-x, 75.30.Vn, 75.60.Ee

The renaissance of interest in charge ordered manganites like Nd0.5Sr0.5MnO3 is due to insulator-metal transition driven by melting of charges under external magnetic fields. Charge ordering (CO) manifests itself as regular ordering of Mn3+:t$`{}_{}{}^{3}{}_{2g}{}^{}`$e$`{}_{}{}^{1}{}_{g}{}^{}`$ and Mn4+: t$`{}_{}{}^{3}{}_{2g}{}^{}`$ ions and ordering of d$`{}_{}{}^{2}{}_{z}{}^{}`$ orbital at Mn3+ sites. The low temperature phase of a perfect CO (Mn3+/Mn4+ = 1:1) system is a CE-type antiferromagnet, but some of the manganites exhibit phase seperation into CO (antiferromagnetic, AF) and non-CO (ferromagnetic, FM) regions. Electron microscopy is one of the powerful techniques to characterize the lattice distortion generated by CO and phase separation. Neverthless, sample in the electron microscope is under the influence of about 2 T (‘self-magnetic field’) and CO in some of the manganites is sensitive to this ‘self-magnetic field’ as found recently in Pr0.7Ca0.3MnO3. Coexistence of ordered regions (correlated to the CO-AF phase) and non-ordered regions (correlated to the FM phase) over wide x range ( 0.01 $``$ x $``$ 0.05) has been reported in Pr0.5Ca0.5Mn1-xCrxO3. In this letter, we report unusual magnetic history dependence of resistivity and nanostructural changes in Pr0.5Ca0.5Mn0.98Cr0.02O3.

Resistivity ($`\rho `$) and magnetization ($`M`$) of polycrystalline Pr0.5Ca0.5Mn0.98Cr0.02O3 were measured in two modes: In ZFC and FC modes, data were taken while warming from 5 K in presence of a known field (H) after cooling in zero-field and under H, respectively from 300 K. In order to obtain the ZFC condition in electron microscopy study, the objective lens of the microscope was cut off at 300 K before introducing the sample and the temperature was decreased to 92 K. After one hour at 92 K, the objective lens was switched on. Numerous crystallites have been characterized using electron diffraction (ED) in order to determine the average value of the charge modulation vector q and lattice imaging for the distribution of fringes. Different zones of different crystallites have been located. Then the temperature was raised to 300 K keeping objective lens switched on, again cooled to 92 K keeping exactly the same current density in order to obtain the FC mode. Exactly the same areas of the same crystallites were re-observed and the results are compared to those recorded under ZFC condition.

The temperature dependence of M for H = 0.01 T, 2 T and $`\rho `$ for H = 0 T and 2 T are shown in Figs. 1(a) and Fig. 1(b) respectively. We use solid curve for the ZFC and dashed curve for the FC modes. The peak in M around 230 K in both H = 0.01 T and H = 2 T signals ordering of charges and orbitals in the paramagnetic phase. A short range FM order is established below TC = 130 K as suggested by the difference in ZFC and FC magnetization at H = 0.01 T. $`\rho `$(0T) starts to decrease not at TC but around 80 K, close to the temperature where the ZFC magnetization also reaches a maximum. The decrease of $`\rho `$(T) below 80 K can be thus considered due to the percolation of FM clusters whose moments are blocked below this temperature. The peak in $`\rho `$(T) shifts to 125 K under H = 2 T and the magnetoresistance (\[$`\rho _{ZFC}`$(0 T)-$`\rho _{FC}`$(2 T)\] /$`\rho _{FC}`$(0T) ) at 100 K is as large as 73 %. The surprising result is the deviation of the ZFC and FC-$`\rho `$(T) below $``$ 135 K and lower values of the FC-$`\rho `$(T). The FC-$`\rho `$(T) exhibits a peak at Tp ( = 128 K) which is higher than that of ZFC-$`\rho `$(T) (Tp = 116 K). Under H = 2 T, TC increases to 140 K, the maximum in ZFC-M(T) shifts down $``$ 30 K and, FC-M(T) increases without saturation down to 5 K. This suggests that long-range FM order is not reached even at H = 2 T and canting of spins in the CO-AF phase contribute to non-saturation of M. We find that irreversibility in $`\rho `$(T) and M(T) presists even at H = 7 T in this compound but becomes negligibe with increasing Cr content  but, the details are beyond the scope of this letter.

Since $`\rho `$(T) values are lower when it is field–cooled, we can expect charge-ordering to be less stable in the FC state than in the ZFC state. The periodic lattice distortion due charge-orbital ordering generates satellites in Electron Diffraction (ED) patterns. The charge-orbital modulation vector exhibit a value q = 0.5-$`\epsilon `$ , where the incommensurability ($`\epsilon `$) is zero for a perfect 1 : 1 charge ordering. Fig. 2(a) and Fig. 2(b) respectively show the ED patterns at 92 K for the ZFC and FC modes. The average q value of the FC crystallites is close to 0.44 which is significantly lower than the average value q = 0.48 of the ZFC crystallites. This is clearly observed in the  ED patterns through the amplitude of the satellite splitting (see small black arrows). The satellites of the ZFC crystallites are also more intense than the ones observed in the FC crystallites. The corresponding lattice images are compared in Fig. 3(a) and 3(b). In the ZFC case (Fig. 3(a)), in addition to majority fringes of periodicity 2a$`{}_{p}{}^{}\sqrt{2}`$ = 10.8 Å , we also find a few fringes of higher periodicty with 3a$`{}_{p}{}^{}\sqrt{2}`$= 16.2 Å at irregular positions in agreement with average q = 0.48. However, in the FC case (Fig. 3(b), numerous defects are observed in the arrangements of fringes: discommensurations ( for example, area noted D) which are easily seen by viewing the images at grazing incidence and local losses of superstructure (for example, area within open arrows). To our best of knowledge, this is the first time influence of magnetic field on nanostructure has been reported. It is important to recall that artefacts due to bad stabilization of the temperature or hysteresis effect can generate the difference in the q values observed for exactly the same area. We repeated our experiments several times and all our results are identical. The satellites are stable under electron beam with long exposure contrary to what was found in Pr0.7Ca0.3MnO3. This suggests that CO is more stable in the present compound.

The above results unequivocally suggest that magnetic history of the sample affects not only the resistivity but also nanostructure. Why is it so ?. Recently, we have suggested that Cr3+ :t$`{}_{}{}^{3}{}_{2g}{}^{}`$ ions randomly occupy Mn3+:t$`{}_{}{}^{3}{}_{2g}{}^{}`$e$`{}_{}{}^{1}{}_{g}{}^{}`$ sites with their t$`{}_{}{}^{3}{}_{2g}{}^{}`$ spins opposite to the replaced Mn3+ ions. The random substitution of Cr3+ ions brings about eg-orbital deficiency at the original Mn3+-sites and affects cooperative Jahn-Teller distortion of Mn3+O6 octahedras. Hence, Cr3+ ions act as random impurities for spin as well as lattice degrees of freedom. Antiferromagnets with random impurities (‘diluted Ising antiferromagnets’) show long range AF order when they are zero field cooled but break into domains upon field cooling. While the sizes of the domains in theoretical models  are purely determined between volume and magnetic surface energies, spin-charge-lattice coupling also plays important role in manganites. Hence, the long range d$`_{z^2}`$ orbital ordering is also broken by the formation of AF spin domains which also causes orbital domains to form. Such orbital domains lead to discommensuration and dislocations observed experimentally. Field cooling is also found to enhance the size of FM domains which are clearly visible in lattice images as regions without superstructure. The breaking of CO-AF matrix into domains and growth of FM domains lead to low resistivity under field cooling.

In conclusion, we have shown that nanostructure and resistivity of Pr0.5Ca0.5Mn0.98 Cr0.02O3 are strongly magnetic field history dependent. Our study clearly suggests that electron microscopy can be used to tune nanostructure in phase separated manganites.

Acknowledgements:

R. M thanks MNENRT (France) for financial support and acknowledges Prof. T. V. Ramakrishnan for his encouragement.

FIGURE CAPTIONS

: (a) Temperature dependence of magnetization at H = 0.01 T and 2 T in ZFC (solid curve) and FC (dashed curve). (b). $`\rho `$(T) under ZFC (solid curve) and FC (dashed curve) for H = 2 T and in absence of field (H = 0 T) modes. TCO: Charge ordering temperature, TC: ferromagnetic Curie temperature.

:  ED patterns at 92 K underZFC (a) and FC (b) modes. Note that satellites marked by double arrows are clearly splitted and less intense in FC mode than in ZFC mode.

:  lattice images at 92 K under ZFC (a) and FC (b) modes. Regions with loss of superstructures are denoted by open arrows and one region with discommensuration is marked by D.


text:
# Computer Science in Physics

## 1 Introduction

An important part of the area of physics known as “statistical physics” is the study of phase transitions, at which the system converts from one state to another. Most interest has centered on “second order” or “continuous” transitions, in which the property which distinguishes the two phases vanishes continuously as the transition is approached. The disappearance of the magnetization of a ferromagnet, such as iron, as the temperature is increased is generally continuous. At the other type of transition, known as “first order” or “discontinuous”, there is a jump in the properties of the system as the transition is crossed, and also a latent heat. An everyday example of a first order transition is the freezing of water.

We shall focus on magnetic transitions in this talk because (i) they can be represented by simple models amenable to numerical study, and (ii) there are many experimental systems which are describable by these models. One of the major advances in the field has been the realization that behavior in the vicinity of the transition (called the “critical point”) is “universal” . This means that “critical behavior” does not depend on the microscopic details of the system but only on much more basic features such its dimensionality and symmetry. Consequently one can use relatively simple models, which can be readily simulated, to make precise comparisons with experiment. That said, it should be emphasized that universality is much better justified for clean systems than for the systems with disorder which we shall be considering in this talk. One goal of applying sophisticated algorithms from computer science to these problems will be to see if universality also holds for disordered systems.

The simplest model which describes a magnetic transition, known as the Ising model, has a variable at each site on a regular lattice which can point either “up” or “down”. This represents the orientation of the magnetic moment of an atom, and is simplified to only allow two possible orientations. We shall follow standard notation in calling these variables “spins”, labeled $`S_i`$ where $`i`$ denotes a lattice site. It is convenient to denote the up spin state by $`S_i=1`$ and the down spin state by $`S_i=1`$.

Neighboring sites on the lattice interact with each other. If the interaction favors parallel alignment of the spins, it is called “ferromagnetic”, while an interaction favoring anti-parallel alignment is called “anti-ferromagnetic”. The energy (confusingly called the “Hamiltonian” in the physics literature) can therefore be written as

$$E=\underset{i,j}{}J_{ij}S_iS_j,$$
(1)
where the sum is over all nearest neighbor pairs of the lattice (counted once) and the interactions are labeled $`J_{ij}`$. If all the $`J_{ij}`$ are positive then the state of lowest energy (ground state) has all spins parallel (either all $`+1`$ or all $`1`$) and is called a ferromagnet. As the temperature is increased, the net alignment of the spins, the magnetization, decreases and vanishes at a critical temperature $`T_c`$, where thermal noise, which tends to randomize the spins, overcomes the interactions, which tend to make them order.

The situation with all interactions negative is simple if the sites of the lattice can be divided into two sublattices, A and B, such that all the neighbors of A are in B and vice-versa. Such lattices are said to be “bi-partite”, and are the only type that we shall consider here. A square grid is an example of a bipartite lattice. The ground state for a bipartite lattice with negative interactions has all spins $`+1`$ on sublattice A and all spins $`1`$ on sublattice B, or vice-versa. Such a state is called an antiferromagnet. Again the ordering decreases to zero at a critical temperature.

The problems of interest to us have two more ingredients. The first is disorder. The interactions are not all equal but are chosen in some random way. The simplest model for disorder is to pick each interaction from a probability distribution, independently of all the others. The second ingredient is “frustration” or competition between the interactions. For the model in Eq. (1) this can be incorporated by allowing the sign as well as the magnitude of the interaction to be random.

That this leads to frustration can be seen in the “toy” example in Fig. 1. This shows just four sites round a square with one anti-ferromagnetic (negative) interaction and three ferromagnetic (positive) interactions. If the spins along the bottom row and top left corner are oriented in directions which minimize the energy (as shown) the spin at the top right is “frustrated” since it receives conflicting instructions from its two neighbors. It wants to be parallel to both of them which is impossible. It is easy to see that there is frustration if there is an odd number of negative interactions round the square, which is then called a “frustrated square”.

If we extend this toy example of four sites to a large lattice, choosing the sign (and possibly the magnitude) of the interactions at random, determining the ground state is non-trivial, which is not the case if all interactions have the same sign. This type of problem, which has been extensively studied, is called a “spin glass” . There are many experiments on magnetic systems with these features of disorder and frustration but it would take us too far afield to discuss them here. Spin glasses are also considered as prototypes for other systems with frustration and disorder which have many features in common. Examples are neural networks, protein folding models, elastic manifolds in random media, and the “vortex glass” transition in superconductors in a magnetic field. These are discussed in the articles in Ref.

Determination of the ground state in systems with disorder and frustration is an optimization problem, in which the “cost function” that has to be minimized is the energy. As we shall see, algorithms from computer science enable us to calculate the ground state of spin glasses for surprisingly large lattice sizes, at least in certain cases. An excellent introduction to optimization algorithms as applied to problems in physics is the review by Rieger .

Spin glasses and other problems with disorder and frustration are hard because the energy varies in a complicated manner as one moves through configuration space. There are local minima of the energy, which we will call “valleys” in the “energy landscape”, separated by “barriers” (i.e. saddle points). Different local minima can have similar energies but have very different configurations of the spins. At a finite temperature the systems should spend time in different valleys with relative proportions given by the appropriate Boltzmann factors . Only one of the local minima will be the global minimum (ground state). This can be hard to find if there are many minima and/or the global minimum has a small “basin of attraction”. However, it is generally quite easy to find a minimum with energy close to the ground state energy, for example by the method of “simulated annealing” .

The precise value of the ground state energy will depend on the particular choice of the random interactions (remember they were picked from a distribution). In physics we usually look at “intensive” quantities (those which do not depend on size of the system, $`N`$, as $`N\mathrm{}`$) such as the ground state energy per spin. Many intensive quantities are “self-averaging” which means that its value does not depend on the realization of disorder for $`N\mathrm{}`$. However, there are sample-to-sample fluctuations, generally of order $`1/\sqrt{N}`$, for finite-sized systems, so we need to average results over many realizations of disorder. This makes the problem even more computationally challenging than if we just had to solve for one sample, but fortunately averaging over samples is clearly “trivially parallelizable”, so we can easily take advantage of the large-scale parallel machines that are widely available at present, or just run the code on a “farm” of independent workstations.

In this talk I will also discuss another widely studied problem with frustration and disorder, known as the random field model. A magnetic field will prefer a spin to align in one direction rather than the other, and so can be represented in the expression for the energy by terms linear in the spins. Eq. (1) is therefore modified to

$$E=\underset{i,j}{}J_{ij}S_iS_j\underset{i}{}h_iS_i,$$
(2)
where we have allowed for a different field $`h_i`$ on each site. The random field model is obtained if one chooses the $`h_i`$ at random with zero average value, and has the $`J_{ij}`$ unfrustrated (so we could set them all to equal unity). Again there are experimental systems which have been widely studied but which space does not allow me to discuss. For more information see Ref.  and the articles in Ref. . In Eq. (2) frustration comes from competition between the interactions on the one hand, which prefer the spins to be parallel, and the random fields on the other, which prefer the spins to follow the local field direction.

The traditional physics approach to studying problems with frustration and disorder is the Monte Carlo simulation, i.e. a random sampling of the states according to the Boltzmann distribution . However, at low temperatures the system gets trapped in one of the valleys for a long “time” and is only very rarely able to escape over a barrier in the energy landscape to another valley. The probability for escape is exponentially small in the ratio of the barrier height to the temperature. As a result, equilibrium simulations can only be done on very small systems at low temperatures. Some speed up can be obtained from recently developed Monte Carlo algorithms such as parallel tempering  (also known as “exchange Monte Carlo”) but the range of sizes that can be studied is still quite limited.

In this talk I will discuss an alternative approach which uses sophisticated optimization algorithms from Computer Science  to find the exact ground state. The idea will be to “beat the small size limit” of Monte Carlo methods. The advantages of the computer science approach are:

1. It is exact. There are no statistical errors or problems of equilibration.
2. One can study large sizes.

However, there are also some disadvantages. These are:

1. Only the ground state is determined, so one is restricted to zero temperature properties.
2. Only for some models are there efficient algorithms.

In the rest of this talk I will discuss what has been learned from applying optimization algorithms to the spin glass and random field problems, and also describe some prospects for the future.

## 2 The Random Field Model

In this section I will discuss how optimization algorithms have enhanced our understanding of the random field model. The energy is given by Eq. (2) with the sites on a regular lattice, which we take to be a square grid in two-dimensions, a simple cubic grid in three-dimensions, and similarly in higher dimensions. The interactions $`J_{ij}`$ are all set to unity and the random fields are chosen from a symmetric distribution with mean and variance given by

$$[h_i]=0,[h_i^2]=h_R^2,$$
(3)
where the rectangular brackets $`[\mathrm{}]`$ denote an average over the disorder, so $`h_R`$ is the strength of the random field.

In the absence of random fields it is known that there is a non-zero magnetization $`S_i`$ at low temperatures and we say that there is “long range order”. This long range order vanishes continuously at critical temperature. When the random fields are turned on one could ask whether even a small random field prevents the formation of long range order at any temperature or whether a critical field strength is needed to destroy long range order at low temperature. A famous argument due to Imry and Ma  states that for dimension two and lower, the random field always “wins” in the sense that long range order is destroyed by an arbitrarily weak random field, with the system “breaking up” into domains of parallel spins. The domain size diverges as $`h_R0`$ so one recovers long range order for $`h_R`$ strictly zero. However, in dimension, $`d`$, greater than two an arbitrarily small random field does not cause the system to break up into domains on long length scales and long range order is preserved up to a critical field strength.

For $`d>2`$, the phase diagram is that sketched in Fig. 2. For $`h_R=0`$ the ferromagnetic phase disappears at $`T=T_c`$ due to thermal fluctuations, while at $`T=0`$ the ferromagnetic phase disappears at a critical value of the random field, $`h_c`$, due to the disordering effects of the random field. This will be important later.

What aspects of the random field problem are physicists interested in? It turns out that many quantities vary with a power law in the vicinity of the critical point. Denoting by $`\delta `$ the deviation from the phase boundary in Fig. 2, the magnetization varies, for $`\delta `$ small, like

$$S_i|\delta |^\beta ,$$
(4)
where $`\beta `$ is known as a “critical exponent”. Other quantities of interest are the specific heat, $`C`$, and the magnetic susceptibility, $`\chi `$, which vary as

$`C`$ $``$ $`\delta ^\alpha ,`$
$`\chi `$ $``$ $`\delta ^\gamma .`$ (5)
Critical exponents such as $`\alpha ,\beta `$ and $`\gamma `$ are of interest because they are universal , only depending on broad features of the problem such as the dimensionality of the lattice and whether or not there is a random field. They are not expected to depend on the strength of the random field, as long as it is non-zero, or the distribution used for the random fields, as long as it is symmetric.

Physicists would like to the understand universal critical behavior, such as the values of the critical exponents. Since the phase boundary is crossed at zero temperature by varying $`h_R`$ one can investigate the behavior near the phase boundary using optimization (i.e. $`T=0`$) algorithms. Fortunately, determining the ground state energy is equivalent to a Max-Flow problem  which can be solved in polynomial time. This was first realized by Barahona and Anglès d’Auriac et al.  and subsequently used by Ogielski , Sourlas and collaborators , and others.

Here I will just discuss briefly some of the results found in Ref. , which investigated the random field model in $`d=3`$. In their optimized implementation of the Max-Flow algorithm, they studied $`L^3`$ lattices up to $`L=90`$, and found empirically that the CPU time varied as $`L^4`$. This is remarkably efficient, being not much more than the time ($`L^3`$) needed to scan once through the lattice. Ref.  provides strong evidence that the transition is actually discontinuous, corresponding to an exponent $`\beta =0`$. This had been suspected earlier from finite-$`T`$ Monte Carlo simulations  on sizes up to $`L=16`$ but the results of Ref.  are more convincing because they are on much larger systems. Normally a first order transition leads to a latent heat and rather weak fluctuation effects compared with a continuous transition. However, no latent heat is seen in the random field problem and, in other respects, there seem to be large fluctuation effects characteristic of a continuous transition. This dichotomy is not understood.

Several different types of random field distribution were used in Ref. . While they all gave rise to $`\beta =0`$, other quantities, also expected to be universal, seemed to depend on the type of disorder, casting some doubt on the hypothesis that universality holds for random systems. This important question needs further work.

The $`T=0`$ approach cannot easily determine the specific heat exponent $`\alpha `$, which is unfortunate because there is a discrepancy between the experimental value  which is close to zero and results from Monte Carlo simulations, e.g. Ref. , which give $`\alpha 0.5`$.

## 3 The Spin Glass

The energy of the spin glass problem is given by Eq. (1) where the $`J_{ij}`$ are taken from a symmetric distribution with mean and variance given by

$$[J_{ij}]=0,[J_{ij}^2]=1.$$
(6)
One often takes a Gaussian distribution, though another popular choice is is the bimodal distribution, also the called $`\pm J`$ distribution, in which the interactions have values $`+1`$ and $`1`$ with equal probability. The latter distribution has the special feature that there are many ground states (we say that the ground state is “degenerate”). In fact the number of ground states is exponentially large in the number of spins $`N`$ giving rise to a finite ground state entropy.

Two principal questions have been asked about spin glasses:

1. Is there a phase transition at finite temperature $`T_c`$?
2. What is the nature of the spin glass state below $`T_c`$?

For the first question, Monte Carlo simulations and early (unsophisticated) ground state calculations have shown;

* In $`d=2`$ the transition is only at $`T=0`$.
* In $`d=3`$ (and higher) the transition is at finite temperature.

The conclusion for $`d=2`$ is very strong and so is the situation in $`d=4,5,\mathrm{}`$ etc. The case of $`d=3`$ has been the most difficult to resolve and earlier work was not very conclusive, but the most recent simulations  seem rather convincing.

Concerning the second question, we have already noted that, because of the complicated energy landscape, there are large clusters of (carefully chosen) spins which can be flipped with rather low energy cost. Is it possible to quantify this remark? Two principal scenarios have been proposed which differ, mainly, as to the energy of these large-scale excitations. These scenarios are:

* The “droplet model” of Fisher and Huse . In this phenomenological picture a few very plausible assumptions are made. The lowest energy to create an excitation of linear size $`L`$ is assumed to vary as

$$\mathrm{\Delta }EL^\theta ,$$
(7)
where $`\theta (>0)`$ is an exponent. $`\theta `$ can not be negative if $`T_c>0`$ otherwise there would be large scale excitations which cost vanishingly small energy and the system would break up into domains at any finite temperature. We shall see below that this is what happens in $`d=2`$. Note that for a ferromagnet there is a positive energy cost for each interaction on the wall of the excitation and, since the wall area goes like $`L^{d1}`$, one has $`\theta =d1`$ in that case. However, for a spin glass it turns out that $`\theta <d1`$ (in fact it is also true that $`\theta <(d1)/2`$). Hence, there is a near cancellation between the effects of the bonds which were “unsatisfied” before the excitation is flipped and then become “satisfied” (which lower the excitation energy), and the the effects of the satisfied bonds which become unsatisfied (which increase the excitation energy).
* The “replica symmetry breaking” (RSB) picture of Parisi . The Parisi theory is the (presumably) exact solution of an artificial model with infinite-range interactions. The assumption is then made that qualitatively similar behavior also occurs for more realistic models with short range interactions. An important ingredient of the RSB picture is that there are excitations of order the size of the system whose energy does not grow with the size of the system i.e.

$$\mathrm{\Delta }E\mathrm{const}.$$
(8)
This is in contrast to the prediction of the droplet theory in Eq. (7). The cancellation between the effects of the satisfied and unsatisfied bonds on the boundary of the excitation is then even more complete than in the droplet model.

To discuss what has been learned from optimization algorithms it is necessary to distinguish $`d=2`$ from higher dimensions. We first consider $`d=2`$.

We have already noted that a square is frustrated if an odd number of its bonds are negative and the converse, that the square is unfrustrated (i.e. each bond can be satisfied) if there are an even number of negative bonds, is also true. Changing the sign of the bonds in such a way that the frustration remains unchanged has no effect on the ground state energy because it can be compensated for by changing the sign of appropriate spins. Hence, for the $`\pm J`$ distribution, the ground state energy is determined entirely by the location of the frustrated squares. For a distribution in which the magnitude of the bonds is not constant we also need to keep track of the magnitude of the bonds (though not the sign) plus the location of the frustrated squares.

Let us therefore indicate the frustrated squares on the lattice by drawing a cross in their center, as shown in Fig. 3. We indicate the unsatisfied bonds by drawing dashed lines at right angles across them. Clearly the lines must begin and end on frustrated squares and so form “strings” connecting the crosses. The ground state energy is increased relative to that of the unfrustrated system by the energy of the bonds crossed by the strings. For a $`\pm J`$ distribution the ground state energy is therefore determined by minimizing the length of the strings. For a distribution of interactions where the magnitude as well as the sign varies one has to minimize the total “weight” of the string, where the weight of a string segment is equal to the magnitude of the bond which it crosses.

This problem is equivalent to a Minimum Weight Perfect Matching Problem  as first realized by Barahona et al.  which can be solved in polynomial time, i.e. it belongs to the category “P” of optimization problems. To be more precise it is a polynomial algorithm provided the lattice is a “planar graph”, i.e. it can be drawn on a piece of paper with no lines crossing. Unfortunately this rules out periodic boundary conditions which are often imposed to eliminate surface effects which arise from the spins on the surface having a different number of neighbors from the spins in the bulk. With periodic boundary conditions the problem belongs to the class “NP”. However, an efficient “Branch and Cut” algorithm  enables quite large sizes to be studied .

In three dimensions or higher calculating the ground state of spin glass is NP for all boundary conditions. Most work has used “heuristic” algorithms, which are not guaranteed to give the exact ground states, but which, when used carefully, do seem to give the true ground state in most instances. The most effective such approach seems to be the “genetic algorithm” developed for spin glasses by Pal  and subsequently used by Palassini and me  and Marinari and Parisi .

Let us now discuss what has been learned about spin glasses from optimization techniques, starting with $`d=2`$.

First we note that the restriction to $`T=0`$ is more serious than for the random field problem because we cannot go through the critical point. We can, however, learn about low energy excitations by computing the ground state, then perturbing the system in some way, and finally recomputing the ground state. As an example let us start with periodic boundary conditions and then change to anti-periodic boundary conditions in one direction, which simply corresponds to changing the sign of the bonds across one boundary. This induces a domain wall across the system, as shown in Fig. 4, such that all the spins on one side of the wall are flipped.

The wall will have an energy, which could have either sign, and its characteristic scale varies with the size of the system in each direction, $`L`$, as in Eq. (7). Starting with the pioneering work of Bray and Moore , and followed by later studies  for larger sizes (up to about $`30^2`$) using the Branch and Cut method, it has been found that $`\theta `$ is negative in $`d=2`$ with a value of around $`0.28`$. The negative value means that the system will break up into large domains at any finite temperature, so there is no spin glass state except at $`T=0`$. These studies also show that the wall is a fractal with fractal dimension $`d_s`$ about 1.28, greater than 1 so it is not a smooth curve, but also less than 2 so the wall is not “space filling”.

Recently Middleton  has used the Matching algorithm to determine ground states of two-dimensional spin glasses with free boundaries for very large sizes, up to $`512^2`$, using a different approach to generate excitations from which $`\theta `$ and $`d_s`$ can be determined. His results agree well with the other work.

The case of $`d=3`$ is more interesting that that of $`d=2`$ not only because this is the physical dimension but also because there is a finite temperature spin glass state. There is general agreement that $`\theta 0.20`$ in $`d=3`$, starting with the first studies of Bray and Moore  which could only consider sizes up to $`4^3`$, and followed by later work  which could go up to about $`10^3`$ using heuristic optimization algorithms. The positive value indicates that the spin glass state should be stable at low but finite temperatures. Ref.  also find that $`d_s2.68`$.

Subsequently, two papers  have argued that a different value of $`\theta `$, consistent with zero, is obtained from excitations in which the boundary conditions are not changed but a carefully chosen set of spins is flipped, for example by thermal noise. This suggests that the spin glass is actually quite close to the RSB picture. However, the sizes are still quite small, up to around $`10^3`$, and the assertion that there are two values (at least) for $`\theta `$ depending the type of excitation being considered is quite messy, so this claim needs further study.

Given the considerable interest in the spin glass in three-dimensions, it is unfortunate that there are no polynomial algorithms for finding the ground state. It is to be hoped that, in the future, algorithms will be developed which both give exact ground states and can treat larger sizes than the present heuristic algorithms.

## 4 Conclusions

This discussion of the role of optimization algorithms in statistical physics has been very brief. For further information the reader should consult the references. Ref.  is a good place to start.

The main conclusions of this talk are:

1. Algorithms from computer science have “broken the size barrier” for some problems in statistical physics, e.g.

1. The random field model
    2. The spin glass in two-dimensions
2. The application of algorithms from computer science to physics problems works best as a collaboration between computer scientists and physicist, e.g. Ref. .
3. For the future I expect there will be developments in the following areas:

1. More models will be solved.
    2. More efficient algorithms will be developed for NP problems such as the spin glass in three-dimensions. So far, with the genetic algorithm, we can study up to of order $`10^3`$ spins. Surely we do better than this.
    3. Spin glasses will be used to investigate the statistics of “hardness”. For example, given an algorithm for the exact ground state of a three-dimensional spin glass such as branch-and-cut, one can study the distribution of CPU times required to solve the problem for different realizations of disorder. It would be interesting to see how the average CPU time varies with system size and do the same for the typical (e.g. median) CPU time. If the distribution of CPU times is very broad, the average may be dominated by a few rare samples which are extremely “hard” and vary with size in a different way from the typical CPU time. This distinction has been made recently in statistical physics in the study of some quantum systems undergoing phase transitions at zero temperature  but, to my knowledge, does not seem to have been investigated systematically in studies of hardness of NP problems.

Acknowledgments:
I would like to thank Heiko Rieger for educating me about the algorithms mentioned in this talk and many other topics. Much of my own work in this field has been with Matteo Palassini and I would like to thank him for a stimulating collaboration. I am especially grateful to the organizer, Reinhard Wilhelm, for inviting me to speak at the 10th Anniversary Dagstuhl conference, which introduced me to large areas of computer science about which I knew nothing before. My research is supported by the NSF under grant DMR-9713977.