text:
# Point defect energetics in silicon using the LDA+U method

## I Introduction

Modeling techniques empower the investigator with the capability to explore into the details of natural phenomena with great dexterity, aptly complementing their experimental counterparts. Hierarchical multiscale modeling (HMM) is one such technique. We have recently investigated self diffusion in silicon-germanium alloys in considerable detail by developing a database of first principles energetics and using the database to perform kinetic Monte Carlo simulations: a typical HMM scheme Ramanarayanan et al. (2003). The exponential growth in computational processing power makes such schemes more viable; the explosive growth in nanotechnology provides technologically relevant platforms that are amenable to such schemes.

Understandably, the accuracy of the results such as those presented in Ref. Ramanarayanan et al., 2003 hinges to a great extent on the correctness of the energetics database. In that study, we used the popular local density approximation (LDA) to develop the energetics database. However, the (approximately 1 eV) discrepancy Souza and Narayanan (1998) between the theoretical activation energy for Si self diffusion computed using the LDA (or the generalized gradient approximation (GGA)) and the experimental values was an important reason that precluded us from making direct comparisons of our results with experimental studies of self diffusion in SiGe alloys such as those in Refs. Zangenberg et al., 2001; Strohm et al., 2002.

Our present work is directed essentially at finding a means of developing a reliable point defect energetics database (like in Ref. Ramanarayanan et al., 2003) using an as parameter free a technique as possible. (While considerable progress has been made on techniques like the quantum Monte Carlo (QMC) method, it is beyond the current computational capacity to use it to develop a detailed energetics database as required for our HMM scheme.) We focus on the activation energy for self-diffusion in pure Si. For self diffusion in Si, we consider the vacancy mechanism, the interstitialcy mechanism, and the concerted exchange mechanism. These three have been shown to be the most likely diffusion mechanisms in silicon Ural et al. (1998, 1999a). For the interstitialcy mechanism, we consider the hexagonal to split $`110`$ mechanism as this has shown to be the most likely mechanism Needs (1999).

This article is organized as follows: In Sec. II, we provide the details of the LDA and the LDA+U calculations. In Sec. III, we present the results of our calculations. We provide an intuitively appealing explanation for the underperformance of the LDA. We then provide a brief overview of the LDA+U technique and discuss the results obtained using this method for the different diffusion mechanisms. We summarize our article in Sec. IV.

## II Method

Self-consistent electronic structure calculations were performed using the plane-wave pseudopotential code VASP Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) with the projector augmented-wave (PAW) potentials Kresse and Joubert (1999); Blöchl (1994) at a kinetic energy cutoff of 20 Rydberg. A 64 atom supercell was used. Electronic minimization was carried out to a tolerance of $`2.7\times 10^5`$ eV and structures were relaxed until the maximum force on any atom was less than 0.015 eV/Å. Saddle point configuration for the hexagonal to split $`110`$ interstitialcy mechanism was determined using the nudged elastic band method Jónsson et al. (1998). Structural relaxation was performed in two stages: initially with a $`2^3`$ Monkhorst-Pack Monkhorst and Pack (1976) k-point sampling followed by a $`6^3`$ Monkhorst-Pack Monkhorst and Pack (1976) k-point sampling. Energy calculations were done using tetrahedron method with a $`6\times 6\times 6`$ division of Brillouin zone. The lattice constant of systems containing point defects (vacancy or interstitial) was determined by fitting the total energy versus the supercell volume to Murnaghan’s equation of state111The Murnaghan’s equation of state Murnaghan (1944) is a relation between the energy $`E`$ and the volume $`V`$ of the supercell and has the following form: $`E=E_0+[(VV_0)/b]\{[(V^{1a}V_0^{1a})V_0^a]/[(1a)b]\}`$. $`E_0`$, $`V_0`$, $`a`$, and $`b`$ are the fitting parameters..

The plane-wave psedopotential code VASP Kresse and Hafner (1993, 1994); Kresse and Furthmüller (1996a, b) also has the option to perform calculations by the LDA+U method. (A brief overview of the method itself will be presented in Sec. III. Here we merely give the details of the parameters used in our calculations.) We have used the rotationally invariant LDA+U scheme according to Liechtenstein et al. Liechtenstein et al. (1995). The on-site interactions were added for the p-orbitals. The effective on site Coulomb interaction parameter ($`U`$) was set to 0 eV and the effective on site exchange interaction parameter ($`J`$) was set to 4 eV. (We provide justification for this choice in Sec. III.) The LDA was used to perform structural relaxation; the LDA+U was used only to perform energetics calculation on the structure obtained using the LDA. We adopted this approach222We would like to point out that Leung et. al. Leung et al. (1999) take a similar approach: they calculate their structure using LDA and use QMC to compute the energetics. because the LDA gives structutal properties (eg. lattice constant, see Table 1) that are closer to experimental values than the LDA+U.

## III Results and Discussion

Table 1 summarizes the results of both the LDA and LDA+U calculations. It also contains experimental values and values from quantum Monte Carlo simulations where available.

There is consensus on the values obtained within the LDA by different theoretical groups Probert and Payne (2003); Puska et al. (1998); Needs (1999). It can be seen that the LDA gives activation energy that is systematically lower than the experimental values.

### III.1  Local density approximation

One of the main reasons for the underperformance of the LDA is that it tends to overbind. Because of this, the dangling bonds that are created due to the defects, tend to bind with each other and with other orbitals. Because of this overbinding, the energy of the system with the defect is lower than it should be in reality. Thus the difference between the energy of the pure crystal and the one with the defects (which is in fact the formation or the migration energy, depending on the situation) is lower than it should be in reality. The reason for the LDA’s overbinding can be traced to the spherical nature of the exchange hole. This leads to a larger negative value of the exchange energy and hence lowers the energy of the system much more than it should be in reality Koch and Holthausen (2000).

A related reason for the overbinding nature is the well-known underestimation of the band gap by LDA. From our LDA based calculations, we find that, for example, the Si with a vacancy has a band gap of only 0.88 eV in contrast to experimental values of 1.12 eV. (We note that the band gap in pure Si from LDA calculations is even less: 0.56 eV (Figs. 1, 2). The system with the impurity has a higher band gap because of the interaction of the defect with the other atoms in the lattice. This can be seen from a simple tight binding analysis. We consider the band gap of the system with the defect because that is the environment that the defect “perceives”.) Because of the small band gap, the defect states that are (usually) created in the band gap tend to interact more strongly with the valence states. This therefore leads to overbinding too. The band structure and the density of states of the system with a vacancy, hexagonal interstitial and split $`110`$ interstitial are shown respectively in Figs. 3, 4, and 5.

There are a couple of interrelated reasons for the LDA’s underestimation of the band gap. The discontinuity of the one-electron potential for localized states, which is a characterisitc of the exact density functional, is absent in the LDA as was shown by Perdew et. al. Perdew et al. (1982). Because the band gap can be expressed as

$`E_{gap}=E[N+1]+E[N1]2E[N]`$ (1)
where $`E[x]`$ is the energy of the system with $`x`$ electrons, the absence of the discontinuity of one-electron potential in the LDA causes the LDA to give an incorrect band gap. The second reason for the underestimation of band gap within LDA is the absence of self interaction correction. The third reason has to do with the Kohn-Sham approach Kohn and Sham (1965) itself. In the Kohn-Sham theory, there is no direct relationship between the orbital energies and the ionization energies (unlike that in Hartree-Fock theory) except for the highest occupied orbital.

### III.2 Local density approximation with Hubbard type correction: LDA+U

The LDA+U method has been widely applied to metallic systems containing localized electrons. We refer the reader to the review article on LDA+U by Anisimov et. al. Anisimov et al. (1997). Here we provide a brief explanation of the technique along the lines of “rectifying” the deficiencies of the local density approximation that were pointed out in Sec. III.1 . The basic idea of the LDA+U method is to redistribute the electron-electron interaction. Because, as we explained in Sec. III.1, overbinding of the localized electrons associated with the dangling bonds is a primary reason for the underperformance of the LDA, a direct way of rectifying this deficiency would be to add an on-site Coulomb repulsion term while subtracting an average overall electron-electron interaction.

$`H_{LDA+U}=H_{LDA}+{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1;ji}{\overset{N}{}}}U_{ij}{\displaystyle |\psi _i(𝒓_1)|^2𝑑𝒓_1|\psi _j(𝒓_2)|^2𝑑𝒓_2}U{\displaystyle \frac{N(N1)}{2}}`$ (2)
where $`N=_{i=1}^N|\psi _i(𝒓)|^2𝑑𝒓`$ corresponds to the total number of electrons in the system. The parameter $`U`$ is adjusted by trial and error so that the resulting band gap equals the experimentally observed band gap of the system in question. The above technique, which only has the Coulomb term, is generalized to include the exchange term as well. We refer the reader to Refs. Liechtenstein et al., 1995; Solovyev et al., 1994 for more details. The $`UN(N1)/2`$ introduces the discontinuity of the of the one-electron potential. We wish to point out that the LDA+U does not correct for the self interaction.

As mentioned in Sec. II, we have used the scheme as implemented in VASP for our calculations. Ideally, the Coulomb ($`U`$) and exchange ($`J`$) interaction parameters would be chosen so that the band gap for the bulk system would be equal to the experimentally observed band gap. However, because our calculations are based on a supercell geometry, there is an inevitable but unphysical interaction between the defect and the matrix. This alters the band gap from that in the bulk. However, because the system “sees” only this band gap, we have chosen $`U`$ and $`J`$ so that the band gap in the system with the defect is close to the experimentally observed band gap. This, of course, would mean that the parameter be tuned for each type of defect. However, for the sake of demonstration, we have worked with a single set of parameters.

We find that the defect formation and migration energies increase almost linearly with the difference: $`UJ`$. The best agreement between the above energies and the experimental data occurs when the calculated band gap corresponds to the experimental one. (This occurs when we choose $`U=0`$ and $`J=4eV`$ as mentioned in Sec. II.) This supports our argument that the correction of the density functional for the bandgap (discontinuity of the potential on the number of electrons) should improve the results. Table 1 summarizes our energetics calculations. Figs. 1-5 show band structures and densities of states of the various systems considered.

### III.3 Defect structure, symmetry, states and energy

#### III.3.1 Vacancy

There is an inward relaxation of the atoms around the vacancy site. The symmetry of the vacancy determines the impurity states’ symmetry. The singlet s-state is located deep inside the valence band (not shown in Fig. 3) while partly occupied triplet p-states are located in the bandgap (Fig.3). There is a splitting of p- states into twofold degenerate lower band and a single upper band (almost unoccupied). There are a few major differences between the LDA results and the LDA+U results. First, there is the expected increase in the bandgap in the LDA+U compared to the LDA leading to the separation between the vacancy states and valence and conduction bands. Second, the dispersion of vacancy p-states is smaller in the LDA+U calculation. Thus, these states are more localized in the LDA+U calculation. The dispersion cannot be removed in the supercell geometry completely, but the LDA definitely seems to overestimate the itinerancy of electrons in the vacancy states. There is no noticeable Jahn-Teller distortion of the atoms surrounding the vacancy in our 64-atom supercell geometry. This is consistent with the observation by Zywietz et. al. Zywietz et al. (1998) who suggest that a 128-atom supercell is required to observe the distortion. The formation and the migration energy for the vacancy defect from LDA+U calculations are 4.59 eV and 0.42 eV respectively. This is in good agreement with experimental results Ural et al. (1999b); Watkins (2000)

#### III.3.2 Hexagonal interstitial

The self interstitial sites provide four additional electrons to the system and create similar “dangling” bonds. Like the vacancy, they create an almost doubly degenerate band near the bottom of the band gap and a single band near the top of the band gap (Fig. 4). The occupied defect states are very close to the bottom of the band gap, in agreement with the observation made by Needs Needs (1999) about the defect states being quite shallow. We find that the atoms surrounding the interstitial move outwards in contrast to the inward movement reported by Needs Needs (1999). This is probably because, while we have optimized the lattice constant of the system with the defect, Needs has maintained the lattice constant at the experimental value. We suspect that Needs’ approach might cause an unphysical restriction on the relaxation. In addition to opening the band gap, the introduction of the on-site repulsion quite dramatically reduces the mixing of the defect states with the valence states. The formation energy for the hexagonal interstitial using the LDA+U method is 4.61 eV which is in good agreement with experimental Ural et al. (1999a) and quantum Monte Carlo Leung et al. (1999) results.

#### III.3.3 Split $`110`$ interstitial

The split $`110`$ interstitial creates defect states in the band gap similar to the hexagonal interstitial excpet that there is a greater splitting of the two lower bands (Fig. 5). We attribute this to the lower symmetry of the split $`110`$ interstitial compared to the hexagonal interstitial. Another difference is the location of the higher defect state. It is not as close to the top as in the case of hexagonal interstitial. We do not have a simple explanation for this. The introduction of the on-site repulsion has a similar effect with respect to reducing the mixing of the defect states with the bulk states. The defect formation energy using the LDA+U method is 4.70 eV which is in good agreement with the experimental Ural et al. (1999a) and quantum Monte Carlo Leung et al. (1999) results. The migration energy from the hexagonal to the split $`110`$ configurations of the interstitial was computed using the nudged elastic band menthod Jónsson et al. (1998). We get a value of 0.44 eV using the LDA+U method. We are not aware of any experimental result for this specific value.

### III.4 Self diffusion in Si

The contribution of the mechanism to the diffusion process is determined by the activation energy and the prefactor. Because all the three mechanisms considered in this article (viz., vacancy, interstitial, and concerted exchange) involve similar number of atoms, one can, to a first approximation, assume that all of them have similar entropic effects and hence similar prefactors. Thus, based on our calculations of the activation energies using the LDA+U method, the vacancy mechanism \[which has an activation energy of 5.01 eV (sum of formation (4.59 eV) and migration (0.42 eV) energies)\] and the hexagonal to split $`110`$ interstitialcy mechanism \[which has an activation energy of 5.14 eV (sum of formation (4.7 eV) and migration (0.44 eV) energies)\] are quite close in activation energy (within the error of the method), and should contribute equally to self-diffusion. The concerted exchange mechanism (which has an activation energy of 5.82 eV) is less significant in this respect. Experimental results give a similar indication Ural et al. (1999a); Fahey et al. (1989).

## IV Summary

The present work identified the cause for the poor description of point defect energetics by the LDA. It corrected for the deficiencies of the LDA by using the LDA+U method. This gave much better agreement of the calculated activation energies with experimental observations. This method can therefore be used for better description of diffusion in similar semiconductor materials.

###### Acknowledgements.

This work was funded by the United States Department of Energy, Basic Energy Sciences Grant No. DE-FG03-99ER45788. The computations were performed on the NPACI IBM p690, NERSC IBM SP RS/6000 (as part CMSN), and on the PC clusters of the Multiscale Simulation Laboratory, Stanford: multipod and xepod. Grants from NPACI and NERSC are sincerely appreciated.


format: text

text:
# Do probabilistic medium-range temperature forecasts need to allow for non-normality?

## 1 Introduction

We consider the question of how to make probabilistic forecasts of temperature at individual spatial locations. Our method is to take the output from a numerical model ensemble forecast as a predictor, and to use a statistical model to convert this predictor into a prediction of the likely temperature distribution. Our interest is in predictions of the whole distribution of temperatures, rather than more restricted questions such as whether the temperature will cross certain thresholds or whether the temperature will lie in certain categories.

In two previous articles (Jewson et al. (2003a) and Jewson (2003)) we have addressed this question using a statistical model that we will call the *gaussian spread regression* model. This model considers the conditional temperature distribution to be normal, and models the observed temperature in terms of the mean and the standard deviation of the ensemble forecast. The mean of the observed temperature distribution is taken to be a linear transformation of the ensemble mean, and the standard deviation of the observed temperature distribution is taken to be a linear transformation of the ensemble standard deviation. This model can be considered as a generalisation of standard linear regression: the extension is that a linear model, rather than just a constant, is used for the standard deviation.

There are a number of possible criticisms of the two papers cited above. In particular:

* the autocorrelations of the forecasts errors were ignored when fitting the model
* the model ignores any predictive information that might be present in the previous day’s temperatures
* the model assumes that temperature uncertainty (the conditional distribution of temperature) is normally distributed

We are continuing our research in an attempt to address each of these issues. This article tackles the third: the assumption of normality of the conditional temperature distribution. We present a new model which no longer makes this assumption. As with the gaussian spread regression model, this model calibrates the ensemble using linear transformations of the ensemble mean and the ensemble standard deviation. The difference is that the observed temperatures are then modelled using a *kernel density* rather than a normal distribution. Kernel densities allow much more flexible modelling of the probability density of the temperature, and can incorporate skewness and bimodality, both of which may exist in the distribution of forecast temperature uncertainty. We will call the new model *kernel spread regression*.

In Jewson (2003) we showed that there is a detectable relationship between spread and skill in one year of past forecasts from the ECMWF model for London Heathrow. We also showed that the gaussian spread regression model applied to these forecasts performs better than linear regression for in-sample tests. However, in out of sample tests we could not show that the model performed better. We argue that this is because one year of data is rather little for fitting the model parameters accurately, combined with the fact that the predictable variations in uncertainty are rather weak.

For these reasons, we see little point in performing out of sample tests on the models we present below. The currently available data is simply too short to be able to distinguish between these models in such tests, since the results of the models are only subtly different. We will therefore only perform in-sample tests and the goodness of the model will be assessed on the basis of in-sample results. Testing using out of sample tests will have to wait until longer time series of stationary past forecasts are available.

In section 2 we will describe the forecast data to be used in this study. In section 3 we introduce the kernel density, and present some properties of kernel densities that will be useful in interpreting the results of our analysis. In section 4 we present the hierarchy of statistical calibration models that we will consider. The models are: gaussian regression, gaussian spread regression, kernel regression and kernel spread regression. In section 5 we present results from the use of these four models to calibrate some forecast data. Finally in section 6 we summarise our results, discuss areas for future work, and draw some conclusions.

## 2 Data

We will base our analyses on one year of ensemble forecast data for the weather station at London’s Heathrow airport, WMO number 03772. The forecasts are predictions of the daily average temperature, and the target days of the forecasts run from 1st January 2002 to 31st December 2002. The forecast was produced from the ECMWF model (Molteni et al., 1996) and downscaled to the airport location using a simple interpolation routine prior to our analysis. There are 51 members in the ensemble. We will compare these forecasts to the quality controlled climate values of daily average temperature for the same location as reported by the UKMO.

There is no guarantee that the forecast system was held constant throughout this period, and as a result there is no guarantee that the forecasts are in any sense stationary, quite apart from issues of seasonality. This is clearly far from ideal with respect to our attempts to build statistical interpretation models on past forecast data but is, however, unavoidable: this is the data we have to work with.

Throughout this paper all equations and all values have had both the seasonal mean and the seasonal standard deviation removed. Removing the seasonal standard deviation removes most of the seasonality in the forecast error statistics, and partly justifies the use of non-seasonal parameters in the statistical models for temperature that we propose.

## 3 The kernel density

Two of the four models that we consider in this study make use of kernel densities. The kernel density is a very flexible method for deriving a density function from observed data. It works by assigning an identical kernel to each individual observed data point. The fitted density is then the sum of these individual kernels. There is a single free parameter, which is the width of the kernel, known as the *bandwidth*. The kernels are normalised so that the total density is one. A number of different kernel shapes are possible, but for simplicity, and to ensure that our densities never give zero probabilities, we will use gaussian kernels. This gives a probability density of the form:

$`p(x)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{i=M}{}}}p_i(x)`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle \frac{1}{\sqrt{2\pi }\lambda M}}\text{exp}\left({\displaystyle \frac{(xx_i)^2}{2\lambda ^2}}\right)`$
where $`\lambda `$ is the bandwidth, $`M`$ is the number of data points, and $`x_i`$ are the data points themselves. In our case the data points $`x_i`$ will be the individual ensemble members, after certain transformations that are described in section 4.

We now derive an important relationship which links the sample variance of the data points $`x_i`$ to the variance of the fitted kernel distribution $`p(x)`$.

If $`x`$ is the temperature, then the expectation of $`x`$ calculated using the probabilities from the calibrated forecast $`p(x)`$ is:

$$\mu =xp(x)𝑑x$$
(2)
But $`p(x)`$ is given by a kernel density, and so:

$$p(x)=\underset{i=1}{\overset{i=M}{}}p_i(x)$$
(3)
where

$$p_i(x)𝑑x=\frac{1}{M}$$
(4)
and so:

$`\mu `$ $`=`$ $`{\displaystyle xp(x)𝑑x}`$
$`=`$ $`{\displaystyle x\left(\underset{i=1}{\overset{i=M}{}}p_i(x)\right)𝑑x}`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle xp_i(x)𝑑x}`$
$`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle xMp_i(x)𝑑x}`$
The normalisation of $`p_i(x)`$ (equation 4) means that we can consider $`Mp_i(x)`$ to be a density with mean of $`x_i`$, where $`x_i`$ is the $`i`$’th ensemble member, and so:

$$\mu =\frac{1}{M}\underset{i=1}{\overset{i=M}{}}x_i$$
(6)
In other words, the mean of $`x`$ under the density $`p(x)`$ is just the arithmetic mean of the ensemble members $`x_i`$.

Similarly we can consider the variance of $`x`$ under the density $`p(x)`$:

$`var(x)`$ $`=`$ $`{\displaystyle (x\mu )^2p(x)𝑑x}`$
$`=`$ $`{\displaystyle (x\mu )^2\left(\underset{i=1}{\overset{i=M}{}}p_i(x)\right)𝑑x}`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle (x\mu )^2p_i(x)𝑑x}`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle (xx_i+x_i\mu )^2p_i(x)𝑑x}`$
$`=`$ $`{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle [(xx_i)^2+(x_i\mu )^2+2(xx_i)(x_i\mu )]p_i(x)𝑑x}`$
$`=`$ $`{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle (xx_i)^2Mp_i𝑑x}+{\displaystyle \frac{1}{M}}{\displaystyle \underset{i=1}{\overset{i=M}{}}}{\displaystyle (x_i\mu )^2Mp_i𝑑x}+{\displaystyle \frac{2}{M}}{\displaystyle \underset{i=1}{\overset{i=M}{}}}(x_i\mu ){\displaystyle (xx_i)p_i𝑑x}`$
The summand in the first of these terms is just the variance of $`x`$ under the density $`Mp_i`$. This is equal to $`\lambda ^2`$, since $`Mp_i`$ is a normal distribution with variance $`\lambda ^2`$. In the second term, we note that $`Mp_i(x)𝑑x=1`$. In the third term, we note that $`(xx_i)Mp_i(x)𝑑x=0`$ by the definition of $`x_i`$.

Putting this all together:

$$var(x)=\lambda ^2+\frac{1}{M}\underset{i=1}{\overset{i=M}{}}(x_i\mu )^2$$

But the second of these terms is just the sample variance of the ensemble members, and so we have:

$$\text{variance of modelled temperatures}=\lambda ^2+\text{sample variance of ensemble members}$$
(8)
What we see is that the variance of the fitted density is equal to the sample variance plus the bandwidth squared. This relation will help us reconcile results from the gaussian and kernel based models given below.

## 4 The kernel spread regression calibration model

We now describe each of the models we will use in this study in turn. These models are: gaussian regression, gaussian spread regression, kernel regression and kernel spread regression. We will use these four models to compare the benefits of using the ensemble spread as a predictor with the benefit of allowing non-normal, rather than just normal, forecast distributions.

### 4.1 Gaussian regression

The simplest model we will use for observed temperatures is just linear regression: we will call this *gaussian regression* to distinguish it from the kernel regression model described below, which replaces the normal distribution with a kernel density. We will write the gaussian regression model as:

$$T_i=N(\alpha +\beta m_i,\gamma )$$
(9)
$`\alpha ,\beta ,\gamma `$ are the parameters of the model, and $`m_i`$ is the ensemble mean. This model, and its performance on this data set, has been discussed extensively in Jewson (2003).

### 4.2 Gaussian spread regression

The second model we will use to model the observed temperatures is the gaussian spread regression model of Jewson et al. (2003a). We will write this model as:

$$T_i=N(\alpha +\beta m_i,\gamma +\delta s_i)$$
(10)
Relative to the gaussian regression model there is an extra parameter $`\delta `$ which scales the extra predictor $`s_i`$ (the ensemble standard deviation). The performance of this model, and the interpretation of the parameters, has been discussed extensively in Jewson et al. (2003a) and Jewson (2003).

The gaussian regression and gaussian spread regression models are included here as a basis for comparison for the kernel based models.

### 4.3 Kernel regression

The first of the two new models that we will present here is an adaption of standard gaussian regression that uses a kernel density rather than a normal distribution. The motivation for this is that the uncertainty in temperature forecasts is not necessarily normally distributed (indeed, is undoubtedly not *exactly* normally distributed) and so it does not necessarily make sense to apply a gaussian model. In the kernel regression model we predict the mean of future temperatures in exactly the same way as we do with the gaussian model: our prediction is a linear transformation of the ensemble mean. The difference from gaussian regression is that the distribution around the mean is derived from the individual ensemble members using the kernel density. To be consistent with this the parameters of the model are fitted using the likelihood calculated from the kernel density rather than the normal distribution.

We will write this model as:

$$T_i=K(\alpha +\beta m_i,\gamma ,\lambda )$$
(11)
where the notation has been chosen to emphasize the similarity with gaussian regression, but with the extra bandwidth parameter $`\lambda `$. The way this model works in detail is:

* The ensemble members are transformed using the linear transformation $`\alpha +\beta x_i`$.
* The standard deviation of the ensemble members is set to $`\gamma `$
* A kernel of bandwidth $`\lambda `$ is placed around each of the transformed ensemble members
* The final density is given by the sum of these kernels

We note that this model, along with the gaussian regression model, ignores any information in the ensemble standard deviation because the standard deviation of the members is forced to a fixed value of gamma.

### 4.4 Kernel spread regression

The second of the two new models that we present generalises the kernel regression model in the same way that gaussian spread regression generalises gaussian regression. In other words, this model extends the kernel regression model to include the ensemble spread as a predictor for the uncertainty of the predicted temperatures. This model is intended to have all the features of the gaussian spread regression model, but also to have the advantage that it can be used for forecast uncertainty that is not close to normally distributed, perhaps because of skewness or bimodality. We write this model as:

$$T_i=K(\alpha +\beta m_i,\gamma +\delta s_i,\lambda )$$
(12)
### 4.5 Parameter fitting

We fit the parameters of all four of our models by finding those values that maximise the likelihood of the observations given the model for our one year of data. We perform this fitting using an iterative optimisation scheme known as Newton’s method. For the gaussian models, this is very fast. For the kernel models it is significantly slower because evaluation of the probability density of the models involves the summing over kernels for each of the ensemble members. However, fitting all 10 leads still only takes a matter of seconds on a rather slow personal computer (100Mhz CPU).

## 5 Results

We now present the results of fitting our four models to the forecasts and the observed temperature data. We start with the results for the gaussian and kernel regression models and then show results for the spread regression models.

### 5.1 Regression model parameters

Parameters from the fitting of the gaussian and kernel regression models are shown in figure 1. These models share the parameters $`\alpha ,\beta ,\gamma `$, while the kernel regression model has the extra parameter $`\lambda `$. We see that the values of $`\alpha ,\beta `$ are effectively the same for the two models. The $`\gamma `$ parameter, however, is significantly different: it is much lower for the kernel regression than it is for the gaussian regression. The reason for this seems to lie in equation 8. In the gaussian model, the variance of temperature is captured by $`\gamma `$, while in the kernel model the variance is captured through a combination of $`\gamma `$ and the bandwidth. Thus the $`\gamma `$ can be partially offset by the bandwidth in this model. We will explore later, for the spread regression model, whether the resulting calibrated variances are the same or not.

That the bandwidth is offsetting the $`\gamma `$ raises the possibility that the model is overspecified because these two parameters are actually trying to capture the same thing, and that only one would do. It would certainly be possible to fit a kernel based model that had a *fixed* bandwidth…presumably the $`\gamma `$ parameter would then adjust to fit the observed variance. Similarly it would be possible to fit a kernel model that had a fixed value for $`\gamma `$: presumably the bandwidth would then adjust to fit the observed variance. However, it does seem that the presence of two independent parameters is beneficial and that the two parameters can potentially play different roles in the calibration. This is because the bandwidth modulates the amount of non-normality. For instance, the bandwidth can adjust to make the total density more or less unimodal (large values) or multimodel (small values).

### 5.2 Spread regression model parameters

We now consider the parameters from the gaussian and kernel spread regression models. Figure 2 shows that the values of $`\alpha `$ and $`\beta `$ are effectively the same as for the regression models. As with the regression models, the value of $`\gamma `$ for the kernel model is much lower than that for the gaussian model. The value of $`\delta `$ is also different between the two models. The extent to which the different values of $`\gamma `$ and $`\delta `$ lead to the same or different final levels for the mean and variability of the calibrated spread is discussed below.

Figure 3 shows the parameters for the kernel spread regression model alone, but with estimated confidence intervals (estimated in the standard way using the Fisher information). We see that the parameters for the ensemble mean are very well estimated (so much so that the confidence intervals are almost hidden by the data itself) while the parameters for the standard deviation are much more uncertain, particularly the delta parameter.

The value of $`\lambda `$ for the spread regression model is shown in the upper panel of figure 4, along with confidence intervals. This parameter is apparently very well estimated.

### 5.3 Density examples

For illustration we now show, in figure 5, two examples of the densities that are produced by the kernel spread regression. These examples are both from lead 10, and show the densities predicted by the model on the days with the highest and lowest ensemble standard deviations. The same bandwidth $`\lambda `$ is applied to these two cases by the model, and so we would expect that the high standard deviation case might show some multimodality, while the low standard deviation case is unlikely to. This is exactly what we see. We note, however, that the presence of multimodality in the wider of the two densities is not at all an indication that we have detected multimodality in the forecast: it is simply the output of the model on this particular day. The optimum value of the bandwidth is based on all days, and the optimisation process used for fitting the parameters may have ’sacrificed’ the goodness of fit on this particular day in order to benefit the likelihood score based on the entire data set.

### 5.4 Non-normality indicator

If the distribution of forecast uncertainty is frequently significantly non-normal one would expect the bandwidth in the kernel spread regression model to be small relative to the total standard deviation of the calibrated forecast. The lower panel of figure 4 shows the ratio of these two terms. In fact, this ratio is never small, indicating only a low level of non-normality is present in our data-set. The decrease in this ratio with lead time suggests that the non-normality increases slightly at longer leads, as one would expect.

### 5.5 Statistics of the calibrated spread

We now look at the statistics of the uncertainty of the calibrated forecast predicted by the two spread regression models. In particular, we consider the mean level of the uncertainty, and the standard deviation of the uncertainty. These are both shown in figure 6. In the upper panel we see that the mean level of the uncertainty in the two models is effectively the same (the lines are virtually indistinguishable). Even though the values of $`\gamma `$ are different, the value of $`\lambda `$ in the kernel spread regression model is exactly compensating so that the two models produce the same mean level of uncertainty. We also see that the same is broadly true for the variability of the uncertainty (lower panel) although with some small differences between the models. We suspect that these differences are simply due to sampling errors: fitting the $`\delta `$ parameter is much more difficult than fitting the $`\gamma `$ parameter because it is associated with the second moment of the second moment of the temperatures. We conclude that even though the values of $`\delta `$ are different in the two models, the combined effect on the variability of the calibrated spread is the same.

### 5.6 Likelihood scores

Our final set of results, in figure 7, compare the performance of the various models using the log-likelihood. This tells us which of the models does best in representing the observed temperatures, and is the acid-test of whether modelling non-normality confers any benefit. Because the various models do not have the same number of parameters, we have corrected the log-likelihoods using the BIC criterion (although the differences in the numbers of parameters is so small relative to the amount of data this makes no qualitative difference).

The likelihoods from the spread regression models are slightly better than those from the regression models. This is presumably because the uncertainty in the temperature forecast shows a small amount of flow dependence, and this flow dependence has been partly predicted by the forecast model. This issue is discussed in more detail in Jewson (2003).

If the distribution of prediction uncertainty were significantly non-normal (e.g. skewed or multimodal) and the forecasts were predicting this non-normality, we might expect that the likelihood values for the kernel models would be significantly better than those for the gaussian models. However, we do not see this in the results. In fact, the differences between the gaussian and kernel models are tiny. This suggests that non-normality in the distribution of forecast uncertainty is either not particularly strong, or that the forecasts are not capturing it (although we cannot distinguish between these two possibilities).

On very close inspection, we see that the log-likelihoods for the kernel models are *slightly* better than those for the gaussian models at long leads. This corresponds to the idea that non-normality in the distribution of forecast uncertainty is more likely at long leads, as non-linear processes will have had more time to take effect. However, a *much* longer data set of past forecasts would be needed to show that this effect is not just an artefact of sampling error. In any case, the difference would seem to be too small to be important.

## 6 Conclusions

We have extended both linear regression and the gaussian spread regression model of Jewson et al. (2003a) to cope with more general distribution shapes by using kernel densities instead of the normal distribution. This gives us a hierarchy of four models: gaussian regression, gaussian spread regression, kernel regression and kernel spread regression. We apply all four models to the calibration of one year of ensemble forecasts. The parameters of the linear transformation of the ensemble mean are the same in all the models. However, the parameters of the transformation of the standard deviation of the ensemble are different. This appears to be because part of the variance of the observed temperatures is taken up by the bandwidth parameter of the kernel models. Nevertheless, although the parameter values change, the mean and standard deviation of the final calibrated uncertainty is more or less the same in the gaussian and kernel models.

We compare the likelihoods attained by the four models. The spread regression models perform marginally better than the regression models. But the kernel models do not perform better than the gaussian models, except for *very slightly* higher likelihoods at long leads. We conclude that either the forecast uncertainty is not non-gaussian, or if it is, then the forecast does not capture that. Any non-normality in the forecast ensemble does not improve the forecasts, and can be ignored.

Which of the gaussian and kernel models should be used in practice? The gaussian models are certainly the simpler of the two and only involve manipulating the ensemble mean and standard deviation rather than ensemble members. They are easier and faster to fit, and the final distribution is easier to communicate as it can be summarised by the mean and the standard deviation. However, the kernel models are more general. We have only tested one station: there may be other locations for which the kernel model shows a greater benefit, and they are unlikely to do worse. These arguments are not conclusive in either direction, and which model to use seems to be a matter of personal preference at this point <sup>2</sup><sup>2</sup>2We prefer the gaussian models, because of their simplicity.. Only the availability of longer stationary time series of past forecasts could allow us to come to a firm conclusion as to whether one model is really better than the other.

There are still some aspects of the calibration of temperature forecasts that need investigation. One is the need for better modelling of the covariance matrix of forecast errors. Jewson (2003) showed that the residuals from the gaussian spread regression model are certainly not white in time, and this should really be taken into account in the model: this could possibly affect the results significantly. We note that the modelling of such autocorrelated residuals would be much easier in the context of the gaussian models than the kernel models.

The models we have considered are parametric. It is possible that non-parametric calibration models may also be possible, and may have some advantages over their parametric cousins. This would be worth investigating. However, the main benefit of non-parametric models is in the flexible range of distributions that they can represent. Given that our results show no material benefit from reducing the dependence on the normal distribution and using a much more flexible distributional shape, we doubt that non-parametric models would help either.

Finally, we mention the question of extending the kinds of calibration methods and models used in this study to other variables such as wind and precipitation. Neither wind nor precipitation are close to normally distributed, and it would not be appropriate to try and calibrate them using the gaussian models. However, for wind at least, it may be possible that the kernel models would work quite well. We plan to test this just as soon as we can get our hands on some ensemble wind forecasts and wind observations.


text:
# Two formation paths for cluster dwarf galaxies?

## 1. The origin of dwarf galaxies in clusters

A correlation between metallicity and mass (or luminosity) of galaxies has been established in numerous galaxy samples. including luminous early-type galaxies (e.g. Kobayashi & Arimoto 1999, Kuntschner 2000), spirals and irregular galaxies (Zaritsky et al. 1994, Lee et al. 2003), cluster dwarfs (Brodie & Huchra 1991) and the Local group (Mateo 1998). The fact that more luminous galaxies are more metal-rich is believed to be the main driver of a number of observed correlations, including the color-magnitude relation followed by early-type galaxies in clusters (Kodama et al. 1998).

In the last years we have carried out a photometric and spectroscopic study of galaxies in the Coma cluster, obtaining William Herschel Telescope spectra for about 300 cluster members including both giant and dwarf galaxies down to $`M_B14`$. One of the results of the spectral analysis was the discovery of a bimodal distribution in the metallicities of the faint cluster galaxies with signs of recent star formation (Poggianti et al. 2001). We identified a group of dwarfs with metallicity around solar or even supersolar (hereafter, metal-rich) and a group with \[M/H\] around -1.5 (metal-poor). All of these dwarfs have luminosity weighted ages of 3 Gyr or less. The metallicity and age estimates were based on the analysis of metallicity-sensitive and age-sensitive spectral indices in the Lick system.

Given the faintness of these galaxies and the relatively low signal-to-noise of our William Herschel spectra, errorbars on spectral indices and on the age and metallicity estimates were too large to draw any definite conclusion. Two additional results, however, lead us to believe that the metallicity bimodality could be real. First of all, in Poggianti et al. (2001) we obtained a coadded spectrum of the metal-rich and the metal-poor dwarfs separately and we showed that, in addition to the magnesium indices on which the metal-rich versus metal-poor division was based, most of the other metallicity-sensitive indices were indeed stronger in the average spectrum of the metal-poor galaxies. Secondly, the two groups of galaxies also clearly separate in the color-magnitude diagram, with the metal-poor dwarfs being significantly bluer than the metal-rich ones at similar luminosities.

In order to either confirm or disprove the existence of a metallicity spread, we have started a multislit programme with FOCAS at the Subaru telescope to obtain high signal-to-noise spectra. Only a few spectra of sufficient quality could be obtained so far due to weather conditions during the first Subaru run.

Line indices have been measured and are currently being compared with spectrophotometric models. These first measurements are encouraging, confirming that metallicity-sensitive indices are indeed very different, despite the similar magnitudes of these galaxies. However, only the simultaneous analysis of several line indices will establish whether this corresponds to large variations in metal content or/and ages or/and chemical abundance ratios.

At this point the existence of a bimodal distribution of metallicities among faint galaxies in Coma is an intriguing possibility suggested by spectral indices and supported by colors. Should the metallicity spread at faint magnitudes be confirmed by the high signal-to-noise spectra, it would have far reaching consequences for our understanding of early-type galaxy formation and evolution, and the hypothesis of different formation processes for faint galaxies in clusters should be seriously contemplated.

The metal-poor dwarfs could represent a “truly primordial” population, while the metal-rich ones could be the evolved remnants of more massive galaxies that lost some of their mass (and luminosity) due to environmental effects. Alternatively, the metal-rich dwarfs could be in some way related with another class of dwarf galaxies that are known to deviate from the metallicity/luminosity relation, created from the tidal debris of mergers between metal-rich massive galaxies (tidal dwarfs). These and other speculative options are currently open.

## References

Brodie, J.P. & Huchra, J.P., 1991, ApJ, 379, 157

Kobayashi, C., Arimoto, N., 1999, ApJ, 527, 573

Kodama, T., Arimoto, N., Barger, A.J., Aragon-Salamanca, A., 1998, A&A, 334, 99

Kuntschner, H., 2000, MNRAS, 315, 184

Lee, H., McCall, M.L., Kingsburgh et al., 2003, AJ, 125, 146

Mateo, M.L., 1998, ARAA, 36, 435

Poggianti, B. M., Bridges, T. J., Mobasher, T. J. et al., 2001, ApJ, 562, 689

Zaritsky, D. Kennicutt, R. C., Huchra, J.P., 1994, ApJ, 420, 87


text:
# Non trivial behavior of the linear response function in phase ordering kinetics1footnote 11footnote 1Talk delivered by Marco Zannetti

## I Introduction

It is a great pleasure to present this talk on the occasion of the 60th birthday of Francesco Guerra. My interest in relaxation phenomena, in fact, initiated quite a few years ago in Salerno where Nelson stochastic mechanics was a subject very much cultivated under the guide of Francesco. I have benefited a lot from his teaching and from the very stimulating atmosphere he was able to create about the applications of stochastic methods to physics in general.

Coming to the subject of this talk, phase ordering Bray94  is usually regarded as the simplest instance of slow relaxation, quite useful for a first and easy to understand approach to concepts like scaling and aging, which are the hallmarks of glassy behavior Cugliandolo2002 . However, next to the similarities there are also fundamental differences which require to keep phase ordering well distinct from the out of equilibrium behavior in glassy systems, both disordered and non disordered. The main source of the differences is the simplicity of the free energy landscape in the case of phase ordering compared to the complexity underlying glassy behavior. This leads to a simple low temperature state for phase ordering systems as opposed to the complexity of replica symmetry breaking for glasses (at least in the mean field picture of spin glasses). An analogous sharp distinction is believed to exist also in the behavior of the out of equilibrium linear response function. As a matter of fact, in the framework where static and dynamic properties are connected Franz98 , systems may be classified on the basis of the fluctuation dissipation relation Ricci99 .

In the jargon of slow relaxation, phase ordering is frequently referred to as coarsening. So, if coarsening is associated to simplicity and triviality of behavior what is there to investigate about it? Besides the obvious motivation that the basic, paradigmatic cases need to be thouroughly understood, an additional reason, among others, is that in some cases the existence of complex slow relaxation is identified through the exclusion of coarsening. This requires precise and reliable knowledge of all what goes on when relaxation proceeds by coarsening. An example comes from the long standing controversy about the nature of the low temperature phase of finite dimensional spin glasses. One recent argument in favour of replica symmetry breaking is that the observed behavior of the response function is incompatible with coarsening Franz98 ; Ricci99 . This might well be the case. However, for the argument to be sound, the understanding of the out of equilibrium behavior of the response function during phase ordering needs to be up to the level that such a delicate issue demands.

It is the purpose of this talk to present an overview of the accurate investigation of the response function in phase ordering that we have carried out in the last few years. Focusing on the zero field cooled magnetization (ZFC), which is one instance of integrated response function, it will be argued that the response function in phase ordering systems is not as trivial as it is believed to be and, after all, it is not the quantity best suited to highlight the differences between systems with and without replica symmetry breaking. In fact, as we shall see, phase ordering and therefore a replica symmetric low temperature state are compatible with a non trivial ZFC. When this happens there is no connection between static and dynamic properties. Phase ordering systems offer examples of two distinct mechanism for the lack of this important feature of slow relaxing systems, the vanishing of the scaling exponent of ZFC and stochastic instability.

## II Phase ordering

Let us first briefly recall the main features of a phase ordering process. Consider a system, like a ferromagnet, with order parameter (vector or scalar, continous or discrete) $`\varphi (\stackrel{}{x})`$ and Hamiltonian $`[\varphi (\stackrel{}{x})]`$ such that below the critical temperature $`T_C`$ the structure of the equilibrium state is simple. For example, in the scalar case, there are two pure ordered states connected by inversion symmetry. The form of the Hamiltonian can be taken the simplest compatible with such a structure, like Ginzburg–Landau–Wilson (GLW) for continous spins or the nearest neighbors Ising Hamiltonian for discrete spins.

In the following we will be interested in the space and time dependent correlation function

$$C(\stackrel{}{r},t,s)=\varphi (\stackrel{}{x},t)\varphi (\stackrel{}{x}^{},s)\varphi (\stackrel{}{x},t)\varphi (\stackrel{}{x}^{},s)$$
(1)
where the average is taken over initial condition and thermal noise, $`\stackrel{}{r}=\stackrel{}{x}\stackrel{}{x^{}}`$ and $`ts0`$ are two times after the quench. The conjugated linear response function is given by

$$R(\stackrel{}{r},t,s)=\frac{\delta \varphi (\stackrel{}{x},t)}{\delta h(\stackrel{}{x^{}},s)}|_{h=0}$$
(2)
and ZFC is defined by

$$\chi (\stackrel{}{r},t,t_w)=_{t_w}^t𝑑sR(\stackrel{}{r},t,s).$$
(3)
### II.1 Dynamics over phase space: equilibration versus falling out of equilibrium

For a temperature $`T`$ below $`T_C`$, in the thermodynamic limit, the phase space $`\mathrm{\Omega }=\{[\varphi (\stackrel{}{x})]\}`$ may be regarded as the union of three ergodic components Palmer  $`\mathrm{\Omega }=\mathrm{\Omega }_+\mathrm{\Omega }_{}\mathrm{\Omega }_0`$, where $`\mathrm{\Omega }_\pm `$ and $`\mathrm{\Omega }_0`$ are the subsets of configurations with magnetization $`lim_V\mathrm{}\frac{1}{V}_V𝑑\stackrel{}{x}\varphi (\stackrel{}{x})`$ positive, negative and vanishing, respectively. Denoting by $`\rho _\pm [\varphi (\stackrel{}{x})]`$ the two broken symmetry pure states, all equilibrium states are the convex linear combinations of $`\rho _\pm `$. In particular, the Gibbs state is the symmetric mixture $`\rho _G[\varphi (\stackrel{}{x})]=\frac{1}{Z}\mathrm{exp}([\varphi (\stackrel{}{x})]/T)=\frac{1}{2}\rho _+[\varphi (\stackrel{}{x})]+\frac{1}{2}\rho _{}[\varphi (\stackrel{}{x})]`$. The $`\mathrm{\Omega }_\pm `$ components are the domains of attraction of the pure states with $`\rho _+(\mathrm{\Omega }_+)=\rho _{}(\mathrm{\Omega }_{})=1`$ and $`\mathrm{\Omega }_0`$ is the border in between them, with zero measure in any of the equilibrium states.

When ergodicity is broken, quite different behaviors may arise Palmer  depending on the initial condition $`\rho _0[\varphi (\stackrel{}{x})]=\rho ([\varphi (\stackrel{}{x})],t=0)`$. Here, we consider the three cases relevant for what follows, assuming that there are not explicit symmetry breaking terms in the equation of motion:

1. equilibration to a pure state

if $`\rho _0(\mathrm{\Omega }_+)=1`$ or $`\rho _0(\mathrm{\Omega }_{})=1`$, in the time evolution configurations are sampled from either one of $`\mathrm{\Omega }_\pm `$ and $`\rho ([\varphi (\stackrel{}{x})],t)`$ equilibrates to the time independent pure state $`\rho _\pm [\varphi (\stackrel{}{x})]`$ within the finite relaxation time $`t_{eq}\xi ^z`$, where $`\xi `$ is the equilibrium correlation lenght and $`z`$ is the dynamic exponent entering the domain growth law, to be defined shortly. The correlation function is the same in the two ergodic components and, after equilibration, is time translation invariant

$$C_{st}(\stackrel{}{r},ts)=\varphi (\stackrel{}{x},t)\varphi (\stackrel{}{x}^{},s)_\pm M^2$$
(4)
where $`\varphi (\stackrel{}{x})_\pm =\pm M`$ is the spontaneous magnetization. For large distances $`r\xi `$ and time separations $`tst_{eq}`$, the clustering property $`\varphi (\stackrel{}{x},t)\varphi (\stackrel{}{x}^{},s)_\pm `$ $`=\varphi (\stackrel{}{x},t)_\pm \varphi (\stackrel{}{x}^{},s)_\pm `$ is obeyed and the correlations decay to zero, as required by ergodicity.
2. equilibration to the Gibbs state

if $`\rho _0(\mathrm{\Omega }_+)=\rho _0(\mathrm{\Omega }_{})=1/2`$, then configurations are sampled evenly from both disjoint components $`\mathrm{\Omega }_+`$ and $`\mathrm{\Omega }_{}`$. The probability density $`\rho ([\varphi (\stackrel{}{x})],t)`$ equilibrates now to the Gibbs state $`\rho _G[\varphi (\stackrel{}{x})]`$ with the same relaxation time $`t_{eq}`$ as in the relaxation to the pure states. Broken ergodicity shows up in the large distance and in the large time properties of the correlation function. After equilibration, one has

$$C_G(\stackrel{}{r},ts)=C_{st}(\stackrel{}{r},ts)+M^2$$
(5)
from which follows that correlations do not vanish asymptotically or that the clustering property is not obeyed

$$\underset{r\mathrm{}}{lim}C_G(\stackrel{}{r},ts)=\underset{(ts)\mathrm{}}{lim}C_G(\stackrel{}{r},ts)=M^2.$$
(6)
3. falling out of equilibrium over the border Laloux ; Newman

if $`\rho _0(\mathrm{\Omega }_0)=1`$, for the infinite system $`\rho (\mathrm{\Omega }_0,t)=1`$ also at any finite time after the quench. Namely, the system does not equilibrate since in any equilibrium state the measure of $`\mathrm{\Omega }_0`$ vanishes. Phase ordering corresponds to this case. In fact, the system is initially prepared in equilibrium at very high temperature (for simplicity $`T_I=\mathrm{}`$) and at the time $`t=0`$ is suddenly quenched to a final temperature $`T`$ below $`T_C`$. In the initial state the probability measure over phase space is uniform $`\rho _0[\varphi (\stackrel{}{x})]=1/|\mathrm{\Omega }|`$, implying that the initial configuration at $`t=0`$ belongs almost certainly to $`\mathrm{\Omega }_0`$, since with a flat measure $`|\mathrm{\Omega }_0|`$ is overwhelmingly larger than $`|\mathrm{\Omega }_\pm |`$.

The morphology of typical configurations visited as the system moves over $`\mathrm{\Omega }_0`$ is a patchwork of domains of the two competing equilibrium phases, which coarsen as the time goes on. The typical size of domains grows with the power law $`L(t)t^{1/z}`$, where $`z=2`$ (independent of dimensionality) for dynamics with non conserved order parameter Bray94 , as it will be consider here. The sampling of configurations of this type is responsible of the peculiar features of phase ordering. At a given time $`t_w`$ there remains defined a length $`L(t_w)`$ such that for space separations $`rL(t_w)`$ or for time separations $`tt_wt_w`$ intra-domains properties are probed. Then, everything goes as in the case 2 of the equilibration to the Gibbs state, ergodicity looks broken and the correlation function obeys Eq. (5). Conversely, for $`rL(t_w)`$ or $`t/t_w1`$, inter-domains properties are probed, ergodicity is restored (as it should be, since evolution takes place within the single ergodic component $`\mathrm{\Omega }_0`$) and eventually the correlation function decays to zero. However, the peculiarity is that if the limit $`t_w\mathrm{}`$ is taken before $`r\mathrm{}`$, in the space sector ergodicity remains broken giving rise, for instance, to the growth of the Bragg peak in the equal time structure factor.

According to this picture, the correlation function can be written as the sum of two contributions

$$C(\stackrel{}{r},t,s)=C_{st}(\stackrel{}{r},ts)+C_{ag}(\stackrel{}{r},t,s)$$
(7)
where the first one is the stationary contribution of Eq. (4) describing equilibrium fluctuations in the pure states and the second one contains all the out of equilibrium information. The latter one is the correlation function of interest in the theory of phase ordering where, in order to isolate it, zero temperature quenches are usually considered as a device to eliminate the stationary component. It is now well established that $`C_{ag}(\stackrel{}{r},t,s)`$ obeys scaling in the form Furukawa

$$C_{ag}(\stackrel{}{r},t,s)=\widehat{C}(r/L(s),t/s)$$
(8)
with $`\widehat{C}(x,y)=M^2`$ for $`x<1`$ and $`y1`$, while

$$\widehat{C}(r/L(s),t/s)(t/s)^{\lambda /z}h(r/L(s))$$
(9)
for large time separation Bray94 , where $`\lambda `$ is the Fisher–Huse exponent.

## III Zero field cooled magnetization

Let us next consider what happens when a time independent external field $`h(\stackrel{}{x})`$ is switched on at the time $`t_w`$. To linear order the expectation value of the order parameter at the time $`t`$ is given by

$$\varphi (\stackrel{}{x},t)_h=\varphi (\stackrel{}{x},t)_0+𝑑\stackrel{}{x^{}}\chi (\stackrel{}{x}\stackrel{}{x^{}},t,t_w)h(\stackrel{}{x^{}})$$
(10)
and if $`h(\stackrel{}{x})`$ is random with expecations $`\overline{h(\stackrel{}{x})}=0`$, $`\overline{h(\stackrel{}{x})h(\stackrel{}{x^{}})}=h_0^2\delta (\stackrel{}{x}\stackrel{}{x^{}})`$ then one has

$$\chi (\stackrel{}{x}\stackrel{}{y},t,t_w)=h_0^2\overline{\varphi (\stackrel{}{x},t)_hh(\stackrel{}{y})}.$$
(11)
Namely, ZFC is the correlation at the time $`t`$ of the order parameter with the random external field.

Going to the three processes considered above

1. after equilibration in the pure state has occurred and the stationary regime has been entered, the order parameter correlates with the external field via the equilibrium thermal fluctuations. The fluctuation dissipation theorem (FDT) is obeyed

$$\chi _{st}(\stackrel{}{r},tt_w)=\frac{1}{T}\left[C_{st}(\stackrel{}{r},ts=0)C_{st}(\stackrel{}{r},tt_w)\right]$$
(12)
and since $`C_{st}(\stackrel{}{r},tt_w)`$ decays to zero for $`tt_w>t_{eq}`$, over the same time scale $`\chi _{st}(\stackrel{}{r},tt_w)`$ saturates to

$$\underset{t\mathrm{}}{lim}\chi _{st}(\stackrel{}{r},tt_w)=\chi _{eq}(\stackrel{}{r})=\frac{1}{T}C_{eq}(\stackrel{}{r})$$
(13)
which is the susceptibility computed in the final equilibrium state $`\rho _\pm [\varphi (\stackrel{}{x})]`$.
2. As far as ZFC is concerned, there is no difference between the relaxation to the mixed Gibbs state and the relaxation to a pure state. Hence, FDT is satisfied and can be written both in terms of $`C_{st}`$ or $`C_G`$ since, as Eq. (5) shows, they differ by a constant.
3. In the phase ordering process the system stays out of equilibrium, so it useful to write ZFC as the sum of two contributions Bouchaud97

$$\chi (\stackrel{}{r},t,t_w)=\chi _{st}(\stackrel{}{r},tt_w)+\chi _{ag}(\stackrel{}{r},t,t_w)$$
(14)
where $`\chi _{st}(\stackrel{}{r},tt_w)`$ satisfies Eq. (12) and $`\chi _{ag}(\stackrel{}{r},t,t_w)`$ represents the additional out of equilibrium response. In connection with this latter contribution there are two basic questions

i) how does it behave with time

ii) what is the relation between $`\chi _{ag}`$ and $`C_{ag}`$, if any.

### III.1 Scaling hypothesis

Since ZFC measures the growth of correlation between the order parameter and the external field, the first question raised above addresses the problem of an out of equilibrium mechanism for this correlation, in addition to the thermal fluctuations accounting for $`\chi _{st}`$. Restricting attention from now on, for simplicity, to the case of coincident points ($`\stackrel{}{r}=0`$) and dropping therefore the space dependence, the starting point for the answer is the assumption of a scaling form

$$\chi _{ag}(t,t_w)t_w^{a_\chi }\widehat{\chi }_{ag}(t/t_w)$$
(15)
which is the counterpart of Eq. (8) for the correlation function.

The next step is to make statements on the exponent $`a_\chi `$ and on the scaling function $`\widehat{\chi }_{ag}(x)`$ . There exists in the literature an estimate of $`a_\chi `$ based on simple reasoning. What makes phase ordering different from relaxation in the pure or in the Gibbs state is the existence of defects. The simplest assumption is that $`\chi _{ag}(t,t_w)`$ is proportional to the density of defects Barrat98 ; Franz98 ; Ricci99 . This implies

$$a_\chi =\delta $$
(16)
where the exponent $`\delta `$ regulates the time dependence of the density of defects $`\rho _{defect}(t)L(t)^nt^\delta `$, namely

$$\delta =n/z$$
(17)
with $`n=1`$ for scalar and $`n=2`$ for vector order parameter Bray94 .

According to this argument $`a_\chi `$ should be independent of dimensionality. This conclusion is not corroborated by the available exact, approximate and numerical results. On the basis of exact analytical solutions for the $`d=1`$ Ising model Lippiello2000 ; Godreche  and for the large $`N`$ model Corberi2002 , approximate analytical results based on the gaussian auxiliary field (GAF) approximation Berthier99 ; Corberi2001  and numerical results from simulations of the Ising model Corberi2001 ; prl ; preprint  with $`d=2,3,4`$, the following general formula for $`a_\chi `$ has been obtained

$$a_\chi =\{\begin{array}{cc}\delta \left(\frac{dd_L}{d_Ud_L}\right)\text{for}d<d_U\hfill & \\ \delta \text{with log corrections for}d=d_U\hfill & \\ \delta \text{for}d>d_U\hfill & \end{array}$$
(18)
where $`d_L`$ and $`d_U>d_L`$ do depend on the system in the following way

* $`d_L`$ is the dimensionality where $`a_\chi =0`$. In the Ising model $`d_L=1`$, while in the large $`N`$ model $`d_L=2`$. The speculation is that in general $`d_L=1`$ for systems with discrete symmetry and $`d_L=2`$ for systems with continous symmetry, therefore suggesting that $`d_L`$ coincides with the lower critical dimensionality of equilibrium critical phenomena, although the reasons for this identification are far from clear.
* $`d_U`$ is a value of the dimensionality specific of ZFC and separating $`d<d_U`$, where $`a_\chi `$ depends on $`d`$, from $`d>d_U`$ where $`a_\chi `$ is independent of dimensionality and Eq. (16) holds. The existence of $`d_U`$ is due preprint  to a mechanism, i.e. the existence of a dangerous irrelevant variable, quite similar (including logarithmic corrections) to the one leading to the breaking of hyperscaling above the upper critical dimensionality in static critical phenomena. However, $`d_U`$ cannot be identified with the upper critical dimensionality since we have found, so far, $`d_U=3`$ in the Ising model and $`d_U=4`$ in the large $`N`$ model. In the scalar case it may be argued generic  that $`d_U`$ coincides with the dimensionality $`d_R=3`$ such that interfaces do roughen for $`dd_R`$ and do not for $`d>d_R`$. A general criterion for establishing the value of $`d_U`$ is not yet known.

The validity of Eq. (15) with $`a_\chi `$ given by Eq. (18) has been checked, in addition to the cases where analytical results are available, with very good accuracy in the simulations of the Ising model preprint . The values of $`\delta `$, $`d_L`$ and $`d_U`$ obtained for the different systems are collected in Table 1 and the behavior of $`a_\chi `$ as dimensionality is varied is displayed in Fig.1.

At this point, a comment is in order. From Eqs. (14) and (15) follows that the existence of out of equilibrium degrees of freedom, or defects, generates the aging contribution to ZFC. However, when $`a_\chi >0`$ this is bound to disappear as time becomes large, either $`t_w\mathrm{}`$ or $`t/t_w\mathrm{}`$. The vanishing of this contribution is fastest for $`d>d_U`$, where $`a_\chi =\delta `$, and becomes slower and slower below $`d_U`$, as $`d_L`$ is approached. Nonetheless, as long as $`a_\chi >0`$, eventually $`\chi _{ag}`$ disappears and Eq. (13) holds also for phase ordering processes for all $`d>d_L`$. This is no longer true for $`dd_L`$ where $`a_\chi 0`$ and, consequently, $`\chi _{ag}`$ gives a contribution to the response which persists also in the asymptotic time region leading to

$$\underset{t\mathrm{}}{lim}\chi (t,t_w)>\chi _{eq}$$
(19)
where $`\chi _{eq}`$ stands for $`\chi _{eq}(\stackrel{}{r}=0)`$. In the following we will consider $`dd_L`$.

### III.2 Fluctuation dissipation relation

Let us now come to the second question. A very useful tool for the study of slow relaxation phenomena has been introduced by Cugliandolo and Kurchan Cugliandolo93  through the off equilibrium fluctuation dissipation relation. This can be introduced in the following way. Given that $`C(t,t_w)`$ is a monotonously decreasing function of $`t`$, for fixed $`t_w`$ it is possible to invert it and write

$$\chi (t,t_w)=\stackrel{~}{\chi }(C(t,t_w),t_w).$$
(20)
Then, if for a fixed value of $`C(t,t_w)`$ there exists the limit

$$\underset{t_w\mathrm{}}{lim}\stackrel{~}{\chi }(C,t_w)=S(C)$$
(21)
the function $`S(C)`$ gives the fluctuation dissipation relation. In the particular case of equilibrium dynamics, FDT is recovered and $`S(C)=[C(0)C]/T`$. Originally introduced in the study of the low temperature phase of spin glass mean-field models, the fluctuation dissipation relation has been found in many other instances of slow relaxation Crisanti2002 .

Now, in order to search for $`S(C)`$ in the case of phase ordering, let us set $`\stackrel{}{r}=0`$ in Eq. (8) and let us eliminate $`t/t_w`$ between $`\widehat{\chi }_{ag}`$ and $`C_{ag}`$ obtaining

$$\chi _{ag}(t,t_w)t_w^{a_\chi }\stackrel{~}{\chi }_{ag}(C_{ag}).$$
(22)
Then, from Eqs. (14,12,22) one can write the general relation

$$\chi (t,t_w)=\frac{1}{T}\left[C_{st}(0)C_{st}(tt_w)\right]+t_w^{a_\chi }\stackrel{~}{\chi }_{ag}(C_{ag}).$$
(23)
Using the identity $`\left[C_{st}(0)C_{st}(tt_w)\right]=\left[C_{st}(0)+M^2C_{st}(tt_w)M^2\right]`$ and considering that in the time interval where $`C_{st}(tt_w)0`$, i.e. for short times, one can replace $`C_{ag}(t/t_w)`$ with $`M^2`$ or equivalently $`C_{st}(tt_w)+M^2=C(t,t_w)`$, the above equation can be rewritten as

$$\chi (t,t_w)=\stackrel{~}{\chi }_{st}(C)+t_w^{a_\chi }\stackrel{~}{\chi }_{ag}(C_{ag})$$
(24)
where the function $`\stackrel{~}{\chi }_{st}(C)`$ is defined by

$$T\stackrel{~}{\chi }_{st}(C)=\{\begin{array}{cc}\left[C(0)C(t,t_w)\right]\text{for}M^2CC(0)\hfill & \\ \left[C(0)M^2\right]\text{for}C<M^2.\hfill & \end{array}$$
(25)
Therefore, from Eq. (24) we have, first of all, that for phase ordering systems the fluctuation dissipation relation exists only if $`a_\chi 0`$ (i.e. for $`dd_L`$) and that for $`a_\chi >0`$

$$S(C)=\stackrel{~}{\chi }_{st}(C).$$
(26)
For $`a_\chi =0`$ a little more care is needed. Equation (24) yields $`\chi (t,t_w)=\stackrel{~}{\chi }_{st}(C)+\stackrel{~}{\chi }_{ag}(C_{ag})`$. Recalling that $`a_\chi =0`$ occurs at $`d=d_L`$, which coincides with the lower critical dimensionality, in order to have a phase ordering process a quench to $`T=0`$ is required. This, in turn, implies $`C_{st}(t,s)=0`$ and $`C_{ag}(t,s)=C(t,s)`$. Therefore, using Eq. (25) we have

$$S(C)=\chi _{eq}^{}+\stackrel{~}{\chi }_{ag}(C)$$
(27)
where $`\chi _{eq}^{}=lim_{T0}[C(0)M^2]/T`$ is the $`T=0`$ equilibrium susceptibility, which vanishes for hard spins while is different from zero for soft spins. Physical implications of these results are discussed in the next section.

## IV Statics from dynamics

One of the main reasons of interest in the fluctuation dissipation relation is that it may provide a link between static and dynamic properties. This was first found by Cugliandolo and Kurchan Cugliandolo93  for mean-field spin glasses and then established in general by Franz et al. Franz98  for slowly relaxing systems.

Let us first introduce the overlap probability function $`\stackrel{~}{P}(q)`$ in the equilibrium state obtained when the perturbation responsible of $`\chi (t,t_w)`$ is switched off. The question is how is $`\stackrel{~}{P}(q)`$ related to the unperturbed overlap function $`P(q)`$. If $`P(q)=\stackrel{~}{P}(q)`$ the system is stochastically stable Guerra . A milder statement of stochastic stability is that $`\stackrel{~}{P}(q)`$ coincides with $`P(q)`$ up to the effects of a global symmetry which might be removed by the perturbation. In the Ising case, where the perturbation breaks the up-down symmetry, defining

$$\widehat{P}(q)=2\theta (q)P(q)$$
(28)
the system is stochastically stable in the sense that $`\stackrel{~}{P}(q)=\widehat{P}(q)`$. If the system is not stochastically stable, $`\stackrel{~}{P}(q)`$ is not related neither to $`P(q)`$ nor to $`\widehat{P}(q)`$. As we shall see, this is the case of the mean spherical model.

The statement is that

1. if $`S(C)`$ exists
2. if $`lim_t\mathrm{}\chi (t,t_w)=\chi _{eq}`$
3. if the system is stochastically stable

then one has

$$T\frac{d^2S(C)}{dC^2}|_{C=q}=\stackrel{~}{P}(q).$$
(29)
We may now check if this chain of connections applies to phase ordering. In replica symmetric low temperature states, as for instance in ferromagnetic systems, the overlap function is always trivial and we have

$$P(q)=\frac{1}{2}\left[\delta (qM^2)+\delta (q+M^2)\right]$$
(30)
with

$$\stackrel{~}{P}(q)=\widehat{P}(q)=\delta (qM^2).$$
(31)
Computing the derivative in the left hand side of Eq. (29) and using Eqs. (26) and (27), for $`d>d_L`$ we find

$$T\frac{d^2S(C)}{dC^2}|_{C=q}=\delta (qM^2).$$
(32)
Therefore, Eq. (29) is satisfied.

### IV.1 Failure at $`d=d_L`$

Let us, next, go to the case $`d=d_L`$. For the sake of definiteness, we consider the case of the Ising model with $`d=1`$. In order to make compatible the two requirements of having an ordered equilibrium state and a well defined linear response function, instead of taking the $`T0`$ limit it is necessary to take the limit of an infinite ferromagnetic coupling Lippiello2000 . Then, $`P(q)`$ and $`\stackrel{~}{P}(q)`$ are given by Eqs. (30) and (31) with $`M^2=1`$ at all temperatures. On the other hand, for any $`T`$ we also have Lippiello2000

$$T\stackrel{~}{\chi }_{ag}(C)=\frac{\sqrt{2}}{\pi }\mathrm{arctan}\left[\sqrt{2}\mathrm{cot}\left(\frac{\pi }{2}C\right)\right].$$
(33)
This gives

$$T\frac{d^2S(C)}{dC^2}|_{C=q}=\frac{\pi \mathrm{cos}(\pi q/2)\mathrm{sin}(\pi q/2)}{[2\mathrm{sin}(\pi q/2)]^2}.$$
(34)
Hence, it is clear that Eq. (29) is not verified. The reason is that the second of the above conditions required for establishing the connection is not satisfied. In fact, from Eqs. (24) and (33), keeping in mind that the limits $`t\mathrm{}`$ and $`C0`$ are equivalent, we have

$$\underset{t\mathrm{}}{lim}T\chi (t,t_w)=1/\sqrt{2}$$
(35)
which is responsible of Eq. (19) and, therefore, of the violation of condition (2) above Eq. (29), since in this case $`\chi _{eq}=0`$.

### IV.2 Failure by stochastic instability

An interesting example Fusco , where statics cannot be reconstructed from dynamics because the third requirement of stochastic stability is not satisfied, comes from the spherical model. More precisely, one must consider in parallel the original version of the spherical model (SM) of Berlin and Kac Berlin  and the mean spherical model (MSM) introduced by Lewis and Wannier Lewis , with the spherical constraint treated in the mean. These two models are equivalent above but not below $`T_C`$ Kac . The low temperature states are quite different, with a bimodal order parameter probability distribution in the SM case and a gaussian distribution centered in the origin in the MSM case. The corresponding overlap functions are also very different Fusco . Considering, for simplicity, $`T=0`$ one has

$$P(q)=\{\begin{array}{cc}\frac{1}{2}\left[\delta (qM^2)+\delta (q+M^2)\right]\text{for SM}\hfill & \\ \frac{1}{\pi M^2}K_0(|q|/M^2)\text{for MSM}\hfill & \end{array}$$
(36)
where $`K_0`$ is a Bessel function of imaginary argument (Fig.2).

However, after switching off an external field, one finds for both models $`\stackrel{~}{P}(q)=\delta (qM^2)`$. This means that stochastic stability holds for SM but not for MSM.

On the other hand, the relaxation properties are the same in the two models, both above and below $`T_C`$ if the thermodynamic limit is taken before the $`t\mathrm{}`$ limit Fusco . Then, the linear response function is the same for both models and obeys Eq. (24) with $`a_\chi `$ given by Eq. (18), where $`\delta `$, $`d_L`$ and $`d_U`$ are the same as for the large $`N`$ model (Table I). Hence, we have that although Eq. (29) is satisfied for both models, nonetheless statics and dynamics are connected only in the SM case, where $`\stackrel{~}{P}(q)=\widehat{P}(q)`$. Instead, this is not possible in the MSM case where $`\stackrel{~}{P}(q)\widehat{P}(q)`$.

## V Conclusions

We have shown that although the low temperature equilibrium state in phase ordering systems is simple, ZFC in quenches below $`T_C`$ displays non trivial features. Summarising, these are

* The aging component $`\chi _{ag}(t,t_w)`$ obeys the scaling form (15) with an exponent $`a_\chi `$ dependent on dimensionality according to Eq. (18).
* The fluctuation dissipation relation $`S(C)`$ is trivial, in the sense that it is consistent with a replica symmetric equilibrium state, for $`d>d_L`$. However, as $`dd_L`$ from above the non trivial contribution due to $`\chi _{ag}(t,t_w)`$ persists for longer and longer times as $`a_\chi 0`$.
* For $`d=d_L`$ the fluctuation dissipation relation is non trivial, since the aging contribution does not disappear asymptotically. From this follows i) that there is no connection between static and dynamic properties and ii) that on the basis of the behavior of ZFC is not possible to establish that the equilibrium state is replica symmetric.
* The mean spherical model provides an explicit example of a system whose low temperature equilibrium state is not stochastically stable, producing the failure of the link between static and dynamic properties for all dimensionalities.
* Extensive numerical simulations for ferromagnetic systems with scalar, vector, conserved and non conserved order parameter, at various space dimensionalities are under way. Preliminary results do show that the behavior of the response function above illustrated might well be generic for phase ordering systems generic .

## Acknowledgements

This work has been partially supported from MURST through PRIN-2002.


text:
# Diffusion on a solid surface: Anomalous is normal

## Abstract

We present a numerical study of classical particles diffusing on a solid surface. The particles’ motion is modeled by an underdamped Langevin equation with ordinary thermal noise. The particle-surface interaction is described by a periodic or a random two dimensional potential. The model leads to a rich variety of different transport regimes, some of which correspond to anomalous diffusion such as has recently been observed in experiments and Monte Carlo simulations. We show that this anomalous behavior is controlled by the friction coefficient, and stress that it emerges naturally in a system described by ordinary canonical Maxwell-Boltzmann statistics.

PACS: 05.40-a, 68.43.Jk, 68.35.Fx

Diffusion of atoms, molecules, and clusters on solid surfaces occurs in a number of modern technologies involving self-assembled molecular film growth, catalysis, and surface-bound nanostructures nanoala . The study of the motion of small and large organic molecules phd ; prl1 , and of adsorbed metal clusters composed of tens and even hundreds of atoms phystoday ; prl2 , has led to the unexpected observation that, as with single atoms oneatom , long jumps may play a dominant role in these motions.

Theoretical, numerical, and phenomenological discussions of surface diffusion have led to the clear understanding that jumps beyond nearest neighbors are ubiquitous in some parameter regimes theory . However, these studies focus on the fact that the motion is necessarily diffusive (which is the case at very long time scales). The possibility that jumps can be so long as to lead to superdiffusive motion over appreciable intermediate time scales is recognized as an interesting problem, but one in which Lévy walks or flights mike  are invoked as a model input. AlthoughLévy-walk-like behavior is clearly observed in Hamiltonian systems Klafterb ; vega  and in microcanonical simulations vega ; Guantes , it is generally believed that such a fine signature of chaos is fully smeared away by thermal fluctuations. In the present work we show that Lévy-like statistics appear quite naturally over long time scales within the usual Langevin framework for underdamped motion in a periodic or a random potential.

While a detailed analysis of surface diffusion requires extensive calculations (e.g., ab initio, or molecular dynamics), even the most powerful currently available computers can not carry such calculations to anywhere near experimentally relevant time scales gross . Moreover, current experimental probes of the topography of surfaces, scanning tunneling microscopy and atomic force microscopy, are usually carried out at relatively high temperatures, which leads to additional difficulties for first-principles calculations. Therefore, simpler approaches are essential and valuable theory .

We consider a generic model of classical particles moving in a two-dimensional potential, under the action of thermal fluctuations and dissipation, the important control parameter being the friction coefficient. In spite of the simplicity of the model, we find that it is able to reproduce the entire range of experimentally and computationally observed phenomena, from superdiffusion all the way to subdiffusion.

The equation of motion of a particle of mass $`m`$ on the surface is

$$m\ddot{𝐱}=V(𝐱/\lambda )\mu \dot{𝐱}+\xi (t)$$
(1)
where $`\lambda `$ is the characteristic length scale of the potential. The parameter $`\mu `$ is the coefficient of friction, and the $`\xi _i(t)`$ are mutually uncorrelated white noises that obey the fluctuation-dissipation relation $`\xi _i(t)\xi _j(t^{})=2\mu k_BT\delta _{ij}\delta (tt^{}).`$ We first consider the nonseparable periodic potential

$$V(x,y)=V_0\mathrm{cos}\left(\frac{\pi x}{\lambda }+\frac{\pi y}{\lambda }\right)\mathrm{cos}\left(\frac{\pi x}{\lambda }\frac{\pi y}{\lambda }\right),$$
(2)
where $`V_0`$ is the barrier height at the saddle points.

Equation (1) can be rewritten in scaled dimensionless variables, $`r_x=x/\lambda `$, $`r_y=y/\lambda `$, and $`s=\sqrt{V_0/m\lambda }t`$, leaving only two independent parameters, the scaled temperature $`𝒯`$ and the scaled dissipation $`\gamma `$,

$$𝒯=k_BT/V_0,\gamma =\mu \lambda /\sqrt{mV_0}.$$
(3)
We study four properties of the motion of the particle: the mean square displacement, the dependence of the diffusion coefficient on friction, the probability density function of displacements, and the velocity power spectrum.

Normal diffusive behavior is characterized by a linear time dependence of the mean square displacement, $`r^2(s)s`$. Non-diffusive behavior shows a different time dependence, $`r^2(s)s^\alpha `$, with $`\alpha >1`$ ($`<1`$) for superdiffusive (subdiffusive) motion. In Fig. 1 we show typical trajectories obtained for two friction coefficients upon numerical simulation of the equations of motion with $`𝒯=0.2`$ (we use this value throughout). One (left panel) is for a large friction coefficient, and the particle follows typical diffusive motion characterized by short steps of length $`\lambda `$ and frequent changes in direction. The other (right panel) corresponds to a small friction coefficient and clearly shows the preponderance of long ($`\lambda `$) tracks along cartesian coordinates.

The evolution of $`r^2`$ averaged over $`1000`$ particles is shown in Fig. 2 for several friction coefficients. For very long times the motion is diffusive, as expected, but for small $`\gamma `$ and at intermediate times there is clear *superdiffusive* ballistic ($`\alpha =2`$) behavior over several decades in time. This is reflective of the long straight stretches seen in the low-$`\gamma `$ trajectory in Fig. 1. We stress that this behavior has emerged naturally and has not required explicit insertion of any but ordinary thermal fluctuations in the model.

Eventhough the motion of the particle may include long superdiffusive stretches, at long times the motion is necessarily diffusive. The dependence of the diffusion coefficient on the friction for small and for large $`\gamma `$ can be obtained analytically using the approximate relation $`Dl^2/2\tau `$, where $`l^2`$ is the mean square size of a jump out of one well and into another, and $`\tau ^1`$ is the mean jump rate (related to the familiar “mean escape rate”). In the overdamped regime, jumps typically occur from one well to a neighboring well, so $`l^21`$. Familiar Kramers formulas can be used to obtain the mean escape rate Hanggi , with the result

$$D\frac{1}{2\pi }\left(\sqrt{\frac{\gamma ^2}{4}+2\pi ^2}\frac{\gamma }{2}\right)e^{\frac{1}{𝒯}}\frac{\pi }{\gamma }e^{\frac{1}{𝒯}}.$$
(4)
The $`\gamma ^1`$ dependence of $`D`$ arises because $`l^2`$ is independent of $`\gamma `$ while $`\tau \gamma `$. In the underdamped limit, standard results Hanggi  lead to $`\tau 1/\gamma `$ sancho , but now $`l^2\gamma ^2`$ we . One obtains

$$D\frac{\pi 𝒯}{4\gamma }e^{\frac{1}{𝒯}},$$
(5)
that is, again an inverse dependence on friction. The theoretical diffusion coefficient as a function of the friction parameter is shown as the solid curves in the inset of Fig. 2. The symbols are the simulation results. The $`\gamma ^1`$ dependences have been noted in the literature theory , but we have provided explicit forms with *no adjustable parameters*.

The diffusion coefficient characterizes only the very long time asymptotic dynamics. The interesting intermediate dynamics in the low friction regime that gives rise to long stretches of ballistic motion is reflected in the probability distribution function (pdf) of particle displacements $`r`$. This pdf is shown in Fig. 3 for $`\gamma =0.0004`$ and tree different time intervals $`\tau _s`$. For comparison, we also show a typical pdf for high damping ($`\gamma =1`$) at the intermediate time interval. In the high-$`\gamma `$ curve the highest maximum corresponds to no jumps (by far the most likely event at short times). The next is associated with jumps to a nearest neighbor well, and so on. In contrast, the low-$`\gamma `$ curves show a very different behavior, with features strongly resembling those of a Lévy-walk model Klafterb ; Klaftera : a peak at small displacements, a power-law intermediate regime, and a side hump at high displacements. Each of these is a distinct signature of Lévy-walk-like dynamics, but one must be cautious in the detailed interpretation of these components. The persistent small displacement peak is associated with long trapping periods during which a particle does not move at all because its energy is not sufficient to overcome the barrier. The high displacement peak, which moves outward with velocity of order unity, is associated with ballistic motion of those particles that acquire enough energy to move (and lose it very slowly). Genuine Lévy-walk dynamics also exhibit a low displacement peak and a superdiffusive peak separated by a power law behavior, but there are some important differences. First, our distribution reflects ballistic transport in the intermediate regime (in the language of Ref. Klaftera , ballistic transport occurs when $`0<\alpha <1`$ in the Lévy model), whereas the regime where the Lévy model shows the features we have described is associated with sub-ballistic (but still superdiffusive) behavior (again, in the language of Ref. Klaftera , the behavior when $`1<\alpha <2`$). Second, the slope in our power law regime (approximately $`0.7`$) is not related to the exponent $`\alpha `$ in the mean square displacement as it is for the Lévy walk (where the slope is $`4\alpha `$). Third, our side hump is strongly broadened whereas the side hump in the Lévy-walk model is associated with motion at a single constant velocity. In our case the velocity varies according to the equilibrium Maxwell-Boltzmann distribution. Nevertheless, the qualitative features of our distribution track those of the Lévy walk. Note that the existence of the pronounced side hump reflects the fact that the particles performing long steps (“flights”) are those with a velocity in the tail of the Maxwellian distribution.

Long ballistic excursions imply velocity correlations over considerable time intervals. The velocity power spectrum $`S(\omega )=𝐯(\omega )𝐯(\omega )`$ for different values of $`\gamma `$ is shown in Fig. 4. The pronounced peak at $`\omega _0=\pi \sqrt{2}`$ is associated with small oscillations in one well. At lower frequencies, $`\omega \omega _0`$, one observes a power-law growth of $`S(\omega )`$ that corresponds to the persistent time correlations associated with the ballistic excursions. At even smaller frequencies there is a crossover to $`S(\omega )=`$const$`=D`$, indicating full decorrelation and pure diffusion.

Disorder in surfaces occurs due to the presence of vacancies and other defects. We have generated a random potential surface described by a Gaussian distribution with a correlation function $`V(𝐱)V(𝐱^{})=(\epsilon /2\pi \lambda ^2)\mathrm{exp}(|𝐱𝐱^{}|^2/2\lambda ^2)`$ (see ojalvo ; theirvarious ; ourvarious ; romero2  for details). We set the intensity $`\epsilon =100`$ and the characteristic length $`\lambda ^{}=4`$ in our simulations. A typical surface generated with this algorithm whose average height $`V_0`$ equals that of the periodic potential is shown in Fig. 5.

The exponent $`\alpha `$ in $`r^2(s)s^\alpha `$ at intermediate times shows the entire range of behaviors *from subdiffusive to superdiffusive* with changing friction. In Fig. 6 we show $`r^2(s)/4s`$, averaged over $`5000`$ particles, as a function of time, for several values of $`\gamma `$. In the overdamped regime we clearly observe *subdiffusive* behavior ($`\alpha <1`$romero2 , while *superdiffusive* ($`\alpha >1`$) behavior is seen for very small $`\gamma `$. The exponents $`\alpha `$, calculated over the last decade of the time variation of the mean square displacement within our finite simulation times, are plotted in the inset of Fig. 6 as a function of $`\gamma `$. Although the subdiffusive behavior is probably the true asymptotic behavior in the overdamped case romero2 , further theoretical and numerical support are needed to assess whether superdiffusion is the asymptotic behavior in the underdamped case.

We summarize our findings. We have explored the behavior of a particle in a two-dimensional potential described by ordinary Langevin dynamics under conditions of thermal equilibrium. In a periodic potential, in the underdamped regime, the motion of the particle includes a ballistic range that can extend over many decades of time. The pdf of the particle’s displacements under these conditions shows a structure strongly resembling one for Lévy walks. This may explain a number of observations involving superdiffusive motion of organic molecules prl1  and atomic clusters prl2  on surfaces without the need to invoke extraordinary fluctuations beyond the usual thermal description. The long-time behavior is diffusive in all cases, and we have been able to predict theoretically the dependence of the diffusion coefficient on friction over essentially the entire range of values of the friction parameter *with no adjustable parameters*. The situation in a random potential is even more complex, and exhibits a wide range of subdiffusive to superdiffusive regimes. Further analysis of the random potential case, and a more extensive presentation of the periodic problem, will be detailed elsewhere we .

This work was supported by the MCyT (Spain) under project BFM2003-07850, by the Engineering Research Program of the Office of Basic Energy Sciences at the U. S. Department of Energy under Grant No. DE-FG03-86ER13606, and by a grant from the University of California Institute for México and the United States (UC MEXUS) and the Consejo Nacional de Ciencia y Tecnología de México (CoNaCyT). A.H.R. acknowledges support from Millennium Initiative, Conacyt-Mexico, under Grant W-8001. I.M.S. acknowledges the hospitality of the University of Barcelona under the CEPBA grant, as well as partial financial support by the Fonds der Chemischen Industrie.