text:
# Dissipative chaotic scattering

## Abstract

We show that weak dissipation, typical in realistic situations, can have a metamorphic consequence on nonhyperbolic chaotic scattering in the sense that the physically important particle-decay law is altered, no matter how small the amount of dissipation. As a result, the previous conclusion about the unity of the fractal dimension of the set of singularities in scattering functions, a major claim about nonhyperbolic chaotic scattering, may not be observable.

Chaotic scattering  is a physical manifestation of transient chaos , which is due to the existence of nonattracting chaotic invariant sets, i.e., chaotic saddles, in the phase space. As a result, a Cantor set of singularities arises in physically measurable scattering functions relating an output variable after the scattering to an input variable before the scattering . Generally, the dynamics of chaotic scattering may be characterized as either hyperbolic or nonhyperbolic. In hyperbolic chaotic scattering, all the periodic orbits are unstable and there are no Kol’mogorov-Arnol’d-Moser (KAM) tori in the phase space and, as such, the survival probability of a particle in the scattering region typically decays exponentially with time. In nonhyperbolic chaotic scattering , there are both KAM tori and chaotic sets in the phase space. Due to the stickiness effect of KAM tori, a particle initialized in the chaotic region can spend a long time in the vicinity of KAM tori, leading to an algebraic decay  of the survival probability of the particle in the scattering region. A surprising result in nonhyperbolic chaotic scattering is that, because of the algebraic decay, the fractal dimension of the set of singularities in a scattering function is unity .

A physically important issue in the study of nonlinear dynamics is to understand how robust a phenomenon is against perturbations or deviations between the underlying mathematical model and physical reality. In the case of chaotic scattering, most of the theoretical investigations so far have been restricted to Hamiltonian or conservative systems. In a realistic situation, a small amount of dissipation can be expected. Take, for example, chaotic scattering arising in the context of particle advection in hydrodynamical flows . In most existing studies, the condition of incompressibility is assumed of the underlying flow , which allows the problem to be casted in the context of Hamiltonian dynamics as the particle velocities can be related to flow’s stream function in a way that is completely analogous to the Hamilton’s equations in classical mechanics. Real hydrodynamical flows cannot be perfectly incompressible, and the effects of inertia and finite mass of the particles advected by the flow are effectively those due to friction, or dissipation .

The aim of this paper is to study the effect of dissipation on chaotic scattering dynamics. We first consider hyperbolic chaotic scattering and argue that weak dissipations have a negligible effect on the physical observables of chaotic scattering, such as scattering functions. We then focus on nonhyperbolic chaotic scattering and find that, in contrast to the hyperbolic case, the scattering dynamics can be altered by weak dissipation in a fundamental way. The major consequence of dissipation is that it typically converts KAM tori into periodic attractors. As a result, the underlying chaotic saddle can undergo a metamorphic bifurcation to a structurally different chaotic set, playing the role of chaotic invariant set that generates fractal basin boundaries . There is an immediate transformation of the decay law of scattering particle from being algebraic in the Hamiltonian case to being exponential in the dissipative case, no matter how small the amount of dissipation. As a result, the fractal dimension of the chaotic saddle decreases from the integer value in the Hamiltonian case. These findings have striking implications to the study of chaotic scattering: they suggest that the algebraic-decay law, regarded to hold universally in nonhyperbolic chaotic scattering, is apparently structurally unstable against weak dissipations. More importantly, the previously believed integer dimensions of the chaotic saddles  in nonhyperbolic chaotic scattering may not be observable in realistic physical situations where dissipation is present.

We begin by presenting a picture for the formation of the fractal sets in chaotic scattering. Consider the idealized model of the hierarchical construction of Cantor sets in the unit interval. For hyperbolic scattering, an open subinterval in the middle of the unit interval is removed first. From each one of the two remaining subintervals, the same fraction from their middle is removed, and so on. Each step of this construction can be thought of as an iteration of the hyperbolic tent map with slope larger then two. The total length that remains decays exponentially with the number of iterations and the resulting Cantor set has a fractal dimension (e.g., box-counting dimension) smaller than one. For nonhyperbolic chaotic scattering, the same construction applies but the fraction removed at each step decreases with time, say, is inversely proportional to time. This simple reduction of the fraction removed captures the essence of the effect of KAM tori: their “stickiness” to particle trajectories in the phase space . The remaining length decays algebraically with time and, even though the measure of the remaining set asymptotes to zero, the resulting Cantor set has dimension one .

How does dissipation change the above construction? The skeleton of the underlying chaotic saddle is formed by periodic orbits. When the system is hyperbolic, its structural stability guarantees the survival of all periodic orbits under small changes of the system parameters. Accordingly, the structure of the Cantor set in the presence of a small amount of dissipation is expected to be the same as before. When the dynamics is nonhyperbolic, however, qualitatively different behavior can take place. Marginally stable periodic orbits in KAM islands can become stable, turning their nearby phase-space regions into the corresponding basins of attraction . This means that, part of the previous chaotic saddle now becomes part of the basins of the attractors. Most importantly for the scattering dynamics, the converted subset supports orbits in the neighborhood of the KAM islands that otherwise would be scattered after a long, algebraic time. These orbits are solely responsible for the nonhyperbolic character of the scattering in the conservative case. Due to the existence of dense orbits in the original chaotic saddle, the noncaptured part of the invariant set remains in the boundaries of basins of the periodic attractors. Therefore, the invariant set is the asymptotic limit of the boundaries between scattered and captured orbits, rather than those between scattered and scattered orbits as in the conservative case. Chaos thus occurs on the nonattracting invariant set whose stable manifold becomes the boundary separating the basins of the attractors and of the scattering trajectories. Through this simple reasoning we can see that the structure and the meaning of the Cantor set is fundamentally altered: in successive steps we remove a constant instead of a decreasing fraction in the middle of each interval. As a result, the scattering dynamics becomes hyperbolic with exponential decay. The dimension of the Cantor set immediately decreases from unity as a dissipation parameter is turned on. There is now more than one possible outcome: some of the removed intervals correspond to scattered orbits and the others correspond to orbits captured by the attractors. We stress that the appearance of attractors accompanied by a metamorphosis of the chaotic saddle can occur for arbitrarily small dissipation.

We now present numerical support for the effect of weak dissipation on chaotic scattering, particularly the metamorphic transformation of particle decay and fractal dimension in the nonhyperbolic case. Our model is a dissipative version of the two-dimensional area-preserving map utilized in Ref.  to establish the unity of the fractal dimension, a particularly convenient model for studying nonhyperbolic chaotic scattering. The map reads

$$M\left(\begin{array}{cc}x\hfill & \\ y\hfill & \end{array}\right)=\{\begin{array}{cc}\lambda [x(x+y)^2/4\nu (x+y)]\hfill & \\ \lambda ^1[y+(x+y)^2/4]\hfill & \end{array},$$
(1)
where $`\lambda >1`$ and $`\nu 0`$ are parameters. The map is conservative for $`\nu =0`$ and dissipative for $`\nu >0`$. For $`\nu =0`$, almost all orbits started from negative values of $`y`$ are scattered to infinity. In this case, the dynamics is nonhyperbolic for $`\lambda 6.5`$ and hyperbolic for $`\lambda 6.5`$. The computation of the particle decay and fractal dimension in nonhyperbolic scattering requires examining very small scales, which makes the numerical computation a highly nontrivial task . The advantage of using map (1) instead of a continuous flow is that it makes high-precision computation possible. Our results are, however, expected to hold in typical nonhyperbolic systems with KAM tori.

We study map (1) in the nonhyperbolic regime, without and with dissipation. We set $`\lambda =4.0`$. When there is no dissipation ($`\nu =0`$), there is a major KAM island in the phase space, as shown in Fig. 1(a). The fractal boundaries of the basins of scattering trajectories to infinity are also shown, which correspond to the stable manifold of the chaotic saddle in the scattering region. When dissipation is present ($`\nu >0`$), the fixed point in the center of the island becomes an attractor. Dynamically, it happens because the magnitudes of the eigenvalues of periodic orbits associated with islands are one, which are reduced by dissipation in general. The basin of attraction of this attractor “captures” the island itself and orbits close to the stable manifold of the previously existing invariant set, as shown in Fig. 1(b). The intricate character of the basin of attraction with apparent fractal boundaries comes from points of the invariant set that are arbitrarily close to the island for $`\nu =0`$. The newly created basin of attraction contains these points and hence, all their preimages as well. These preimages extend in the phase space along the original stable manifold of the chaotic saddle, which is the reason that the boundaries mimic those of the original basins of scattering trajectories \[Figs. 1(a) versus 1(b)\]. Because of this similarity, the scattering functions and time-delay functions, which are physically measurable, resemble each other in both the conservative and weakly dissipative case, as shown, respectively, in Figs. 1(c) and 1(d), where the time delay of particles launched from the horizontal line $`y=2`$ toward the scattering region is plotted against their $`x`$-coordinates on the line.

To examine the decay laws of the scattering particles, we approximate the survival probability of a particle in the scattering region by $`R(n)`$, the fraction of a large number of particles still remaining in the scattering region (defined by $`\sqrt{x^2+y^2}<r`$) at time $`n`$, which are initiated in subregions close to the boundaries of the scattering basins. For convenience, we choose $`r=100`$ and choose initial conditions from the horizontal line at $`y_0=2`$. When the dynamics is nonhyperbolic and conservative, the decay of $`R(n)`$ with time is exponential for small $`n`$ and algebraic for large $`n`$, as shown in Fig. 2(a) for $`\lambda =4.0`$: $`Re^{\alpha n}`$ for $`n250`$ and $`Rn^\beta `$ for $`n250`$, where $`\alpha 0.08`$ and $`\beta 1.0`$. In the presence of a small dissipation, the time decay becomes strictly exponential and with the same decay rate of the exponential regime of the conservative case, $`\alpha 0.08`$, as shown in Fig. 2(b) for $`\nu =0.001`$. The original algebraic decay in the conservative case is destroyed by the dissipation because orbits with points close to the island, and that otherwise would be stuck, are captured by the periodic attractor. The decay rate in general changes under further increases of the dissipation. For $`\nu =0.01`$, for instance, we obtain $`\alpha 0.06`$. In the hyperbolic region the time decay is always exponential. For $`\lambda =8.0`$, for example, the decay rate $`\alpha `$ remains essentially constant and equal to $`0.9`$ in the range $`0\nu 0.01`$.

The uncertainty algorithm  can be used to compute the fractal dimension $`D`$ of the set of intersection points between the stable manifold of the chaotic saddle and a line from which scattering particles are initiated. As above, we choose the line $`y_0=2`$. In the absence of attractors, $`D`$ is the dimension of the set of singularities in scattering functions. In Ref. , it is argued that $`D=1`$ when map (1) is nonhyperbolic and conservative. A technical point about the numerical evaluation of the dimension in this case is that the result converges slowly to unity, and the convergence rate is determined by the reduction of length scales in the computation . When a small amount of dissipation is present, $`D+1`$ becomes the dimension of the boundaries between scattered and captured basins. The numerical convergence of $`D`$ is in this case faster and essentially independent of the size of the interval under consideration. For $`\lambda =4.0`$ and $`\nu =0.01`$, we obtain $`D0.8`$, a well-convergent value as the length scale is reduced over six orders of magnitude. The dimension is much less sensitive to the presence of dissipation if the dynamics is hyperbolic. For $`\lambda =8.0`$, for instance, we obtain $`D0.44`$ for both $`\nu =0`$ and $`0.01`$.

Dissipation can also lead to several coexisting attractors for some values of $`\lambda `$ in the originally nonhyperbolic region. Periodic attractors are created through saddle-node bifurcations as $`\lambda `$ is varied. These attractors then undergo period doubling cascades until the accumulation point where chaos appears. For small dissipation, however, the chaotic interval in the parameter space can be so small that it is difficult to detect chaotic attractors numerically. In fact, for $`\nu `$ on the order of $`0.01`$ or smaller, the dynamics of map (1) is dominated by low-period periodic attractors. Periodic attractors of high periods either have small basins of attraction or exist in small intervals in the parameter space. Since periodic attractors result from the stabilization of periodic orbits in KAM islands, the parameter regions in which these attractors exist are approximately the same as those of the corresponding islands. In addition, the sizes of the basins are of the same order of the sizes of the original islands in the phase space .

In summary, our qualitative and quantitative examinations indicate that weak dissipation, no matter how small, can fundamentally alter the nature and dynamics of nonhyperbolic chaotic scattering. The algebraic-decay law, commonly believed to hold in such a case, is typically converted into an exponential-decay law in the metamorphic sense that the conversion can be induced by arbitrarily small amount of dissipation. A consequence of such a metamorphosis is that the previously claimed  unity of the fractal dimension of the set of singularities in scattering functions may not be physically meaningful. To our knowledge, there has been no previous attempt to address the effect of dissipation in open Hamiltonian systems, but this is a physically important issue of nonlinear dynamics that deserves further attention.

This work was supported by AFOSR under Grant No. F49620-98-1-0400. AEM acknowledges financial support from FAPESP.


format: text

text:
# A Model for Granular Texture with Steric Exclusion

## Abstract

We propose a new method to characterize the geometrical texture of a granular packing at the particle scale including the steric hindrance effect. This method is based on the assumption of a maximum disorder (entropy) compatible both with strain-induced anisotropy of the contact network and steric exclusions. We show that the predicted statistics for the local configurations is in a fairly agreement with our numerical data.

In modeling granular media from a microscopic point of view, one is naturally led to consider a fundamental question: how do the properties of the spatial organization of particles control the behaviour at the macroscopic scale? A major difficulty is that the effects arising from the geometrical frustration of particles prevail due to hard core exclusions, and thus the standard methods developped for multibody systems fail as such to provide a satisfactory description of many subtle structural features of a granular medium.

The most basic effect of mutual exclusions is that a particle may be simultaneously touched by only a few neighbouring particles. The coordination number (the number of touching neighbours) cannot exceed 6 in a 2D assembly of particles of nearly the same size and 12 in a 3D packing. These local environments fluctuate strongly in space both in terms of the coordination number and the angular positions of the neighbours.

Another basic observation is that the relative angular positions of the particles i.e. the directions of the contact normals are not isotropically distributed .

The adequate description of local environments including these geometrical features, i.e. a combination of excluded-volume constraints, disorder and anisotropy, is a crucial step in tracing back the failure and flow properties of granular materials to the particle scale. Consider, for instance, a granular packing in statistical equilibrium. The force balance on a particle involves the angular positions of the few neighbours belonging to the local environment. The fact that two neighbours may not occupy two angular positions separated by an angle smaller than a finite angle $`\delta \theta `$ (approximately $`\frac{\pi }{3}`$ for particles of nearly the same size, see figure 1) reduces drastically the accessible equilibrium states of the particle and hence that of the whole packing. This effect is even more critical at incipient failure and for larger anisotropies of the contact directions. While the texture anisotropy has been considered by several authors in the past , the key role of steric exclusions has been almost completely disregarded in micro-macro models of the quasi-static rheology of granular materials.

In this paper, we propose a model that allows to construct local environments incorporating both the steric exclusions and generic disorder of granular media from the global geometrical texture described in terms of the directions of interparticular contacts or the corresponding anisotropy. The basic assumption of our model is a maximum directional disorder consistent both with the texture anisotropy and with steric constraints. We apply our approach to a two-dimensional medium and solve the resulting equation. The predicted statistics of local environments will then be compared with raw data from numerical simulations.

Let us consider the simple model of a 2D granular assembly of rigid disks. In simulations or experiments in 2D geometry, a weak polydispersity is necessary to avoid local crystalline order. In our theoretical description, as we shall see below, the geometrical disorder is explicitly taken into account through the coordination number and the contact directions. We may thus assume that particles are of the same size, and the steric hindrance will be characterized by a single exclusion angle $`\delta \theta =\pi /3`$.

The directional organization of the contact network is described by the probability density of contact normals $`p(\theta )`$. This is the probability that a given particle has a contact along the direction parametrized by $`\theta `$, the polar angle of the contact normal $`\stackrel{}{n}`$; figure 1. This distribution is induced by the relative motions of particles, and so it essentially reflects the deformation history or the dynamics of the preparation process. The function $`p(\theta )`$ is $`\pi `$-periodic. It is often assumed that it can be approximated by its truncated Fourier transform

$$p(\theta )=\frac{1}{\pi }\left\{1+a\mathrm{cos}2(\theta \theta _p)\right\},$$
(1)
where $`a`$ represents the anisotropy of the texture and $`\theta _p`$ is the average contact direction. This functional form is reasonable for a simple deformation history, but it may considerably deviate from this simple form otherwise . In order to bypass a particular fitting form such as equation 1 for characterizing the anisotropy, one may more generally construct the “fabric tensor” $`𝑭\stackrel{}{n}\stackrel{}{n}`$ where the brackets denote averaging over all contacts in a representative element of volume, and $``$ is the dyadic product. Then, the texture anisotropy $`a`$ is defined to be $`a=2|F_1F_2|`$ where $`F_1`$ and $`F_2`$ are the eigenvalues of $`𝑭`$.

A richer description of the texture is provided by the probability density $`f_z(\theta )`$ that a particle with $`z`$ contacts has one contact along $`\theta `$. This function is defined over the range $`[0,2\pi ]`$. Like $`p(\theta )`$, it can be Fourier expanded, but the corresponding anisotropy coefficients now will depend on $`z`$. At leading order in $`\theta `$, we have

$$f_z(\theta )=\frac{1}{2\pi }\left\{1+a_z\mathrm{cos}2(\theta \theta _z)\right\}.$$
(2)
It is important to point out here that the steric exclusions do not allow for arbitrary global distributions $`f_z`$. For $`f_z`$ defined by Equation 2, the anisotropy $`a_z`$ could take any values in the range $`[0,1]`$ ($`a_z1`$ implies negative values of $`f_z`$ for some directions, whereas negative values can be avoided by turning $`\theta _z`$ into $`\theta _z+\pi /2`$). Nevertheless, the steric exclusions impose an upper bound $`a_{max}`$ on the anisotropy. Indeed, at most one contact can be found within an angular sector of $`\pi /3`$ around a particle, i.e.

$$_{\theta \frac{\pi }{6}}^{\theta +\frac{\pi }{6}}zf_z(\theta )𝑑\theta 1$$
(3)
which together with the expression of $`f_z`$ given Eq. (2) yields

$$a_{max}(z)=\frac{4\pi }{\sqrt{3}}\left(\frac{1}{z}\frac{1}{6}\right).$$
(4)
This means that $`a_{max}`$ decreases with $`z`$ and becomes zero for $`z=6`$. In order to evaluate $`f_z`$ for a granular sample, the subset of particles with $`z`$ contacts are to be singled out. Our numerical simulations show that $`a_z`$ indeed decreases with $`z`$.

None of the above functions however directly provides an exhaustive representation of the local environment of a particle. A complete description requires the coordination number $`z`$ of the particle and the directions $`\theta _1,\mathrm{},\theta _z`$ of the contact normals around it. Those angles should fulfill the constraint of steric exclusion around particles, so that an interval of width $`\pi /3`$ centered on the normal direction of any contact is excluded. In this way, the particle environments should be described through multicontact probability distributions $`g_z(\theta _1,\mathrm{},\theta _z)`$. This distribution is $`2\pi `$-periodic and its integration over all contact directions but one, $`\theta _k`$, gives back the global distribution $`f_z(\theta _k)`$ defined above:

$$_{_z(k)}g_z(\{\theta _i\})d\{\theta _i\}_{ik}=f_z(\theta _k),$$
(5)
where $`_z(k)`$ is the domain corresponding to integration from $`0`$ to $`2\pi `$ over all the $`\theta _i`$ except $`\theta _k`$.

Here, we would like to construct the local distributions $`g_z(\{\theta _i\})`$ from the global distribution $`f_z(\theta )`$. Of course, the solution is not unique since $`g_z`$ contains a much richer information than $`f_z`$. The point, however, is to get a solution that incorportes the required local information with no bias towards a particular solution. In the absence of local steric constraints, the most unbiased situation implies that

$$g_z(\theta _1,\mathrm{},\theta _z)=\underset{i=1}{\overset{z}{}}f_z(\theta _i).$$
(6)
In our case, this solution is wrong since the steric exclusions require

$$g_z(\{\theta _i\})=0\mathrm{if}|\theta _i\theta _j|<\pi /3,$$
(7)
for $`ij`$ which is obviously not satisfied by the above solution.

The solution that we propose is still to resort to a similar least biased estimate, provided the above steric constraints are taken into account. Operationnally, this translates into the maximization of an entropy functional $`S[g_z]`$ under constraints. According to Shannon’s entropy, we have

$$S[g_z]=_{𝒟_z}g_z(\{\theta _i\})\mathrm{log}[g_z(\{\theta _i\})]d\{\theta _i\},$$
(8)
where $`𝒟_z=([0,2\pi ])^z`$ is the integration domain. This functional should be maximized over the set of functions that fulfill both the steric exclusions, Eq. (7), and the normalization conditions, Eq. (5). We note that the entropy formalism, as a basic statistical tool, has been previously applied to granular materials in order to characterize other aspects of a random packing such as the coordination number and the void ratio .

Equation (7) implies that the maximization of $`S`$ has to be restricted to the region $`𝒜_z`$ of $`𝒟_z`$ that is allowed by steric exclusions. Indeed, the maximization applies only to the lacking information, and here the value of $`g_z`$ is already known ($`g_z=0`$) for the angles belonging to the excluded region $`_z=𝒟_z𝒜_z`$. Hence, the above restriction allows to account for the steric constraints.

On the other hand, the normalization conditions (Eq 5.) can be imposed through Lagrange multipliers. This leads to the maximization of the functional

$$T[g_z]=S[g_z]\underset{i=1}{\overset{z}{}}_0^{2\pi }\lambda _i(\theta _i)C_i[g_z]𝑑\theta _i$$
(9)
over $`𝒜_z`$, where the $`\lambda _i(\theta _i)`$ are the Lagrange multipliers and

$$C_i[g_z]=f_z(\theta _i)_{_z(i)}g_z(\{\theta _i\})d\{\theta _k\}_{ki}.$$
(10)
The solution is determined by setting the functional derivative of $`T[g_z]`$ to zero. Let us remark that here the only correlations among the angles are the steric exclusions themselves. This means that the angles $`\theta _i`$ can be considered as independant variables when they are restricted to the allowed domain $`A_z`$. As a consequence, the solution is given by a product of identical functions $`h`$ of a single angle in the allowed domain. In this way, the general solution is written as follows

$$g_z(\{\theta _i\})=\left(\underset{ij}{}G(\theta _i\theta _j)\right)\underset{k=1}{\overset{z}{}}h_z(\theta _k),$$
(11)
where $`G(\theta )`$ is a $`2\pi `$-periodic function such that $`G(\theta )=0`$ for $`|\theta |<\pi /3`$, and $`G=1`$ otherwise. These prefactors take care of the steric exclusions. Note that without these exclusions, we have $`G=1`$ for all angles and the solution (6) is recovered with $`h_z(\theta )f_z(\theta )`$.

The function $`h_z`$ is determined by the normalization condition (5) which becomes now the following implicit equation in terms of $`h`$:

$$_{𝒜_z(i)}\underset{k=1}{\overset{z}{}}h_z(\theta _k)d\{\theta _j\}_{ji}=f_z(\theta _i),$$
(12)
where $`𝒜_z(i)=_z(i)𝒜_z`$. This equation can be used for a numerical computation of $`h_z`$ for a given $`f_z`$. We used a simple iterative scheme over the functions $`u_n(\theta )`$ satisfying

$$u_{n+1}(\theta _1)=ϵu_n+(1ϵ)\frac{f_z(\theta _1)}{_{A_z(1)}u_n(\theta _2)\mathrm{}u_n(\theta _z)𝑑\theta _2\mathrm{}𝑑\theta _z}.$$
(13)
where $`ϵ`$ is a relaxation parameter.

The solution $`h_z`$ of equation (12) is simply a fixed point of the latter recursion. Our iterative process converges smoothly to the solution when the latter exists and, moreover, its convergence domain matches exactly the set of admissibles values of the anisotropy according to equation (4). We show one example of solution in Figure 2: The function $`h_3(\theta )`$ was computed according to the above procedure with $`f_3(\theta )`$ given by Eq. (2) for $`z=3`$ and $`a_3=0.3`$. We see that $`h_3`$ has the same monotony as $`f_3`$, but it deviates from a simple sinusoidal form.

Of course, it is desirable to compare the multicontact distributions $`g_z`$ predicted by the above approach with direct data from experimental or numerical observations. However, for a sufficient statistical precision requires very large samples. In fact, given a granular sample, the subset of particles with $`z`$ contacts are considered. Assuming that in a 2D system of disks, the number of particles is the same for each of the prevailing values 3, 4 and 5 of $`z`$, which is actually not the case since there are typically more particles with four contacts, then only nearly one third of the particles is available for the evaluation of $`g_z`$. On the other hand, subdividing the interval $`[0,2\pi ]`$ into $`n`$ angular sectors and requiring on average $`m`$ events in each elementary box in $`[0,2\pi ]^z`$, one would need a sample of $`m(n^z)`$ contacts, or approximately $`(2m/z)n^z`$ particles. For $`z=4`$, $`n=20`$ and $`m=100`$, this amounts to $`8.10^6`$ particles, which is practically inaccessible with the present-time computation power.

Nervertheless, we still may get a reliable comparison by focussing on the distributions $`p_d(\mathrm{\Delta }\theta )`$ of the difference $`\mathrm{\Delta }\theta =|\theta _1\theta _2|`$ between contact directions extracted from pair distributions $`p_2(\theta _1,\theta _2)`$. The latter was calculated from numerical samples and from the theoretical estimation of $`g_z`$. To do so, we performed numerical simulations of a system of 4000 particles with bi-periodic boundary conditions by means of the molecular dynamics method. We used a viscous-regularized Coulomb friction law and a linear spring-dashpot model for particle interactions . The coefficient of friction was $`0.5`$ and the particle radii were uniformly distributed between $`R_1`$ and $`R_2`$ with $`R_2=2R_1`$. The system was subjected to simple shearing and thanks to bi-periodic boundary conditions we were able to keep shearing in the steady state for very large strains without any wall effect. Figure 3 shows the function $`f_4(\theta )`$ corresponding to the subset of particles with four contacts. The global steady-state anisotropy, $`a_z`$, calculated from the fabric tensor was about $`0.13`$.

The distribution $`p_d`$ was calculated over the cumulated data from several well-separated snapshots of the steady state for the set of particles with coordination number $`z`$. For the theoretical evaluation of $`p_d`$ corresponding to these particles, we used the numerical distributions $`f_z(\theta )`$ for the calcultation of $`g_z`$ following our approach, from which the distributions $`p_2`$ and $`p_d`$ were extracted.

In figures 4 and 5, we have plotted theoretical and numerical $`p_d`$ as a function of $`\mathrm{\Delta }\theta `$ for $`z=4`$ and $`z=5`$, respectively. Although only the geometrical constraints are taken into account for the theoretical construction of $`g_z`$ (no force considerations), we see that the predicted distributions $`p_d`$ both for $`z=4`$ and $`z=5`$ fit fairly well the numerical data. The fit is nevertheless better for $`z=5`$ than for $`z=4`$. On the other hand, the fit for $`z=3`$ (not shown) is much less satisfactory. In fact, in going from $`z=3`$ to $`z=5`$, the requirement of force balance becomes less and less stringent, whereas steric exclusions dominate more and more the behaviour. Let us also remark that for $`z=4`$ a small peak appears on the numerical curve at $`\mathrm{\Delta }\theta =2\pi /3`$, as can be seen in figure 4. This peak, which is absent from the theoretical curve, is probably the signature of a local order induced by the rather small polydispersity of our sample.

These results suggest that, if the steric exclusions are correctly taken into account, a satisfactory evaluation of the statistics of local environments in terms of the functions $`g_z`$ can be obtained at least for $`z=4`$ and $`z=5`$. The steric exclusions control to a large extent the local distributions (figures 4, 5) and reduce the range of admissible anisotropies (equation 4) .

A natural extension of this model, at the cost of additional parameters, is to include the force balance as well as the steric exclusions. This would allow to tackle the case of $`z=3`$ more accurately than we did in this paper. The idea is to include both contact directions and contact forces in the model. The forces are coupled to a local stress tensor in the same way as the contact directions were coupled to the fabric tensor. Then, the entropy of local distributions is maximized regarding constraints arising from steric exclusions, but also the force balance, and constitutive laws such as Coulomb friction and Signorini’s condition (positivity of normal forces).

Apart from its importance in itself as a general model for the granular texture including steric exclusions, the characterization of the texture in terms of local environments is a crucial step in order to relate the macroscopic behaviour to the particle scale. For example, the range of admissible stresses can be studied by considering representative configurations generated according to the multicontact statistics of local environments . A similar procedure based on entropy formalism may be applied as well to other relevant microscopic configurations such as the cells composed of contiguous particles that carry the most local strains at the microscopic strains.


text:
# Zinc Impurities in the 2D Hubbard model

## Introduction

Chemical substitution in high-$`T_c`$ superconductors provides a powerful probe into the complex nature of both the superconducting and normal state of these materials. Experiments substituting different impurities for Cu show that nonmagnetic impurities (Zn, Al) are just as effective in suppressing superconductivity as magnetic dopants (Ni, Fe) Xiao . Based on Anderson’s theorem And  these results are taken as strong evidence for an unconventional pairing state described by an anisotropic order parameter with nodes on the Fermi surface. Indeed, recently the characteristic fourfold symmetry of the $`d_{x^2y^2}`$-wave order parameter was observed in the spatial variation of the local density of states near the Zn impurity using scanning tunneling microscopy Pan .

Nuclear magnetic resonance experiments show that even a nonmagnetic impurity substituted for Cu induces an effective magnetic moment residing on the neighboring Cu sites NMR . In addition, the bulk susceptibility in impurity substituted underdoped cuprate superconductors shows Curie-Weiss like behavior irrespective of the magnetic structure of the dopants Xiao . This indicates the existence of free magnetic moments in the impurity doped $`\mathrm{CuO}_2`$ planes. The formation of these moments with substitution of nonmagnetic impurities is explained to arise from breaking singletsXiao ; Finkel  in an antiferromagnetically correlated host by removing Cu spins.

In this letter we focus on the suppression of superconductivity by Zn impurities and the change of the bulk magnetic susceptibility to a Curie-Weiss like behavior in underdoped systems. Despite a wide variety of theoretical studies addressing these issues, a complete understanding of the effects of Zn doping has not been achieved. These approaches are mostly based on the description of impurities embedded in a BCS host BCS\_Zn  or on generalizations of the Abrikosov-Gorkov equations for nonmagnetic impurities in unconventional superconductors using phenomenological pairing interactions AG\_Zn . In order to capture effects like moment formation by substitution of nonmagnetic impurities, it seems necessary to describe correlations and the scattering from impurities on the same footing within a microscopic approach.

The most widely applied technique developed to describe disordered systems is the Coherent Potential Approximation (CPA) CPA . The CPA shares the same microscopic definition DCA\_Jarrell1  as its equivalent for correlated clean systems, the Dynamical Mean Field Approximation (DMFA) DMFA . Both approaches map the lattice system onto an effective impurity problem embedded in a host that represents the remaining degrees of freedom. This single-site approximation neglects interference effects of the scattering off different impurity sites (crossing diagrams) and correlation effects become purely local. Therefore, it inhibits a transition to a state described by a nonlocal ($`d`$-wave) order parameter and thus is not appropriate for our investigations here.

The Dynamical Cluster Approximation (DCA) DCA\_Hettler1 ; DCA\_Hettler2 ; DCA\_Maier1 ; DCA\_Jarrell1  systematically incorporates nonlocal corrections to the CPA/DMFA by mapping the lattice system onto an embedded periodic cluster of size $`N_c`$. For $`N_c=1`$ the DCA is equivalent to the CPA/DMFA and by increasing the cluster size $`N_c`$ the length-scale of possible dynamical correlations can be gradually increased while the DCA solution remains in the thermodynamic limit. In the clean limit the DCA applied to the Hubbard model has been shown to describe the essential low energy physics of the cuprates DCA\_Jarrell2 ; DCA\_Maier2 ; DCA\_Jarrell3 : It captures the antiferromagnetic phase near half filling and the transition to a superconducting phase with $`d_{x^2y^2}`$-wave order parameter at finite doping. In the normal state it exhibits non-Fermi liquid behavior in form of a pseudogap in the density of states and a suppression of spin excitations at low temperatures in the underdoped regime.

In this letter we study the two-dimensional (2D) Hubbard model including a potential scattering term according to the chemistry of Zn impurities using the DCA. We will show that potential scattering by Zn impurities strongly suppresses superconductivity as well as spin fluctuations and changes the magnetic susceptibility to a Curie-Weiss like temperature dependence in the underdoped region.

## Formalism

The correlated electrons in a $`\mathrm{CuO}_2`$ plane doped with Zn impurities are described by the 2D Hubbard model with diagonal disorder

$$H=t\underset{<ij>,\sigma }{}c_{i\sigma }^{}c_{j\sigma }+U\underset{i}{}n_in_i+\underset{i\sigma }{}ϵ_in_{i\sigma },$$
(1)
where we used standard notation. The disorder induced by the Zn sites occurs in the local orbital energies $`ϵ_i`$ which are independent quenched random variables distributed according to some specified distribution $`P(ϵ_1,\mathrm{},ϵ_N)=_{i=1}^NP(ϵ_i)`$ where $`N`$ is the number of lattice sites. For a concentration $`x`$ of Zn impurities we use the binary alloy distribution

$$P(ϵ_i)=x\delta (ϵ_i+V/2)+(1x)\delta (ϵ_iV/2),$$
(2)
where $`V`$ is the energy difference between the Cu and Zn 3$`d`$ orbitals. For simplicity we use site-independent values for the hopping integral $`t`$ and the Coulomb repulsion $`U`$.

The microscopic derivation of the DCA algorithm was discussed in detail for correlated systems in DCA\_Hettler2 ; DCA\_Jarrell3  and for disordered systems in DCA\_Jarrell1 . The DCA has a simple physical interpretation for systems where the intersite correlations have only short spatial range. The corresponding self-energy may then be calculated on a coarse grid of $`N_c=L_c^D`$ selected $`𝐊`$ points only, where $`L_c`$ is the linear dimension of the cluster of $`𝐊`$ points. Knowledge of the momentum dependence on a finer grid may be discarded to reduce the complexity of the problem. To this end the first Brillouin zone is divided into $`N_c`$ cells of size $`(2\pi /L_c)^D`$ around the cluster momenta $`𝐊`$. The propagators used to form the self-energy are coarse grained or averaged over the momenta $`𝐊+\stackrel{~}{𝐤}`$ surrounding the cluster momentum $`𝐊`$

$`\overline{G}(𝐊,i\omega _n)`$ $`=`$ $`{\displaystyle \frac{N_c}{N}}{\displaystyle \underset{\stackrel{~}{𝐤}}{}}G(𝐊+\stackrel{~}{𝐤},i\omega _n)`$ (3)
$`=`$ $`{\displaystyle \frac{N_c}{N}}{\displaystyle \underset{\stackrel{~}{𝐤}}{}}{\displaystyle \frac{1}{i\omega _nϵ_{𝐊+\stackrel{~}{𝐤}}\mathrm{\Sigma }(𝐊,i\omega _n)}}.`$
In Eq.(3) we make use of the fact that the approximation of the lattice self-energy $`\mathrm{\Sigma }(𝐤)`$ by the cluster self-energy $`\mathrm{\Sigma }(𝐊)`$ optimizes the free energy DCA\_Hettler2 ; DCA\_Maier1 . $`\mathrm{\Sigma }(𝐊)`$ summarizes the single-particle effects of the interaction term (second term) and the disorder term (third term) of the Hamiltonian, Eq.(1). The dispersion $`ϵ_{𝐊+\stackrel{~}{𝐤}}`$ denotes the spectrum of the kinetic part (first term) of Eq.(1). In order to avoid overcounting of self-energy diagrams on the cluster, the cluster excluded propagator

$$𝒢(𝐊)^1=\overline{G}(𝐊)^1+\mathrm{\Sigma }(𝐊)$$
(4)
is used as bare propagator for the cluster problem. To diagonalize the disorder part it is convenient to perform a Fourier transform to cluster real space $`𝒢_{ij}=1/N_c_𝐊e^{i𝐊(𝐗_i𝐗_j)}𝒢(𝐊)`$ where the disorder is diagonal. The inverse bare cluster propagator for a particular disorder configuration $`\{ϵ_i\}`$ is then given by

$$[𝒢_{ϵ_1,\mathrm{},ϵ_N}^1]_{ij}=[𝒢^1]_{ij}ϵ_i\delta _{ij},$$
(5)
and is used to initialize a QMC simulation to calculate the effects of the Coulomb interaction $`U`$. The QMC result for the cluster Green function $`G_{ϵ_1,\mathrm{},ϵ_N,ij}^c=G_{ij}^c[U,𝒢_{ϵ_1,\mathrm{},ϵ_N,ij}]`$ thus depends on the particular disorder configuration $`\{ϵ_i\}`$. The disorder averaged cluster Green function is then obtained from the individual results by

$$G_{ij}^c=G_{ϵ_1,\mathrm{},ϵ_N,ij}^c,$$
(6)
where the average $`<\mathrm{}>=_{i=1}^N\mathrm{d}ϵ_iP(ϵ_i)(\mathrm{})`$ is to be taken for a system of $`N_c`$ sites. Note that in principal QMC calculations for all possible disorder configurations $`\{ϵ_i\}`$ have to be carried through. However, contributions from configurations with $`m`$ Zn impurities on the cluster ($`mN_c`$) are weighted by a factor of $`x^m(1x)^{(N_cm)}`$ to the integral Eq.(6). Therefore it seems reasonable for small concentrations $`x1`$ and cluster sizes $`N_c`$ to consider only those configurations with none or only a single Zn impurity on the cluster and to neglect configurations with more than one Zn impurity. Then the disorder averaged cluster Green function $`G_{ij}^c`$ is obtained as a weighted sum from the two configurations as

$$G_{ij}^c=xN_cG_{1,ij}^c+(1xN_c)G_{0,ij}^c,$$
(7)
where $`G_{1/0,ij}^c`$ denotes the QMC result for the cluster Green function for the configuration with one or zero Zn impurities on the cluster, respectively. The prefactors follow with $`(1x)^{N_c}1xN_c`$, $`xN_c(1x)^{N_c1}xN_c`$ for $`x<<1`$ in linear approximation in the impurity concentration $`x`$. The disorder averaged cluster Green function $`G_{ij}^c`$ is then transformed back to cluster reciprocal space to calculate a new estimate of the cluster self-energy

$$\mathrm{\Sigma }(𝐊)=𝒢(𝐊)^1G^c(𝐊)^1,$$
(8)
which is then used to repeat the steps starting from Eq.(3) until convergence is reached.

## Results

We study the superconducting instability and magnetic properties of the Hamiltonian Eq.(1) in the intermediate coupling regime at 5% hole doping ($`n=0.95`$) for different Zn impurity concentrations $`x`$. We perform calculations for $`N_c=4`$, the smallest cluster size that allows for a transition to a $`d_{x^2y^2}`$-wave superconducting state while preserving the lattice translational and point group symmetries. We set the hopping integral $`t=0.25\mathrm{eV}`$ and the Coulomb repulsion $`U=W=2\mathrm{e}\mathrm{V}`$, where $`W=8t`$ is the bandwidth of the non-interacting system. We choose $`V=5\mathrm{e}\mathrm{V}>W+U`$ to simulate the nonmagnetic closed shell ($`d^{10}`$) configuration of the impurity site.

We searched for superconductivity with $`s`$-wave, extended $`s`$-wave, $`p`$\- and $`d`$-wave order parameters. As we found for the clean system ($`x`$=0) DCA\_Jarrell2 ; DCA\_Jarrell3 , only the $`d_{x^2y^2}`$-wave pair field susceptibility for zero center of mass momentum diverges at low temperatures. Fig.1 illustrates the inverse pair field susceptibility $`P_d^1`$ as a function of temperature $`T`$ for the impurity concentrations $`x=0`$ (circles), $`x=0.05`$ (squares) and $`x=0.10`$ (diamonds).

The corresponding critical temperature $`T_c`$ is then calculated by extrapolating $`P_d^1(T)`$ to zero using the function $`P_d^1(TT_c)^\gamma `$ (represented by the solid lines in Fig.1). $`T_c`$ is rapidly suppressed, and the critical exponent $`\gamma `$ increases, with increasing Zn content $`x`$. Consistent with experiments Xiao  this falloff of $`T_c`$ is linear in $`x`$ resulting in a critical impurity concentration $`x_c0.10`$, beyond which the instability to the superconducting phase disappears. Compared to experiments ($`x_c0.05`$) our calculation predicts a result about twice as high for the critical doping. However, considering the simplified description of the problem using a 2D Hubbard model with the complexity of the Zn impurities reduced to a model of diagonal disorder, we expect agreement only on a qualitative level.

In order to study the correlation between the suppression of superconductivity and possible moment formation by Zn substitution we study the bulk ($`𝐪=(0,0)`$) and the antiferromagnetic ($`𝐪=(\pi ,\pi )`$) spin susceptibilities. Fig.2 shows the bulk magnetic susceptibility $`\chi (T)`$ as a function of temperature for the same impurity concentrations $`x`$ used in Fig.1.

As we demonstrated in DCA\_Jarrell2 , the clean system exhibits an anomaly in $`\chi (T)`$ at a characteristic temperature $`T^{}`$, below which spin excitations are suppressed and $`\chi (T)`$ falls with decreasing temperature. $`T^{}`$ thus marks the onset of singlet formation due to strong short-ranged antiferromagnetic spin correlations. In agreement with experiments Xiao , the substitution of impurities leads to a dramatic change of the low temperature behavior of the magnetic spin susceptibility. Despite the nonmagnetic nature of the impurity scattering term in our simulation, $`\chi (T)`$ displays an upturn at low temperatures as $`x`$ is increased. Moreover at $`x=0.10`$, $`\chi (T)`$ shows Curie-Weiss like behavior; i.e. it can be fit by the function $`\chi (T)=\chi _0\chi _1T+C/(T\mathrm{\Theta })`$, with $`\mathrm{\Theta }0`$. This clearly indicates the formation of free magnetic moments due to the breaking of singlet bonds by the substitution of nonmagnetic impurities.

The inverse antiferromagnetic spin susceptibility $`\chi _{AF}(T)`$ for the same parameter set as in Fig.2 is illustrated in Fig.3. For the clean system (x=0) 5% hole doping marks the critical doping, beyond which no transition to an antiferromagnetic state can be found DCA\_Jarrell2 . As can be seen from the extrapolation of the circles in Fig.3 the corresponding Neél temperature $`T_\mathrm{N}=0`$. As the concentration of impurities is increased we notice that for $`T0.05`$ $`\chi _{AF}`$ decreases with $`x`$. However at lower temperatures ($`T0.05`$) $`\chi _{AF}^1`$ for finite impurity concentration $`x`$ falls below the curve of the clean system and eventually goes through zero indicating a transition to the antiferromagnetic state. Therefore it can be inferred that a finite concentration of impurities, i.e. spin vacancies, acts to enhance antiferromagnetic spin correlations.

Our results for the site and time averaged local magnetic moment $`T\chi _{ii}`$ shown in the inset to Fig.3 support this conjecture: At low temperatures this quantity increases with impurity concentration $`x`$ indicating that local moments get stabilized. It is important to note that recent NMR measurements on YBCO Jul  show a similar enhancement of antiferromagnetic spin correlations around Zn impurities at low temperatures.

Similar results were also obtained in Bulut ; Dagotto  where an antiferromagnetic Heisenberg model with a spin vacancy was studied. In the latter report, the enhancement of antiferromagnetic spin correlations was ascribed to pruning of singlets by the nonmagnetic impurity in an resonating-valence-bond picture (RVB). In the RVB picture, the short-ranged order in the doped model is described as a collection of nearest-neighbor singlet bonds which fluctuate as a function of time and space. When a Zn impurity is introduced into the system, the singlet fluctuations on adjacent sites are suppressed since these sites now have one less neighbor to form short-ranged bonds with. The suppression of singlet fluctuations on the sites adjacent to the Zn dopant also enhances moment formation on these sites.

We find further evidence of the suppression of spin fluctuations by Zn impurities in the decrease of the antiferromagnetic exponent $`\gamma `$ towards the mean-field result of one with increasing impurity content $`x`$. Concomitant with this behavior the corresponding exponent $`\gamma `$ for the pair-field susceptibility (cf. Fig.1) increases towards one. From this observation we conclude that the suppression of superconductivity in Zn doped systems originates in the suppression of antiferromagnetic spin fluctuations that mediate pairing.

## Summary

In this letter we have used the 2D Hubbard model with binary diagonal disorder to study the effects of nonmagnetic (Zn) impurities on high-temperature superconductivity. We find that Zn impurities strongly suppress superconductivity. As evidenced by the increasing mean-field character of the pairfield and antiferromagnetic susceptibilities with Zn substitution, spin fluctuations that mediate pairing get suppressed. Consistent with experiments, $`T_c`$ decreases linearly with impurity concentration. With increasing Zn content we further find a change of the magnetic susceptibility to Curie-Weiss like behavior indicative of the existence of free magnetic moments.

Our results can be understood within the RVB picture. The nonmagnetic impurities break local singlets and thus generate unpaired spins. A concomitant suppression of antiferromagnetic spin fluctuations that mediate pairing results from the pruning of RVB states Dagotto . Consequently, superconductivity gets strongly suppressed with increasing Zn content.

## Acknowledgements

We acknowledge useful conversations with W. Putikka, D.J. Scalapino, Y. Wang and M. Vojta. This work was supported by NSF grant DMR-0073308. This research was supported in part by NSF cooperative agreement ACI-9619020 through computing resources provided by the National Partnership for Advanced Computational Infrastructure at the Pittsburgh Supercomputer Center.


text:
# Evidence for sharper than expected transition between metastable and unstable states

## 1 Introduction

Metastability has traditionally been studied by looking for droplets . In this framework density fluctuations are produced from the initially homogeneous matrix following a sudden quench from a stable equilibrium state to a state underneath the coexistence curve. After a time-lag a quasi-equilibrium is established and a size distribution for the fluctuations (droplets) develops, yielding droplets up to a critical size .

Behind this approach to metastability is a geometric interpretation, fostered by experimental observations of droplets . The experimentally observed droplets are, however, in size much larger than the critical droplets so that the evidence for nucleation is indirect . Also in computer simulations one can follow the development of the droplet . Close to the coexistence curve the droplet picture seems to hold to a good degree , while deeper into the two-phase region considerable doubt has been raised as to the validity of the nucleation theory .

To take the point of view of density fluctuations  is by no means the only possible. Here we want to follow a different point of view. What we want to consider are the fluctuations of a macroscopically available observable and the ensemble of these fluctuations . To be specific and to link to the simulation results that we present below, consider a non-conserved order parameter. We start out from a given equilibrium state and suddenly bring the system into a non-equilibrium state underneath the coexistence curve. The order parameter will follow this change. Again, after some time–lag the order parameter will fluctuate around some quasi-equilibrium value and suddenly change to the stable equilibrium state. We now consider the ensemble of such trajectories of the order-parameter. We avoid to take an average of the order-parameter in the non-equilibrium situation so as not wash out any of the abrupt changes that take place while the system seeks its way out of the quasi-equilibrium state into the stable equilibrium state. For each of the trajectories we define a lifetime and consider the distribution of the lifetimes and also consider the average lifetime! We show below that the average lifetime dramatically decreases at particular values depending on the imposed parameters giving rise to a possible dynamic definition of a spinodal.

So far the existence of a spinodal in the sense of a line of second order phase transition inside the two-phase region has been denied for systems with short range interactions. Only for systems with infinite range interaction (mean field theory) this concept has been thought to make sense . The concept is based on the existence of a free energy and is purely static. Here we side step this and look at metastability in a dynamical way and find a clear distinction between metastable and unstable states.

## 2 Methodology

The trajectories discussed above can be obtained naturally using the Ising model together with the Monte Carlo Method . The Monte Carlo method defines the transition probabilities between states starting from an initial state.

The Hamiltonian of the Ising model for a simple cubic lattice $`L^3`$ is defined by the

$`_{\text{Ising}}(s)=J{\displaystyle \underset{<ij>}{}}s_is_j+h{\displaystyle \underset{i}{}}s_i,s_i=\pm 1`$ (1)
where $`ij`$ are nearest-neighbour pairs of lattice sites. The exchange coupling $`J`$ is restricted in our case to be positive (ferromagnetic). $`s_i`$ is called a spin and the sum over all lattice sites is the magnetization $`M`$ (we define $`m`$ to be the magnetization per spin). $`h`$ is a dimensionless magnetic field.

The dynamics of the system is specified by the transition probabilities of a Markov chain that establishes a path through the available phase space. Here we use a Monte Carlo algorithm  to generate Markov Chains. We used the Metropolis transition probabilities

$$P(s_is_i^{})=\mathrm{min}\{1,\mathrm{exp}(\mathrm{\Delta })\}.$$
(2)
Time $`t`$ in this context is measured in Monte Carlo steps per spin. One Monte Carlo step (MCS) per lattice site, i.e. one sweep through the entire lattice, comprises one time unit. Neither magnetization nor energy are preserved in the model which makes possible to compute both quantities as a function of temperature, applied field and time.

Typically we started the simulation runs with a magnetization of $`1`$, i.e. in equilibrium and a predefined starting point for the random numbers for the Monte Carlo process. We then turn on the temperature and the applied field. This brings us instantly below the coexistence curve into the two-phase region. Due to the transitions induced by the Metropolis transition probabilities the magnetization (and the energy as well) develops in time and yields a trajectory of magnetization (energy) values

$$m(0),m(1),\mathrm{},m(n).$$
(3)
The magnetization was traced to a specific number of Monte Carlo steps. After the prescribed number of steps were reached, the system was brought back to the magnetization $`1`$ and a new quench performed with a new and different starting point for the random numbers. The statistics was gathered over these trajectories and quenches with each trajectory kept on file.

## 3 Relaxation paths and their statistics

Let us for the moment fix the temperature ($`T=0.59`$) and analyse the distribution of lifetimes. Figure 1 shows a typical relaxation path. After fluctuating around some value the magnetization drops sharply. At this point we say that the lifetime of the quasi-stable state has been reached. After the sharp drop the magnetization reached quickly its value in equilibrium on the other side of the phase diagram.

For every quench depth, i.e., applied magnetic field, we can compile a statistics on the lifetimes. This is shown in Figure 2. Plotted there is the probability for the occurrence of a specific lifetime as a function of the quench depth. Note first that the distribution of lifetimes gets narrower as we quench deeper into the two-phase region. Closer to the coexistence curve, i.e. for shallow quenched, the average lifetime increases but the width of the distribution broadens. This is of particular interest for experiments. There we also should expect a huge variety of lifetimes after a quench making it difficult to interpret the data. Even a quench close to the coexistence curve can lead to a very short lifetime and thus to an interpretation that the state may have been unstable. In passing we note that the distribution is not Poissonian.

We can plot the average lifetime, at fixed temperature and plot this as a function of the quench depth. This is shown in Figure 3. While the average lifetime of the states close to the coexistence curve must be very large, we expect the lifetime to drop as we quench deeper into the two-phase region. The plot shows that the average lifetime decreases rapidly at some "critical quench depth". The inset in Figure 3 shows that the average lifetime decreases exponentially fast close to that field and then saturates. This saturation depends on the precise nature of the Monte Carlo move but must always be present since it always takes a minimum number of steps to proceed from one state to another.

Clearly, the location of the "limit of stability"-line depends on the kind of Monte Carlo moves that we implement. But, the line will always be present for all local move algorithms. For all natural systems we expect this result to carry over. Indeed for a pure $`CO_2`$-system a similar result has been seen . There the stability line for states in the two-region was also much closer to the coexistence curve then expected.

## 4 Finite-Size Effects

Figure 4 shows the finite-size effects for the average lifetime. We have performed simulations with the linear system sizes of $`L=32`$, $`L=42`$ and $`L=96`$ for the temperature $`T/T_c=0.58`$. First we note that the pseudo spinodal shifts towards the coexistence with increasing system size. A simple extrapolation gives a 30% shift with respect to the $`L=32`$ result. We further note the decrease in the average lifetime for those quench depths that extend beyond the pseudo spinodal point. Similar finite-size effects have also been seen by Novotny et. al. .

The inset in Figure 4 shows the extrapolation of the value of the limit of stability to the infinite system size.

In Figure 5 we have compiled the results of our findings in a phase diagram. The new limit of stability is far away from the mean field spinodal curve.

## 5 Conclusions

Relaxation paths and their statistics yield a fresh outlook on metastability. Here the emphasis is on the dynamics of the single events and not on the average behaviour of the system. Our point of view is that averaging in the non-equilibrium situation can wash away or hide the phenomenon. Here we have shown that the statistics of the relaxation paths can lead to a sharp distinction between meta- and unstable states.

A surprising finding is that the distribution of the lifetime broadens significantly. Opposite to the intuition we found that close to coexistence curve the distribution is broad and very flat. Close to the limit of stability the distribution is very narrow and sharply peaked.

Clearly further work needs to be done, cementing the evidence for the sharp border, in the sense defined above. As pointed out above, simulations using different transition probabilities would be helpful to investigate the dependence of this dynamic phenomenon on the imposed dynamics, at least for those systems that do not have an intrinsic dynamics.

It would also be of interest to perform molecular dynamics simulations of systems using realistic interaction potential. There the question would also be how far the metastable region really extends and compare this to experimental evidence for nucleation or spinodal decomposition in the light of the relaxation trajectories.

It is also unclear at the moment how to extend the definition to systems with a conserved order parameter.

## 6 Acknowledgment

This work was funded by a grant from BASF AG and a grant from the BMBF. C.E.C would like to thank the DAAD and FAPERJ for financial support. We are also very grateful for discussions with K. Binder, E. Hädicke, B. Rathke, R. Strey and D. Stauffer.