text:
# Synchronization in the non-extensive Kuramoto model

## Abstract

We study the onset of synchronization in square lattices of limit cycle oscillators with long-range coupling by means of numerical simulations of the Kuramoto model. In this regime the critical coupling strength depends on the system size and interaction range reflecting the non extensive behavior of the system, but an adequate scaling removes the dependency and collapses the long-range synchronization curves with the one resulting from a system with uniform coupling.

05.45.Xt,05.70.Fh

Coupled limit-cycle oscillators have been extensively used to model the collective behavior of systems in physics , chemistry and biological sciences . One of the simplest models to describe those systems is the Kuramoto model for the time evolution for the phases of the oscillators. As shown in , it possible to reduce any system of $`N`$ weakly coupled limit cycle oscillators near the Hopf bifurcation, to the set of equations

$$\dot{\theta }_i=\omega _i+\underset{j}{\overset{N}{}}\mathrm{\Gamma }_{ij}(\theta _j\theta _i),(i=1,2,\mathrm{},N)$$
(1)
where $`\theta _i`$ and $`\omega _i`$ are the phase and natural frequency of the $`ith`$ oscillator respectively. The coupling term $`\mathrm{\Gamma }_{ij}(\theta _i\theta _j)`$ is a periodic function that may be obtained from the original equations of motion, regardless of the dimensionality of the oscillating field. The case of uniform, or mean field, coupling has been considered with detail by Kuramoto  for the simplest form of the coupling term,

$$\mathrm{\Gamma }_{ij}=N^1K\mathrm{sin}(\theta _j\theta _i),(K>0)$$
(2)
where the factor $`N^1`$ appears in order to ensure that the typical strength of the net local field experienced by each oscillator is independent of the total number $`N`$. Kuramoto considered the possibility of a collective oscillatory motion when the natural frequencies are picked from a given distribution $`g(\omega )`$, as would happen in a real system. As the coupling modifies the frequencies, acting as an external oscillatory force, it is necessary to define asymptotic frequency of the $`ith`$ oscillator after the transient regime as

$$\stackrel{~}{\omega }_i=\underset{T\mathrm{}}{lim}\frac{1}{T}\left[\theta _i(t+T)\theta _i(t)\right].$$
(3)
The existence of a collective behavior with phase and frequency entrainment depends on the coupling strength and on the choice of $`g(\omega )`$ . In what follows we consider only the case where $`g(\omega )`$ is Gaussian. Kuramoto characterized the degree of synchronization in the asymptotic regime by means of two order parameter like quantities: $`R`$ for the frequencies, and $`\sigma `$ for the phases. Frequency entrainment may be quantified simply by analyzing the formation of frequency clusters and defining

$$R=\underset{N\mathrm{}}{lim}\frac{N_s}{N},$$
(4)
where $`N_s`$ is the size of the largest cluster of oscillators with a common asymptotic frequency $`\stackrel{~}{\omega }`$. The complex quantity

$$\sigma =\frac{1}{N}\underset{j=1}{\overset{N}{}}\mathrm{exp}(i\theta _j)$$
(5)
is appropriate for the phase ordering as it measures the concentration of phases on a certain value. $`R`$ and $`\sigma `$ behave as equilibrium order parameters in the sense that, for $`K`$ below some critical value $`K_c`$, oscillations are incoherent leading to $`R`$ and $`\sigma `$=0. For strong enough coupling, or $`K>K_c`$, collective oscillation takes place and $`R,\sigma >0`$.

For the case of uniform coupling, phase and frequency entrainment appear simultaneously at the critical value

$$K_c=\frac{2}{\pi g(\omega _0)},$$
(6)
where $`\omega _0`$ is the average value of the distribution of natural frequencies .

The synchronization in lattices of oscillators with nearest neighbor coupling has also been studied with equilibrium statistical mechanics tools. The existence of a critical dimension for synchronization was questioned in  as macroscopic synchronization would require $`R`$ to be finite in the limit $`N\mathrm{}`$. This issue was clarified analytically in  and numerically in  and. The behavior of the synchronization curve $`R\times K`$ for different systems sizes in two dimensional lattices clearly showed that there is a limiting curve corresponding to the thermodynamical limit  and, although far from equilibrium, the system is extensive. As opposed to the mean field case, there is no need for correction factors in the coupling term.

Here we focus on systems with long range distance dependent coupling of the form

$$\mathrm{\Gamma }_{ij}=\frac{K}{r_{ij}^\alpha }\mathrm{sin}(\theta _j\theta _i),K>0,$$
(7)
where $`r_{ij}`$ is the distance between oscillators $`i`$ and $`j`$. The limit $`\alpha =0`$ corresponds to uniform coupling, and $`\alpha =\mathrm{}`$ to nearest neighbor coupling. From the results obtained in those limits, one can expect that there is a value of $`\alpha `$ below which the system is non-extensive due to the long range interactions. In Hamiltonian systems the limit of non-extensivity can be found by requiring the upper value of the energy to be finite as $`N\mathrm{}`$. Limit cycle oscillators are open systems and energy is a meaningless quantity so, instead, we establish the finiteness condition for $`(1/N)_i\dot{\theta }_i`$ in the asymptotic regime. This condition may be established by requiring that the upper value of the coupling term divided by $`N`$ must be finite, as in the mean field coupling case. With this, we end up with a condition identical to the one found for extensivity of the Ising model with distance dependent coupling , for $`\alpha <d`$ the system is non-extensive. For any value of $`d`$ and $`\alpha `$, it is possible to incorporate the dependence on $`\alpha `$ and $`N`$ of the coupling term into the value of $`K`$ with the definition of a scaled coupling strength

$$K^{}(\alpha ,N)=KN^{},$$
(8)
where

$$N^{}(\alpha ,N)=_1^{N^{1/d}}𝑑rr^{d1}r^\alpha =\frac{N^{1\alpha /d}1}{1\alpha /d}$$
(9)
For $`\alpha >d`$ (8) provides a correction for finite size effects, and for $`\alpha <d`$ it removes the dependency of the critical coupling strength on $`\alpha `$ and $`N`$. For the mean-field coupling ($`\alpha =0`$) and $`N^{}=N`$, for $`\alpha =d`$ $`N^{}`$ behaves like $`\mathrm{ln}N`$.

In order to numerically study the non-extensive region, we used two dimensional square lattices of coupled oscillators with natural frequencies sorted from the Gaussian distribution $`g(\omega )=(1/\sqrt{2\pi })\mathrm{exp}(\omega ^2/2)`$. Equation (1) with the coupling (7) was then integrated using the Euler method with d$`t=0.1`$. The first 3000 iterations where discarded as a transient, and after 6000 iterations past the transient, we calculated the asymptotic frequencies as defined in (3). We used the Hoshen-Kopelman method  to analyze the frequency pattern, two sites were considered to belong to same cluster if the frequency difference between them was less than $`2\pi /T`$, where $`T`$ is the time elapsed after the transient. $`R`$ was then associated to the largest cluster, and averaged over 10 initial conditions for the natural frequencies.

The non-extensivity is evident if we consider the behavior of $`R`$ as a function of $`K`$ for $`\alpha <d`$. For fixed $`\alpha `$, the critical value of $`K_c`$ decreases with increasing $`N`$ as can be seen in Fig. (1.a), as more terms contribute to the coupling. The shift in $`K_c`$ may be corrected with the use of (8). In fact, if we plot $`R`$ as a function of $`K^{}`$, there is no dependence on the system size, the curves collapse and a unique value of $`K_c^{}`$ appears. Figure (1.b) shows the behavior of $`R`$ as a function of $`K^{}`$ for fixed $`\alpha `$.

The dependence of $`K_c`$ on $`\alpha `$ is evident if we plot $`R`$ for fixed $`N`$ and different values of $`\alpha `$. Figure (2.a) shows the curves for $`\alpha =0.25`$, 0.50, 0.75, and 0 (mean-field coupling without the correction factor $`1/N`$). Consistently with dependence on system size, as $`\alpha `$ decreases, more terms are included in the coupling, leading to smaller values of $`K_c`$. Again, (8) fixes the dependence on $`\alpha `$ providing a perfect collapse even with the mean-field curve as can be seen in Figure (2.b). Notice that the value of $`K_c^{}`$ in Figures (1) and (2) is in good agreement with (6), our frequency distribution leads to $`K_c^{}=1.59`$.

In summary, we have extended the recently proposed equilibrium non-extensive scaling to a system of coupled limit-cycle oscillators far from equilibrium, by requiring the finiteness of the long-range coupling term $`K/r^\alpha `$. With this requirement it was possible to define a universal coupling strength $`K^{}`$, independent of $`\alpha `$ and system size, allowing the collapse of the curves $`R\times K^{}`$ for different sizes and exponents, including the case of uniform coupling. This result allows the study of more general coupling forms, for example to include lattice deformation effects, for any range of the coupling term.


format: text

text:
# Heterogeneous Interfacial Failure between Two Elastic Blocks

## I Introduction

The failure of interfaces under stress is a problem that has obvious important technological relevance. In addition, from a more fundamental point of view, this problem exhibits very interesting features. It is the aim of this paper to bring out some of these features by means of a numerical model based on a discretization of the original problem.

During more than a decade, failure processes in different contexts have caught the attention of the physics community. For a considerable longer period, the mechanics community has been involved in the study of such phenomena. In order to place the present study in its proper context, we need to go back to 1926 with the study of Peirce on what today is known as the global load sharing fiber bundle . This consists of $`N`$ parallel fibers, each with its own breaking threshold and connected in such a way that when a fiber fails, the load it was carrying would be dsitributed equally among all the surviving fibers. In 1945, Daniels published a very thorough study of this model, which today forms the starting point of any excursion into this field . The model has since these early days been generalized in many directions, one of which consists in replacing the “democratic” load-sharing rule by different local ones. One much studied variant is the local load sharing model, where the load on the failing fiber is distributed equally among the nearest surviving fibers . The average load–deformation characteristics of the global load sharing model was calculated by Daniels . Corresponding work on the local load sharing model may be found in Refs. . There has also been several studies of time-dependent phenomena in connection with the two variants of the fiber bundle model, see Ref. . (This paper in addition contains a very thorough review of the literature in this field.) There are also a number of studies “on the market” that may be placed between the two extremes that the global and local load sharing models consitute. Among them, we find the early study by Newman and Gabrielov , who constructed a hierarchically connected fiber bundle. Other work on hierarchical fiber bundle models may be found in Refs. .

Much work by the physics community has gone into studying network models, of which the fuse model is the most well known . This model consists of a network of electrical fuses where their burn-out thresholds have been drawn from some probability distribution. This model may be regarded as yet another generalization of the fiber bundle model, however, this time along the axis on which we find chains of fiber bundles . Among the several interesting questions that have been studied in connection with the fuse model, we mention the question of whether the breakdown process has the character of a second or first order phase transition . Central to this question is the question of the distribution of fuses that burn out simultaneously or — equivalently in the fiber bundle model — the number of fibers that fail simultaneously. This question was first raised and solved analytically in the context of the global load sharing fiber bundle  and then for the local load sharing model . The same question was first studied in connection with the fuse model in Ref. .

The particular problem we study here, elastic interfacial failure, has been addressed in the literature earlier by Delaplace et al. . The system consists of two elastic media that have been welded together, thus sharing a common interface. In general, the media can have different elastic constants. However, for the sake of simplicity and without loss of generality, we assume one of the media to be infinitely stiff while the other is elastic. We can view this simplification as an effective representation of the original system since it does not change the physics. Furthermore, we assume that the “soft” medium is uniform with respect to its elastic properties; it has the same elastic properties everywhere. However, the local strength of the glue — defined as the maximum local load it may sustain without failing — varies from point to point along the interface. This is the source of disorder in the system. In real systems this disorder would typically be correlated. In this first attack on the problem, we assume the disorder to be uncorrelated. Our main interest here is to understand how correlations develop due to the failure process. In Section II we describe the numerical model in detail.

The two joined media are subjected to a progressive uniform load perpendicular to the glue interface. Local failures will develop in the interface which changes the stress field on the remaining intact interface. These changes in the stress field will compete with the local strength of the glue to determine where the next failure happens. Sometimes, a local failure will occur due to the glue being particularly weak at that point on the interface, other times failure will occur due to enhancements in the local stress field. This competition leads to the development of spatial correlations both in the stress field and in the failure patterns, and in the failure process itself.

The two media can be pulled apart by controlling (fixing) either the applied force or the displacement. The displacement is defined as the change in the distance between two points, one in each medium, positioned far from the glued interface. Clearly, the line connecting these points is perpendicular to the average position of the interface. In our case, the pulling is accomplished by controlling the displacement. As the displacement is increased very slowly, glued points will fail. Sometimes the failed regions are very small, other times, the failed region is larger. Such events, when a large area fails “intantaneously” compared to the time scale at which the displacement is changed, are called bursts. One of the quantities of interest to us is the burst distribution  as the failure process evolves. We show in Section III A that this distribution follows a power law.

In Section III B, we investigate the scaling properties of the load and displacement of the system at the point when the failure process becomes unstable. This is the point at which any further increase of either load or displacement will lead to a catastrophic burst where all remaining glue fails. This point defines the strength of the interface, and the question we pose is how this scales with the system size.

We then investigate the geometrical properties of the failed regions at the point when catastrophic failure occurs in Section III C. We find that the area distribution of the failed regions follow a power law.

We present our conclusions and outlook for further work in Section IV.

## II Model

The system described in the Introduction is continuous: Two media, one elastic, the other infinitely stiff, are glued together thus forming a common interface. In order to treat this problem numerically, the continuum problem is replaced by a discrete one. We use a discrete model for this. In Section II A we describe the discrete model, while the numerical algorithms are discussed in Section II B.

### A Description of Model

We discretize the glued interface by replacing it with two two-dimensional square $`L\times L`$ lattices with periodic boundary conditions. The lower one represents the hard, stiff surface and the upper one the elastic surface. The nodes of the two lattices are matched (i.e. there is no relative lateral displacement). The glue is modelled by springs connecting opposing nodes in the two lattices. These harmonic springs all have the same spring constant (set to unity) but breaking thresholds randomly drawn from a uniform distribution between zero and one. As the two glued media are separated, the forces carried by the springs increases. When the force carried by a spring reaches its breaking threshold, it breaks irreversibly and the forces redistribute. The springs are thus broken one by one until the two media are no longer in mechanical contact. As this process is proceeding, the elastic body is of course deforming to accomodate the changes in the forces acting on it.

The equations governing the system are as follows. The force, $`f_i`$, carried by the $`i`$th spring is given by Hooke’s law,

$$f_i=k(u_iD),$$
(1)
where $`k`$ is the spring constant, $`D`$ is the displacement of the hard medium, and $`u_i`$ is the deformation of the elastic medium at site $`i`$. All unbroken springs have $`k=1`$ while a broken spring has $`k=0`$. The quantity $`(u_iD)`$ is, therefore, the length spring $`i`$ was stretched. In addition, a force applied at a point on an elastic surface will deform this surface over a region whose extent depends on its elastic properties. This is described by the coupled system of equations,

$$u_i=\underset{j}{}G_{i,j}f_j,$$
(2)
where the elastic Green function, $`G_{i,j}`$ is given by

$$G_{i,j}=\frac{1s^2}{\pi Ea^2}_{a/2}^{+a/2}\frac{dx^{}dy^{}}{|(xx^{},yy^{})|}.$$
(3)
In this equation, $`s`$ is the Poisson ratio, $`e`$ the elastic constant, and $`|\stackrel{}{i}\stackrel{}{j}|`$ the distance between sites $`i`$ and $`j`$. The indices $`i`$ and $`j`$ run over all $`L^2`$ sites. Strictly speaking, this Green function applies for a medium occupying the infinite half space. However, with a judicious choice of elastic constants, we may use it for a finite medium if its range is small compared to $`L`$, the size of the system.

By combining equations (1) and (2), we obtain

$$(𝐈+\mathrm{𝐊𝐆})\stackrel{}{f}=𝐊\stackrel{}{D},$$
(4)
where we are using matrix-vector notation. $`𝐈`$ is the $`L^2\times L^2`$ identity matrix, and $`𝐆`$ is the Green function represented as an $`L^2\times L^2`$ dense matrix. The constant vector $`\stackrel{}{D}`$ is $`L^2`$ dimensional. The diagonal matrix $`𝐊`$ is also $`L^2\times L^2`$. Its matrix elements are either 1, for unbroken springs, or 0 for broken ones. Of course $`𝐊`$ and $`𝐆`$ do not commute (except initially when there are no broken springs).

Once equation (4) is solved for the force, $`\stackrel{}{f}`$, equation (2) easily yields the deformations of the elastic surface.

### B Numerical Method: Fourier Acceleration

Equation (4) is of the familiar form $`𝐀\stackrel{}{x}=\stackrel{}{b}`$. Since the Green function connects all nodes to all other nodes, the $`L^2\times L^2`$ matrix $`𝐀`$ is dense which puts severe limits on the size of the system that may be studied. There are direct, time consuming methods to deal with such matrices, see Ref. . However, as we shall see, this problem may be circumvented and much more efficient methods may be employed such as the Conjugate Gradient algorithm (CG) .

The simulation proceeds as follows: We start with all springs present, each with its randomly drawn breakdown threshold. The two media are then pulled apart, the forces calculated using CG, and the spring which is the nearest to its threshold is broken, i.e. the matrix element corresponding to it in the matrix $`𝐊`$ is zeroed. Then the new forces are calculated, a new spring broken and so on till all springs have been broken and the media separated.

However, there are two problems that render the simulation of large systems extremely difficult. The first is that since $`𝐆`$ is $`L^2\times L^2`$ dense matrix, the number of operations per CG iteration scales like $`L^4`$. Even more serious is the fact that as the system evolves and springs are broken, the matrix $`(𝐈+k𝐆)`$ becomes very ill-conditioned.

To overcome the problematic $`L^4`$ scaling of the algorithm we note that the Green function is diagonal in Fourier space. Consequently, doing matrix-vector multiplications using FFTs the scaling is much improved and goes like $`L^2\mathrm{ln}(L)`$. Symbolically, this can be expressed as follow:

$$(𝐈+\mathrm{𝐊𝐅}^\mathrm{𝟏}\mathrm{𝐅𝐆})𝐅^\mathrm{𝟏}𝐅\stackrel{}{f}=𝐊\stackrel{}{D},$$
(5)
where $`𝐅`$ is the FFT operator and $`𝐅^\mathrm{𝟏}`$ its inverse ($`𝐅^\mathrm{𝟏}𝐅=1`$). Since $`𝐈`$ and $`𝐊`$ are diagonal, operations involving them are performed in real space. With this formulation, the number of operations/iteration in the CG algorithm now scales like $`L^2\mathrm{ln}(L)`$.

To overcome the runaway behavior due to the ill-conditioning we need to precondition the matrix . This means that instead of solving equation (5), we solve the equivalent problem

$$𝐐(𝐈+\mathrm{𝐊𝐅}^\mathrm{𝟏}\mathrm{𝐅𝐆})𝐅^\mathrm{𝟏}𝐅\stackrel{}{f}=\mathrm{𝐐𝐊}\stackrel{}{D},$$
(6)
where we simply multiplied both sides by the arbitrary, positive definite preconditioning matrix $`𝐐`$. Clearly, the ideal choice is $`𝐐_\mathrm{𝟎}=(𝐈+\mathrm{𝐊𝐆})^1`$ which would always solve the problem in one iteration. Since this is not possible in general, we look for a form for $`𝐐`$ which satisfies the following two conditions: (1) As close as possible to $`𝐐_\mathrm{𝟎}`$, and (2) fast to calculate. The choice of a good $`𝐐`$ is further complicated by the fact that as the system evolves and springs are broken, corresponding matrix elements of $`𝐊`$ are set to zero. So, the matrix $`(𝐈+\mathrm{𝐊𝐆})`$ evolves from the initial form $`(𝐈+𝐆)`$ to the final one $`𝐈`$. We were not able to find a fixed $`𝐐`$ that worked throughout the breaking process.

We therefore chose the form

$$𝐐=𝐈(\mathrm{𝐊𝐆})+(\mathrm{𝐊𝐆})(\mathrm{𝐊𝐆})(\mathrm{𝐊𝐆})(\mathrm{𝐊𝐆})(\mathrm{𝐊𝐆})+\mathrm{}$$
(7)
which is nothing but the Taylor series expansion of $`𝐐_\mathrm{𝟎}=(𝐈+\mathrm{𝐊𝐆})^1`$. For best performance, the number of terms kept in the expansion is left as a parameter since it depends on the physical parameters of the system. It is important to emphasize the following points. (a) As springs are broken, the preconditioning matrix evolves with the ill-conditioned matrix and, therefore, remains a good approximation of its inverse throughout the breaking process. (b) All matrix multiplications involving $`𝐆`$ are done using FFTs. (c) The calculation of $`𝐐`$ can be easily organized so that it scales like $`nL^2\mathrm{ln}(L)`$ where $`n`$ is the number of term kept in the Taylor expansion, equation (7).

We therefore have a stable accelerated algorithm which scales essentially as the volume of the system. For example, for a $`128\times 128`$ system, and taking $`n=5`$, the CG algorithm always converged in four or five iterations with the prescribed precision of $`10^{12}`$.

## III Results

We now present the results of our numerical simulations. We show in Fig. 1 a representation of the failure process. Each elementary square represents a spring (a bond), and the gray scale indicates when a particular spring failed: The darker the color, the earlier the failure. In this particular case, the elastic constant $`e`$ was set to 10. There are no apparent spatial correlations between the failing bonds in this figure. However, we show in Fig. 2 the distance between successively failing bonds for the same disorder realization of Fig. 1. We see clearly in this figure that early in the process there is no localization effect: Bonds tend to break far apart, the location being determined by the strength of bonds, i.e. early failure is disorder driven. However, halfway into the breakdown process, localization clearly develops. In Fig. 3, we show the corresponding plot for $`e=100`$. In this case localization never develops for this size system and destribution of thresholds.

On the other hand, if the threshold distribution is much narrower than $`[0,1]`$ used above, localization can develop early. For example, we show in Fig. 4 the fracture graymap (like Fig. 1) for a uniform threshold distribution in the interval $`[9.5,10.5]`$. We clearly see the fracture starting towards the center of the figure and spreading out in a spiral till finally the symmetry is broken and the system ruptures along one of the lattice directions.

Fig. 5 shows the force-displacement curve for a system with elastic constant $`e=10`$. Whether we control the applied force, $`F`$, or the displacement, $`D`$, the system will eventually suffer catastrophic collapse. However, this is not so when $`e=100`$ as shown in Fig. 6. In this case, only controlling force will lead to catastrophic failure. In the limit when $`e\mathrm{}`$, the model becomes the global load-sharing fiber bundle model , where $`F=(1D)D`$. In this limit there are no spatial correlations and the force instability is due to the the decreasing total elastic constant of the system making the force on each surviving bond increase faster than the typical spread of threshold values. No such effect exists when controlling displacement $`D`$. However, when the elastic constant, $`e`$, is small, spatial correlations in the form of localization do develop, and these are responsible for the diplacement instability which is seen in Fig. 5. In other words, the localization clearly visible in Fig. 2 starts to develop when the system is near the peak of its force-displacement curve, and dominates when the system is on the negative slope branch of that curve.

### A Burst Distribution

We now turn to the study of the burst distribution. We define the size of a burst, $`\mathrm{\Delta }`$, in our model as the number of bonds that fail simultaneously while the force $`F`$ is held constant. In the global load-sharing fiber bundle model it has been shown that the burst distribution is given by

$$N(\mathrm{\Delta },D)=\frac{1}{\mathrm{\Delta }^\tau }n\left(\mathrm{\Delta }^\sigma (xx_c)\right)$$
(8)
where $`x_c`$ is the damage, i.e. the density of broken bonds, at which the model fails catastrophically, $`n`$ is a crossover function which approaches a constant when the argument approaches zero, and which falls off as $`\mathrm{exp}(y^2)`$ as the argument $`y`$ is large. Furthermore,

$$\tau =\frac{5}{2}\mathrm{and}\sigma =\frac{1}{2}$$
(9)
independent of the threshold distribution.

We show in Figs. 7 and 8 the burst distribution for $`e=10`$ and 100. In both cases we find that the burst distribution follows a power law with an exponent $`\tau =2.6\pm 0.1`$. We may argue that the exponent is the same as the one found in the global loading fiber bundle, Eq. (9), $`\tau =5/2`$ in the following way. The characteristics, $`F=F(x)`$ must have a quadratic maximum somewhere. For $`e=100`$, such a maximum exists in the middle of failure process as seen in Fig. 8, whereas for $`e=10`$, the system only approaches such a maximum near total failure, see Fig. 7. Assuming that the fluctuations about the average characteristics are brownian — which can be shown analytically in the limit $`e\mathrm{}`$  — near the maximum the probability to find a burst of size $`\mathrm{\Delta }`$ is proportional to $`\mathrm{\Delta }^{3/2}\mathrm{exp}[\mathrm{\Delta }(xx_c)^2]`$. This result comes from a mapping onto the Gambler’s ruin problem . Furthermore, in order to guarantee that the burst is not a burst within an even larger burst, the starting point of the burst must be the highest point on the characteristics that has occured so far in the failure process. This condition may also be mapped onto the Gambler’s ruin problem, and leads to an extra factor $`(xx_c)`$ in the probability for a burst to occur. The probability to find a burst of size $`\mathrm{\Delta }`$ throughout the failure process is then the integral over $`x`$ as $`x`$ approaches $`x_c`$, $`^{x_c}𝑑x(xx_c)\mathrm{\Delta }^{3/2}\mathrm{exp}[\mathrm{\Delta }(xx_c)^2]`$ which is proportional to $`\mathrm{\Delta }^{5/2}`$.

As can be seen in Figs. 7 and 8, the numerical data are consistent with the expected value $`\tau =5/2`$.

### B Strength Scaling

The load-displacement curves for different system sizes $`L`$ coincide when we use the reduced variables $`F/L^a`$ and $`D/L^b`$ where $`a=2.0`$ and $`b=0.0`$, as seen in Fig. 9. We expect the exponent $`a=2`$ since $`F/L^2`$ is the normal stress on the surface. In the case of an infinitely stiff system, we expect $`b=0`$. The elastic system studied here behaves in the same way as long as the elastic constant, $`e`$, is also scaled with $`L`$. For example for $`L=128`$ we took $`e=10`$, for $`L=64`$ $`e`$ should take half that value in order to reproduce the physics. This is easy to understand considering the dependence of the Green fucntion, Eq. 3, on the elastic constant.

### C Spatial Damage Distribution at Failure

As the failure process proceeds, there is an increasing competition between local failure due to stress enhancement and local failure due to local weakness of material. As we saw above, when we control the displacement, $`D`$, and $`e`$ is sufficiently small (for example $`e=10`$), catastrophic failure eventually occurs due to localization. The onset of this localization, i.e. the catastrophic regime, occurs when the two mechanisms are equally important. One may suspect that criticality due to self organization  occurs at this point. In order to test whether this is the case, we have measured the size distribution of broken bond clusters at the point when $`D`$ reaches its maximum point on the $`FD`$ characteristics, i.e. the onset of localization and catastrophic failure. The analysis was performed using a Hoshen-Kopelman algorithm . We show the result in Fig. 10, for 56 disorder realizations, $`L=128`$ and $`e=10`$. The result is consistent with a power law distribution with exponent $`1.6`$, and consequently with self organization.

## IV Conclusion

We have studied numerically the failure of the glued interface between an elastic and an infinitely stiff blocks of material. To this end we have developed a new, stable and accelerated algorithm which scales $`L^2\mathrm{ln}(L)`$ which enabled us to study much bigger systems than previously possible.

Our main physical results are: (1) the distribution of simultaneously failing glue (bursts) is a power law with exponent $`2.6\pm 0.1`$, which is consistent with the burst distribution found in the global load sharing fiber bundle problem. (2) The point of catastrophic failure scales as $`L^2`$ in force and $`L^0`$ in displacement. (3) The area distribution of failed regions (clusters) at the onset of catastrophic failure when displacement is the control parameter is consistent with a power law with an exponent equal to $`1.6`$. This hints at self organization.

In addition, we saw that for large $`e`$, e.g. $`e=100`$, the system does not suffer catastrophic failure, and there is no localization. On the other haand, smaller values of $`e`$, e.g. $`e=10`$, resulted in catastrophic failure due to localization. By doing the simulations for various values of $`e`$ we estimate that failure due to localization starts to occur for $`e3540`$. As we will see below, these values of $`e`$ obtained for $`128\times 128`$ systems should be scaled appropriately when the size of the system is changed.

The disorder in our system was uncorrelated. As mentioned above, it is realistic to introduce correlations as exist for example in fracture surfaces. This can be done by generating spring breaking thresholds that have the desired correlations. Furthermore, we have used a flat distribution for the disorder. One can use other distributions, e.g. $`r^\alpha `$, where $`r`$ is a uniformly distributed random number and $`\alpha `$ an exponent that can be negative. It is known from random fuse models of fracture that the breakdown process depends on the value of $`\alpha `$. It is not clear how these issues will modify our current results. This work is in progress.

Another work in progress is to study the propagating fracture front when the above glued media are ripped apart by pulling only on one side of the $`L\times L`$ square system. As before the breaking thresholds can be correlated or uncorrelated. These results will also be compared with the results of experiment currently underway.

Finally, we have chosen to introduce disorder into the breaking thresholds of the springs. However, we can just as easily introduce it into the spring constants themselves, again with or without correlations.

###### Acknowledgements.

We thank F. A. Oliveira and H. Nazareno and the ICCMP of the Universidade de Brasília for warm hospitality and support during the initial phases of this work. This work was partially funded by the CNRS PICS contract $`\mathrm{\#}753`$ and the Norwegian research council, NFR. We also thank NORDITA for its hospitality and further support. A. H. thanks the Niels Bohr Institute for their support and hospitality during his sabbatical year.


text:
# Separation of time and length scales in spin-glasses: temperature as a microscope

## 1 Facts and puzzles

Although spin-glasses are totally useless pieces of material, they constitute an exceptionally convenient laboratory frame for theoretical and experimental investigations . Theoretical concepts and experimental protocols relevant for more ‘useful’ glassy materials (polymer and molecular glasses, foams and pastes, etc.) have been elaborated and tested on spin-glasses . There are at least two reasons for this: (a) the theoretical models are conceptually simpler (although still highly non trivial) and (b) the use of very sensitive magnetic detectors allows one to probe in details the a.c and d.c spin dynamics of these systems down to very small external fields. The corresponding mechanical measurements in other glassy systems are much more difficult to control, although some recent progress have been made , in particular concerning the measurement of age-dependent structure factors .

The aging dynamics of spin glasses has therefore been studied in glory details recently, and has revealed an extremely rich phenomenology . The most striking aspect is the role of small temperature changes, that we summarize as follows: (a) Superactivated behaviour: time scales grow faster than what would be expected from simple thermal activation when the temperature is decreased; (b) Rejuvenation and memory: after a small negative temperature jump (within the glass phase) the system behaves as if it had been quenched from above the glass transition temperature $`T_c`$ (rejuvenation). However, a perfect memory of the time spent at the initial temperature is somehow kept, as clearly demonstrated by the now well known ‘dip’ imprinting experiments (see Fig. 1) ; and (c) Weak cooling rate dependence: the a.c. susceptibility hardly depends on the thermal history, or actually only on the cooling rate over the very last Kelvins before reaching the final temperature. The cooling rate at higher temperature, in particular when crossing $`T_c`$, is irrelevant, in strong contrast with Random Field like systems . This absence of cooling rate dependence is in fact another manifestation of rejuvenation.

It is now well established, both experimentally  and numerically  that a certain ‘coherence’ length is growing in an aging spin-glass. This was first predicted in the context of the ‘droplet model’ , but is presumably of much more general validity: larger length scales take a longer time to evolve. This is expected to be true even if the basic tenets of the droplet model turn out to be incorrect. The coherence length is found to grow as a power of time: $`\mathrm{}t^{1/z}`$, with an apparent exponent $`1/z`$ linear in temperature . This suggests an activated behaviour over barriers that grow as the logarithm of the length $`\mathrm{}`$. Indeed, writing $`t(\mathrm{})\mathrm{exp}[\mathrm{\Delta }(\mathrm{})/k_BT]`$ with $`\mathrm{\Delta }(\mathrm{})=\mathrm{\Delta }_0\mathrm{log}\mathrm{}`$ leads to $`1/z=k_BT/\mathrm{\Delta }_0`$. This is confusing because (a) as mentioned above, experiments suggest superactivated behaviour; (b) barriers should grow with $`\mathrm{}`$ faster than excitations energies, which are thought to grow as $`\mathrm{}^\theta `$ with $`\theta =0.2`$ in the droplet model. (Note however that recent numerical simulations suggest that non compact excitations indeed correspond to $`\theta =0`$ ); and (c) the exponent $`z`$ determined numerically or experimentally can be written as $`z(T)=z_cT_c/T`$, where $`z_c6`$ is the critical exponent that governs the dynamics of the spin-glass at the critical point. This coincidence suggests that the system is somehow affected by critical fluctuations (as also suggested, on the basis of different arguments, in ).

The rejuvenation effect (and absence of cooling rate dependence) brings still more confusion. Activated barrier crossing is obviously easier at higher temperature: so why waiting longer around $`T_c`$ does not help to equilibrate the system, as happens in e.g. random field like systems? (see the e.g. the discussion in ). This should be even more crucial if the dynamics is superactivated. A way out of this contradiction is to invoke chaos in temperature : if new patterns need to equilibrate when the temperature is changed, it is clear that the time spent at a higher temperature does not help much to equilibrate the system at the final temperature. However (a) no sign of chaos has been found in most recent static numerical studies of the 3d Edwards-Anderson model  (at variance with earlier studies ), nor in the theoretical analysis of the sk model ; (b) no rejuvenation in the dynamics upon temperature changes has been found either in numerical simulations ; and (c) the coherence length that has grown at one given temperature seems to carry on growing (although at a different rate) at another temperature, in contradiction with the chaos idea. This continuity has been established both numerically  and experimentally . ‘Chaos’ should furthermore be compatible with memory: growing new patterns at a given temperature should not erase the patterns grown at higher temperatures if one wants to account for the ‘multi-dip’ experiment shown in Fig. 1. A possible scenario for this was suggested in .

In this paper, we wish to develop a consistent qualitative picture for the dynamics in spin-glasses that allows one to resolve the above apparent contradictions. This picture, as the original droplet model , heavily relies on three basic ideas: (i) times grow as the exponential of the energy barriers; (ii) the energy barriers grow as a power of the length scales involved in the dynamics; and (iii) the energy barriers vanish at the critical temperature $`T_c`$. These ingredients are enough to understand that changing the temperature for a given observation time corresponds to changing the length scale at which the system is probed. Therefore, the objects contributing to the aging dynamics are different at different temperatures. This gives a precise content to the ‘hierarchical’ picture advocated in many papers . Although the present picture is similar in spirit to the droplet model, we will assume (but this is not crucial) that the low energy metastable states exist on arbitrary length scales (i.e. $`\theta =0`$). We will also introduce the idea of a fixed landscape rejuvenation, which concerns short length scales and is distinct from the large scale chaos that appears in the context of the droplet model. We reanalyze previous experimental data within this framework, which strongly suggest that energy barriers vanish at $`T_c`$, and allow us to extract estimates of the exponent $`\psi `$ that relates energy barriers and length scales.

## 2 Basic ingredients

Let us consider a large scale low lying excitation in a spin glass. This excitation is made of a large connected cluster of spins that is flipped with respect to the ground state. The surface of this cluster can be thought of as a ‘domain wall’ which happens to occupy a very favorable position since the overall energy of this surface is very small. Therefore, this wall is ‘pinned’ by the disorder and tends to adopt some special conformation. If $`\theta =0`$, this wall has no overall tendency to disappear with time, and will be present in the system even after very long times. If $`\theta >0`$, conversely, these large scale walls tend to fraction in smaller and smaller bubbles before disappearing in the equilibrium state. However, in some exceptional circumstances, the energy of these droplets is smaller than $`k_BT`$, and these walls survive in equilibrium. In the droplet picture, this occurs with probability $`k_BT/\mathrm{}^\theta `$ .

Now, there are many conformations of the domain wall which have approximately the same energy. One can flip clusters of spins that touch the domain wall at a small cost, corresponding to a local modification (‘blister’) of the conformation of the wall. These excitations can occur on all length scales: one can create blisters within blisters, etc. The situation here is not specific to spin-glasses but is also true for a domain wall in a disordered ferromagnet that can adopt many different metastable configuration. An important difference is that in a disordered ferromagnet, these domain walls have a positive energy and tend to disappear with time: this is the coarsening phenomenon. The evolution of the system is in this case a slow but irreversible march towards order . The evolution in phase space is biased by the fact that these walls cost energy.

It is useful to decompose the conformations of these pinned ‘domain walls’ on different length scales. Identical conformations on length scale $`\mathrm{}_n`$ may differ by the presence or absence of blisters of smaller sizes $`\mathrm{}_{n1},\mathrm{}_{n2},\mathrm{}`$, where $`\mathrm{}_n=b^n\mathrm{}_0`$, and $`b`$ an arbitrary factor, say $`b=2`$. Here again, we follow ideas developed in the context of the droplet model for spin-glasses  or pinned domain walls . In the droplet model, the time needed to evolve the conformation on scale $`\mathrm{}_n`$ is taken to be:

$$t_n=t(\mathrm{}_n)\tau _0\mathrm{exp}\left(\frac{\mathrm{{\rm Y}}\mathrm{}_n^\psi }{k_BT}\right)$$
(1)
where $`\tau _0`$ is a microscopic time, $`\mathrm{{\rm Y}}`$ a typical energy setting the scale of energy barriers between conformations, and $`\psi `$ the so-called barrier exponent. This form has an immediate consequence: the time needed to evolve the system on scale $`\mathrm{}_n`$ is extremely long compared to the time needed to evolve the system on scales $`\mathrm{}_{n1},\mathrm{}_{n2},\mathrm{}`$. This means, as emphasized in , that on a time scale $`t_n`$, all excitations on scales $`\mathrm{}_n^{}`$ with $`n^{}>n`$ are essentially frozen, whereas all excitations on scales $`\mathrm{}_n^{}`$ with $`n^{}<n`$ are essentially equilibrated. Note also that short length scales are ‘slaved’ to large length scales: when a large length scale flips over, all the smaller length scales have to re-equilibrate in a new environment. In this sense, the dynamics is hierarchical (see a related discussion in .

We will choose to write Eq.(1) in a generalized form, more appropriate to describe the vicinity of the spin-glass transition111In the following, all length scales are expressed in lattice size units.:

$$t_n=t(\mathrm{}_n,T)\tau _0\mathrm{}_n^{z_c}\mathrm{exp}\left(\frac{\mathrm{{\rm Y}}(T)\mathrm{}_n^\psi }{k_BT}\right),$$
(2)
with $`\mathrm{{\rm Y}}(T)=\mathrm{{\rm Y}}_0[T_cT/T_c]^{\psi \nu }`$, $`\nu `$ being the critical exponent governing the divergence of the equilibrium correlation length $`\xi (T)`$ at $`T_c`$: $`\xi (T)=|T_cT|^\nu `$. Therefore, the term in the exponential can be rewritten as: $`(\mathrm{{\rm Y}}_0/k_BT)[\mathrm{}_n/\xi (T)]^\psi `$. Since $`\mathrm{{\rm Y}}_0`$ is expected to be of the order of $`k_BT_c`$, one sees that as long as $`\mathrm{}_n`$ is smaller than $`\xi (T)`$, barriers are small compared to $`k_BT`$ and the exponential term in Eq. (2) can be set to $`1`$. This leads to the usual critical dynamics (non activated) relation:

$$t_n\tau _0\mathrm{}_n^{z_c}.$$
(3)
Conversely, for length scales larger than $`\xi `$, the exponential becomes the dominant factor, and one recovers (to logarithmic accuracy) the Fisher-Huse relation (1) with the microscopic time $`\tau _0`$ replaced by the typical critical time scale $`\tau _0\xi ^{z_c}`$. As we shall find below, the experiments are typically in a crossover regime where barriers are larger than, but comparable to $`k_BT`$. Therefore, we keep the full form of Eq. (2), which correctly interpolates between the two regimes, to describe the intermediate regime.

Several interesting consequences of Eq. (2) are worth discussing.

* Fig. 2 shows $`\mathrm{log}_{10}t_n`$ versus $`\mathrm{}_n`$ for a choice of parameters suggested by the experiments on AgMn: $`\psi =1.5`$, $`\mathrm{{\rm Y}}_0/k_BT_c=2.`$, $`z_c=5`$ and $`\nu =1.3`$, and for different values of $`T/T_c`$. The thick horizontal line corresponds to $`t/\tau _0=10^{15}`$, corresponding to an experimental time scale of $`1000`$ seconds. One sees that the associated length scales are very modest, in the range 10 to 100, and change appreciably when the temperature is changed only slightly. The horizontal dotted line corresponds to numerically accessible time scales. We see that in that case, length scales are extremely small and do not separate at all with temperature.
* Let us fix a certain length scale $`\mathrm{}`$ and change the temperature from $`T_1`$ to $`T_2=T_1\mathrm{\Delta }T`$. The associated time scale $`t_{n1}`$ at temperature $`T_1`$ is changed into $`t_{n2}`$ given by:

$$\frac{t_{n2}}{\tau _{\mathrm{}}}=\left(\frac{t_{n1}}{\tau _{\mathrm{}}}\right)^\beta \beta =\frac{T_1}{T_2}\left(\frac{T_cT_2}{T_cT_1}\right)^{\psi \nu }$$
(4)
with $`\tau _{\mathrm{}}=\tau _0\mathrm{}^{z_c}`$. This leads to superactivation effects. For example, with $`T_1=0.9T_c`$ and $`T_2=0.8T_c`$, $`\psi \nu =2`$, one finds $`\beta 4`$ ! As soon as $`t_{n1}\tau _{\mathrm{}}`$, the value of $`t_{n2}`$ is astronomically large. This means that the separation of time scales when the temperature is lowered is extremely fast, and to a good approximation, those length scales that are aging at temperature $`T_1`$ become completely frozen at a slightly lower temperature. In this sense, temperature acts as a microscope. If $`\mathrm{{\rm Y}}`$ was temperature independent, one would find $`\beta =T_1/T_2=1.125`$ for the above choice of parameters. The separation of time scales would then only be mild.
* If one fits Eq. (2) with a temperature dependent power law over a restricted range of times, one finds an effective exponent:

$$z_{eff}=\frac{d\mathrm{log}t_n}{d\mathrm{log}\mathrm{}_n}=z_c+\psi \frac{\mathrm{{\rm Y}}(T)\mathrm{}_n^\psi }{k_BT},$$
(5)
that grows when the temperature is lowered. We will come back to this point later.

## 3 Old experiments revisited

### 3.1 Small temperature jumps

We now revisit two sets of experiments and show that the results are indeed compatible with Eq. (2). The first set of experiments concerns trm relaxation. One first measures an array of curves corresponding to the standard protocol: cool the system from above $`T_c`$ to $`T_1`$, leave a (small) magnetic field on for a certain waiting time $`t_{w1}`$, and observe the relaxation at time $`t_{w1}+t`$. One finds curves that scale approximately as $`t/t_{w1}`$ once the stationary (fast) initial part has been subtracted. Then, a second protocol is followed: cool now to $`T_2=T_1\mathrm{\Delta }T`$, and wait for a time $`t_{w2}`$. Then heat up the system to $`T_1`$ and simultaneously cut the field. What effective waiting time $`t_{w1}^{eff}`$ should one choose to match this second type of trm with a standard, isothermal one ? The experimental results show that as long as $`\mathrm{\Delta }T`$ is small enough, it is always possible to find such a $`t_{w1}^{eff}`$ so that the curves match perfectly. This means that the objects involved in the dynamics are exactly the same for the two temperatures. For larger $`\mathrm{\Delta }T`$’s, the curves are distorted and such a perfect matching is impossible. Typically, experiments were performed with very small $`\mathrm{\Delta }T=0.02`$ K ($`\mathrm{\Delta }T/T_g0.2\%`$). Obviously, since the time spent at $`T_2`$ corresponds to a smaller time spend at $`T_1`$, one finds $`t_{w1}^{eff}<t_{w2}`$, see Table. It was furthermore shown in  that the correspondence between $`t_{w1}^{eff}`$ and $`t_{w2}`$ could not be understood in terms of simple thermal activation. More precisely, one finds that $`\mathrm{log}t_{w1}^{eff}/\tau _0<(T_2/T_1)\mathrm{log}t_{w2}/\tau _0`$, corresponding to superactivation (unless $`\tau _0`$ is chosen to be unphysically small). This was interpreted in  as indicating a divergence of the corresponding barrier at smaller temperatures. Note however that a barrier involving a finite number of spins cannot diverge at any temperature. Another interpretation is that barriers actually vanish at $`T_c`$, as in the droplet model. This is reasonable since the ‘domain walls’ (whatever their precise nature) become more and more loosely defined and can no longer be pinned at $`T_c`$ (this is also true in a disordered ferromagnet – but not in a random field system: see ).

We have therefore reanalyzed the very clean experimental data on Ag Mn of ref.  by postulating Eq. (2). We fix $`z_c=5`$ and $`\nu =1.3`$ to reasonable values and determine $`\psi `$ and $`\mathrm{{\rm Y}}_0/k_BT_c`$ such that the length scale $`\mathrm{}_1`$ corresponding to the time $`t_{w1}^{eff}`$ and the length scale $`\mathrm{}_2`$ corresponding to the associated time $`t_{w2}`$ are as close as possible: since the shape of the trm’s are made to coincide, the corresponding length scales should also coincide. Therefore, we choose $`\psi `$ and $`\mathrm{{\rm Y}}_0/k_BT_c`$ such that the mean squared relative difference:

$$^2=\underset{i=1}{\overset{N}{}}\left(\frac{\mathrm{}_1^i\mathrm{}_2^i}{\mathrm{}_2^i}\right)^2,$$
(6)
summed over all experiments reported in Table 1, is as small as possible. The best values are found to be in a ‘crescent’ in the $`\psi ,\mathrm{{\rm Y}}`$ plane with for example $`\psi =1.5`$ and $`\mathrm{{\rm Y}}_0/k_BT_c=2.`$, corresponding to a rather small root mean square relative error of $`=0.48\%`$, or $`\psi =2`$ and $`\mathrm{{\rm Y}}_0/k_BT_c=0.3`$, corresponding to an error of $`=0.40\%`$ (the experimental error on $`t_{w1}^{eff}`$ corresponds to a relative error on $`\mathrm{}`$ on the order of $`0.5\%`$). The typical length scales obtained are shown in Fig. 2. The value $`\psi =2`$ corresponds with the upper bound proposed by Fisher and Huse for an Ising spin glass. A similar range of values ($`\psi =1.3`$ for $`\mathrm{{\rm Y}}_0/k_BT_c=1`$) is found for the insulating CrIn compound , using the same procedure but on less precise data (larger $`\mathrm{\Delta }T`$). The point here is not to claim a very good precision on the value of $`\psi `$, but rather to show that the results are compatible with the idea that the barriers continuously vanish at $`T_c`$. Note that more recent experiments on the role of small temperature jump on Ising spin-glasses confirm the present analysis, although the separation of time scales is much ‘milder’: the dynamics for $`T>0.6T_c`$ involve both activated events over temperature dependent barriers and critical dynamics that increases the effective value of the ‘trial time’ $`\tau _{\mathrm{}}=\tau _0\mathrm{}^{z_c}`$. 222The fact that a rather large value of this trial time was needed to account for the experiments was also noticed in . However, the value of $`\psi `$ for this Ising-like sample is significantly smaller ($`\psi 0.30.5`$) than is the less anisotropic sample reported above, a somewhat counter-intuitive result.

### 3.2 A time dependent length probed by magnetic field

We now turn to another set of more recent experiments , which exploits the fact that a small magnetic field acts as to reduce the energy barriers, suggested in . If the number of spins involved in a re-conformation is $`N`$, one can expect that a field $`H`$ will perturb the barriers by an amount proportional to $`N\chi H^2`$, where $`\chi `$ is the magnetic susceptibility. The typical relaxation time of the trm aged for a time $`t_w`$ is therefore multiplied by $`\mathrm{exp}(\alpha N(t_w,T)\chi H^2/k_BT)`$, where $`\alpha `$ is a numerical factor. By measuring this reduction factor for different waiting times $`t_w`$, one can estimate the typical number of spins involved in the dynamics as a function of the waiting time. Finally, by writing $`N(t_w,T)\mathrm{}^3`$, one has access to a time dependent coherence length $`\mathrm{}(t_w,T)`$. As mentioned in the introduction, this dependence can be fitted, for three different types of spin glasses, by a power law with a temperature dependent exponent $`z(T)`$. Note that the estimated number of spins, in the range $`10^510^6`$, is compatible with the length scales reported in Fig 2.

Here, we want to reanalyze the data of  in the light of Eq. (2), and show that again, these experimental results suggest that barriers vanish at $`T_c`$. In order to do so, we have plotted the quantity:

$$G(t_w,T)=\left(\frac{\mathrm{log}t_w/\tau _0\frac{z_c}{3}\mathrm{log}N(t_w,T)}{\frac{T_c}{T}N(t_w,T)^{\psi /3}}\right)^{1/\psi \nu }$$
(7)
as a function of $`T/T_c`$, for different spin-glasses. If Eq. (2) is correct, one should observe $`G(t_w,T)=G_0(1T/T_c)`$, where $`G_0`$ is a numerical factor. The vanishing of $`G(t_w,T)`$ is direct manifestation of the vanishing of the barriers. The results are shown in Fig 3. We have kept the values $`z=5`$, $`\psi =1.5`$ suggested above. A linear fit through the points is very reasonable; the most interesting point is that this linear fit is found to be $`G_{fit}(t_w,T)=0.58(1.025T/T_c)`$, very close to what is expected from Eq. (2). Note that the extrapolated value of $`T/T_c`$ is even closer to one if $`z_c`$ is chosen to be equal to $`6`$. Therefore, the power-law dependence of $`\mathrm{}`$ reported in  (and also in the numerical work of ) might actually be an effective power law, as suggested by Eq. (5) above, which naturally matches the critical dynamics exponent when $`TT_c`$.

In summary, we have shown in this section that two completely independent sets of experiments can be interpreted consistently within the framework of section 2: time scales and length scales are related by Eq. (2). The most important aspect is the fact that energy barriers vanish at $`T_c`$. This means that (a) dynamics is superactivated and time scales separate extremely fast in spin-glasses, and (b) the critical point does significantly affect the dynamics by slowing down the ‘microscopic’ frequency.

## 4 Fixed landscape rejuvenation and Memory

### 4.1 Qualitative ideas

Let us now come back in more details to the temperature cycling experiments of  and show how these can be qualitatively interpreted within the general picture of section 2. The crucial observation is that for a given time scale $`t_{w1}`$ and temperature $`T_1`$, the aging dynamics is dominated by a characteristic length $`\mathrm{}_1`$ such that $`t(\mathrm{}_1,T_1)=t_{w1}`$. Larger length scales are essentially frozen and do not evolve on the time scale of the experiment, while shorter length scales are fully equilibrated (for a given larger scale conformation) and only contribute to the stationary part of the response function. Note that this picture is obviously a caricature because the energy barriers in disordered systems are expected to fluctuate in space. Therefore, some larger length scales might locally see an exceptionally low barrier, and vice versa.

Now, let the temperature change from $`T_1`$ to $`T_2=T_1\mathrm{\Delta }T`$. Because of the separation of time scales, a relatively small $`\mathrm{\Delta }T`$ is sufficient to freeze completely the dynamics on scale $`\mathrm{}_1`$ (which will therefore retain the memory of the stay at $`T_1`$) and to slow down the initially fast dynamics on shorter length scales, selecting a particular one $`\mathrm{}_2`$ to be in the experimental time window. Since at $`T_1`$ this length scale is equilibrated, the different conformations appear in the course of time with their Boltzmann weights. As the temperature is changed, these Boltzmann weights are modified and the system has to evolve towards a new state. This is true even if the (free)-energy landscape does not significantly evolve between the two temperatures: this is what we call fixed landscape rejuvenation. Take for example a simple two-level (for example two conformations of a domain wall differing by a ‘blister’ of scale $`\mathrm{}`$), with an energy difference $`E(\mathrm{})`$. The population difference between the two levels change significantly (say by at least 10%) if the change of temperature is such that:

$$\mathrm{\Delta }T\frac{\mathrm{tanh}(\frac{E(\mathrm{})}{2k_BT})}{T}=\frac{\mathrm{\Delta }TE(\mathrm{})}{2k_BT^2}\frac{1}{\mathrm{cosh}^2(\frac{E(\mathrm{})}{2k_BT})}0.1$$
(8)
The minimal value of $`\mathrm{\Delta }T`$ for such a rearrangement to occur is therefore $`\mathrm{\Delta }T^{}/T0.15`$ which is corresponds to $`E(\mathrm{})2.4k_BT`$. (For multi-level systems, $`\mathrm{\Delta }T^{}`$ can be much smaller than this, see below and ). It is furthermore easy to show that when a two-level system is driven out-of-equilibrium by a rapid change of temperature, an excess low-frequency dissipation follows. The out of phase a.c. susceptibility is indeed given by:

$$\chi ^{\prime \prime }(\omega ,t_w)\frac{\mathrm{exp}[t_w/t(\mathrm{},T_2)]}{\omega t(\mathrm{},T_2)}\omega t(\mathrm{},T_2)1$$
(9)
where $`t(\mathrm{})`$ is the relaxation time of the two-level system. Therefore, cooling a disordered system can induce a strong increase in the out of phase susceptibility (rejuvenation) if there are metastable states such that $`t(\mathrm{}_2,T_2)\omega ^1`$ and $`E(\mathrm{}_2)k_BT`$.

Lower frequencies therefore probe larger length scales $`\mathrm{}_2`$. In the droplet model, $`E(\mathrm{})`$ is typically of order $`\mathrm{{\rm Y}}\mathrm{}^\theta `$, and the probability to observe an ‘active’ droplet of energy of order $`k_BT`$ is $`(k_BT/\mathrm{{\rm Y}})\mathrm{}^\theta `$. If $`\theta `$ is positive, we expect that rejuvenation should asymptotically disappear as $`\omega 0`$ since the probability of observing a large droplets goes to zero. However, since $`\mathrm{}_2`$ depends logarithmically on frequency and $`\theta `$ is small ($`0.2`$), this will in practice never happen and things will look very much as if $`\theta =0`$. In other words, the fact that $`\mathrm{}`$ cannot much exceed $`100`$ in experiments and that $`\theta `$ is so small, means that $`\mathrm{}^\theta 0.4`$! The influence of large scale, low energy excitations (of order $`k_BT`$), is therefore dominant in real spin-glasses. This feature is actually the basic outcome of mean-field models.

As emphasized above, the strong separation of time scales enables one to observe simultaneously rejuvenation and memory if $`\mathrm{\Delta }T`$ is large enough: the length scales that one observes at $`T_1`$ are totally frozen at temperature $`T_2`$ and therefore resume aging, unaffected by the long stay at smaller temperatures. Of course, the states on scale $`\mathrm{}_2`$ are now out of equilibrium at $`T_1`$; however, the time needed for them to equilibrate is very short, since it is given by:

$$t(\mathrm{}_2,T_1)=\tau _\mathrm{}_2\left(\frac{t(\mathrm{}_2,T_2)}{\tau _\mathrm{}_2}\right)^{1/\beta }=\tau _0^{11/\beta }\mathrm{}_2^{z_c(11/\beta )}\omega ^{1/\beta }$$
(10)
Taking for example $`\beta =2`$, $`\mathrm{}_2=10`$, we find that $`\omega t(\mathrm{}_2,T_1)=10^3(\omega \tau _0)^{1/2}10^3`$ for $`\omega =1`$ Hz. Therefore, these length scales have equilibrated far before the first oscillation of the a.c. field has taken place.

It is interesting to remark that in the scenario, memory cannot be observed after a positive temperature cycle $`T_1T_1+\mathrm{\Delta }T`$, if $`\mathrm{\Delta }T`$ is large enough to induce rejuvenation. This is because upon heating, length scales larger than $`\mathrm{}_1`$ will now un-freeze and evolve. Since the shorter length scales are slaved to the larger length scales, the length scale $`\mathrm{}_1`$ will itself have to re-equilibrate completely at higher temperatures – this erases the memory of the stay at $`T_1`$. As will be discussed below, this is not necessarily true in the ‘chaos’ scenario. Note that the ‘slaving’ of small length scales to large length scales also explains qualitatively the ‘second noise spectrum’ experiments of Weissmann , which reflect the fact that when a large length scale ‘jumps’, all the smaller length scales have to re-equilibrate in a new environment, uncorrelated with the first .

Before describing a more realistic model, let us summarize the above discussion by a schematic figure in the plane $`\mathrm{},\mathrm{\Delta }T`$, showing how the different length scales evolve during the temperature cycle. This figure allows to account qualitatively for the temperature cycling experiments, in particular the fact that memory is not perfect if $`t_{w2}`$ is large enough or $`\mathrm{\Delta }T`$ small enough (see also ). If $`\mathrm{\Delta }T<\mathrm{\Delta }T^{}`$, no rejuvenation is expected and aging continues, although at a slower rate. Rejuvenation and memory are made possible by the strong length and time scale separation as a function of temperature for a fixed time scale. As shown in Fig. 2 this is not true in numerical simulations, where one cannot probe the multi-scale dynamics of the system.

Note finally that the above discussion is not restricted to spin-glasses, and can be applied to describe the aging dynamics of pinned domain walls in ferromagnets , ferroelectrics , and glassy polymers  where indeed similar effects have been observed experimentally. There is a significant difference, though, which is due to domain coarsening. The dynamics of the system progressively gets rid of the domain walls; this corresponds to a cumulative aging effect, which is cooling rate dependent. As time evolves, the fraction of spins that belong to domain walls and contribute to aging dynamics systematically decreases , as observed in these systems .

### 4.2 More sophisticated models

The simple two-level picture discussed above is obviously not sufficient to explain in details the experimental data. For example, the aging behaviour of $`\chi ^{\prime \prime }`$ is well fitted by a power-law: $`\chi _{AG}^{\prime \prime }(\omega ,t_w)(\omega t_w)^b`$, with $`b0.2`$. This is quite different from Eq. (9), although the latter correctly predicts a strong increase of $`\chi ^{\prime \prime }`$ at low frequencies. On a given length scale $`\mathrm{}`$, there are actually many metastable conformations of the domain walls, between which the system jumps. One expects that the energy barriers between these metastable states are distributed, over a certain energy scale $`\mathrm{{\rm Y}}(T)\mathrm{}^\psi `$. A simple model is to assume that energy barriers are independent random variables – this is the trap model . There is a limit where this model has a sharply defined behaviour: when the distribution of barriers is exponential, of mean $`T_c(\mathrm{})`$, then there is a well defined transition temperature which separates an equilibrium phase for $`T>T_c(\mathrm{})`$ where all configurations contribute more or less equally to the equilibrium partition function, from an aging phase for $`T<T_c(\mathrm{})`$ where the partition sum is dominated by a few metastable states only . This means that a small temperature change for $`T_c(\mathrm{})+\mathrm{\Delta }T/2`$ to $`T_c(\mathrm{})\mathrm{\Delta }T/2`$ completely changes the way the configuration space is explored. Correspondingly, a strong rejuvenation effect is expected upon crossing $`T_c(\mathrm{})`$ (see  for details). Furthermore, the realized probability to find the system in a metastable state of lifetime $`\tau `$ after a waiting time $`t_w`$ is, for $`T<T_c(\mathrm{})`$, given by :

$$P_{\mathrm{}}(\tau ,t_w)\frac{t_w^{x_{\mathrm{}}1}}{\tau ^x_{\mathrm{}}}\left(\frac{\tau }{t_w}\right)x_{\mathrm{}}=\frac{T}{T_c(\mathrm{})}<1,$$
(11)
where $`(u)`$ is a cut-off function, decaying as $`1/u`$ for large arguments. One can therefore obtain the a.c. susceptibility by averaging the simple two-level contribution Eq. (9) over the realized distribution of trapping times:

$$\chi ^{\prime \prime }(\omega ,t_w)=_{\omega ^1}^{\mathrm{}}𝑑\tau P_{\mathrm{}}(\tau ,t_w)\frac{\mathrm{exp}[t_w/\tau ]}{\omega \tau }(\omega t_w)^b,$$
(12)
with $`b=1x_{\mathrm{}}`$ . Therefore, the introduction of a large number of metastable states allows one to obtain both a strong rejuvenation effect at the transition temperature, and a realistic form for the decay of the a.c. susceptibility. Since there is a continuum hierarchy of length scales, each of which corresponding to a different transition temperature $`T_c(\mathrm{})`$, the problem is actually a multilevel trap model of the type studied in . The contributions of the different length scales are intertwined, and one expects rejuvenation to occur at all temperatures provided the corresponding $`\mathrm{}`$ is somewhat larger than the microscopic length. 333The a.c. susceptibility is however dominated at long times $`t_w`$ by the length scales such that $`x_{\mathrm{}}1`$, therefore explaining why $`b=1x_{\mathrm{}}`$ is found to be small ($`1/f`$ noise) and nearly temperature independent. The coexistence of strong rejuvenation and memory in this model has been numerically demonstrated very convincingly in , and confirm the above qualitative discussion.

There are many possible variants of the above trap model . The basic results still hold for more general barrier height distributions, with a weakly time dependent exponent $`x_{\mathrm{}}`$ . One can also consider long-range correlated energy landscapes, such as the one dimensional Sinai model where barriers typically grow with the position of the representative point . Many results on aging have been obtained for this model , which has a typical self-similar ‘valley within valley’ structure. Preliminary numerical studies also show that rejuvenation and memory effects are also present in the Sinai model . The mechanism at work in these landscape models is very close to the hierarchical picture first advocated by experimentalists to explain rejuvenation and memory : a reduction of temperature reveals finer details of the energy landscape within which the system must equilibrate.

### 4.3 More on ‘temperature chaos’

#### 4.3.1 Rejuvenation….

As mentioned in the introduction, the strong rejuvenation effect and the absence of cooling rate dependence have also been interpreted in terms of ‘temperature chaos’. Within the droplet model, the argument suggesting this behaviour is the following: the free energy of an excitation of length $`\mathrm{}`$ is small as a result of the compensation of its energy and entropy, both of them being much larger than $`\mathrm{{\rm Y}}\mathrm{}^\theta `$. Therefore, a small change of temperature $`\mathrm{\Delta }T`$ should ruin this subtle compensation on large length scale. This can be interpreted in terms of a complete re-shuffling of the dominant configuration beyond a certain overlap length $`\mathrm{}^{}`$ that diverges when $`\mathrm{\Delta }T0`$. The physical mechanism is actually related to the discussion of the previous subsection: since the shorter length scales have to reorganize when the temperature changes, the free energy of the larger length scales will be strongly affected. So the description of large length scales in terms of simple two (or multi) level systems with a fixed energy landscape is inappropriate: the ‘landscape’ itself is temperature dependent.

Using the value of the exponents in the droplet model in three dimensions, one naively estimates $`\mathrm{}^{}\mathrm{{\rm Y}}(T)/\mathrm{\Delta }T`$. Using the above results on $`\mathrm{{\rm Y}}(T)`$, one finds that for $`T=0.7T_c`$, $`T_c=15`$ K and $`\mathrm{\Delta }T=1`$ K, one should have $`\mathrm{}^{}1`$; therefore ‘chaos’ effects should be observable both numerically and experimentally. Since no sign of chaos was detected numerically for small systems , this suggests that numerical prefactors are perhaps large, and that actually the overlap length is, for practical purposes, larger than the numerically or experimentally relevant length scale. This scenario is supported by recent precise numerical studies of temperature chaos in a simpler problem – the pinning of a one dimensional domain wall (the directed polymer problem). In this case, the overlap length indeed scales as expected with $`\mathrm{\Delta }T`$, although rather large sizes are needed to observe the effect. The heuristic scaling arguments for chaos  can furthermore, for this simplified problem, be made more precise , and suggest that the decorrelation beyond the overlap length only decays as a power law rather than an exponential . The need to go to very large length scales appears to be true also for the 2D spin-glass . Therefore, even if temperature chaos appears to be absent in the sk model, it is possible that this effect exists in three dimensions, although perhaps on very large length scales.

On the other hand, as we discussed above, a hierarchical landscape rejuvenation appears sufficient to account for most of the experimental data. The main distinction is the existence of a characteristic overlap length. In this respect, one reason to believe that ‘strong’ temperature chaos is perhaps not relevant is that rejuvenation appears to be a small scale, rather than large scale, phenomenon. When the temperature is changed sufficiently to induce partial rejuvenation, the aging a.c. susceptibility is made up of two contribution, a short time, rejuvenation part, and a long time contribution which is the continuation of aging at the first temperature, with a shifted effective age to account for the modification of time scales with temperature. This suggests that the length scale $`\mathrm{}_1`$ built at $`T_1`$ actually continues to grow at $`T_2`$, although by definition, since some rejuvenation is observed, one should be in a regime where $`\mathrm{}^{}<\mathrm{}_1`$. The same effect is seen in zfc (zero field cooled) experiments and in numerical simulations : after a short transient, the length scale that grew at $`T_1`$ continues to grow (albeit at a different rate) at $`T_2`$.

Another argument against the relevance of an overlap length is provided by the quantitative analysis of the effect of a small temperature cycle on the a.c. susceptibility. If the length scales $`<\mathrm{}^{}`$ are unaffected, one expects that the effective initial age $`t^{}`$ of the a.c. susceptibility at $`T_1\mathrm{\Delta }T`$ is such that $`t(\mathrm{}^{},T_1\mathrm{\Delta }T)=t^{}`$ (provided that the size of the domains grown at $`T_1`$ is larger than the overlap length, i.e. $`\mathrm{}_1>\mathrm{}^{}`$). We have extracted $`t^{}`$ from a set of experiments where $`t_{w1}`$ and $`T_1`$ are fixed, and $`\mathrm{\Delta }T`$ is varied. For CrIn, we find that $`t^{}`$ behaves as $`\mathrm{\Delta }T^a`$, with $`a1`$ for $`T/T_c=0.72`$ and $`a1.5`$ for $`T/T_c=0.84`$. This is incompatible with the assumption that $`\mathrm{}^{}\mathrm{\Delta }T^1`$, which would lead to an extremely fast divergence of $`t^{}`$ for small $`\mathrm{\Delta }T`$, according to Eq. (2). A possibility would be that $`\mathrm{}^{}`$ behaves as $`\mathrm{\Delta }T`$ to a small negative power, but this violates bounds obtained within the framework of the droplet model .

#### 4.3.2 …and memory ?

Strong chaos does obviously explain complete rejuvenation for $`\mathrm{\Delta }T`$ sufficiently large. In order to be compatible with memory, one should argue, as above, that the length scale $`\mathrm{}_2`$ growing at $`T_2`$ remains somewhat smaller than $`\mathrm{}_1`$; the fast separation of time scales is then used to account for memory – see the discussion in . On the other hand, the magnetic field is known to lead to very strong ‘chaos’-like effects. For example, a small magnetic field cycle is sufficient to rejuvenate completely the system. Furthermore, chaos with magnetic field in the sk model was obtained long ago . An interesting idea is then to perform a simultaneous temperature and magnetic field cycle, with a $`\mathrm{\Delta }T`$ such that rejuvenation and perfect memory are observed in the absence of magnetic field changes. In the temperature chaos scenario, adding a magnetic field should not change anything since the overlap length is already very small. In the fixed landscape scenario, the magnetic field should have an effect since the large scale structures built at $`T_1`$ will couple to the field (as was discussed in subsection 3.2) and therefore speed up their dynamics at $`T_2`$. Therefore, some loss of memory should be observed, as indeed suggested by the experiments: see Figure 5.

Following this line of thought, we have also performed an experiment to test a very spectacular prediction of : a double temperature cycling where $`T_1T_2=T_1\mathrm{\Delta }T`$ is followed, after a certain time $`t_{w2}`$ at $`T_2`$ by a second small quench $`T_2T_3=T_2\mathrm{\Delta }T^{}`$. If the system is re-heated rapidly to $`T_1`$ without stopping at $`T_2`$, the dramatic prediction of  is that the memory effect should be destroyed. Only if the system is allowed to ‘take its breath’ at $`T_2`$ will memory be preserved. The problem is that ‘re-heating rapidly’ cannot be achieved experimentally, since the fastest achievable temperature ramps allows the system to stay at each temperature during a time which is huge in microscopic units. This might be enough (again due to the fast separation of time scales) to allow for memory conservation ! Our idea was then to use the magnetic field to prevent the system from tracing back its previous history, which is the crucial ingredient, in the chaos scenario, to preserve the memory (see the detailed discussion in ). We have therefore applied a magnetic field while the system is re-heated: the target state at $`T_2`$ is therefore completely scrambled. In spite of this, the loss of memory when returning at $`T_1`$ is the same as for a single quench procedure. Therefore, the spectacular effect predicted within the ‘strong chaos’ scenario fleshed out in  is not observed.

Finally, note that in the chaos scenario, it should be possible to choose $`\mathrm{\Delta }T`$, $`t_{w1}`$ and $`t_{w2}`$ in a positive temperature cycling such that $`\mathrm{}_2<\mathrm{}_1`$. In this case, memory should be preserved in a positive cycling experiment, even if some rejuvenation takes place at $`T_2`$. To our knowledge, this has never been observed, but this might again be due to the fact that the time scales separate very quickly. Taking the same value of parameters as in Fig. 2, we find that for $`T_2=0.77T_c`$, $`t_{w2}=10`$ sec., $`\mathrm{}_220`$. One would need to wait at least $`t_{w1}10^5`$ seconds (1.5 day) at $`T_1=0.7T_c`$ to reach the same length and allow memory to be preserved. Note that the fundamental assymmetry between small positive and small negative temperature jumps is probably the most clear cut difference between the ‘chaos’ scenario and the hiererchical landscape scenario. It would therefore be crucial to find an experimental situation where rejuvenation and memory in a positive $`\mathrm{\Delta }T`$ cycle should in principle be observed. It would be interesting to study this issue in the Ising-like sample studied in , where the separation of time scales is milder.

## 5 Conclusion – Open problems

In this paper, we have summarized the different puzzles raised by aging experiments of spin-glasses and their different interpretations. We try to reconcile the ‘real space’, droplet like pictures and the hierarchical pictures that have been proposed in the past. The basic ingredient is a strong separation of the time scales that govern the dynamics of the system on different length scales. Changing the temperature changes the length scale at which the system is observed, thereby allowing rejuvenation (that concerns short length scales) and memory (stored in long length scales) to coexist. We have shown that previous experiments can be reanalyzed in terms of vanishing energy barriers at the spin-glass transition, which in turn leads to the ‘super-activated’ behaviour observed in several experiments. We have argued that the power-law dependence of the coherence length with time might actually reflect a slow crossover form critical to activated dynamics. Finally, we have tried to distinguish between a hierarchical landscape rejuvenation, which is already present in simple multi-level systems like the Random Energy Model , from a more sophisticated ‘strong’ chaos that arises because large length scales free energies are renormalized by small length scale fluctuations . We have argued that most experiments can be accounted for without invoking the existence of an overlap length. Some specific predictions of the strong chaos scenario have been tested and our not borne out by our new experimental results. Nevertheless, we believe that this effect should exist on sufficiently large length scales, but these are perhaps out of reach both from numerical and experimental possibilities. It should finally be noted that rejuvenation and memory effects have also been observed in other, very different systems, such as pmma , where the relevance of temperature chaos is not clear, whereas a scenario based on multiscale dynamics is plausible.

From a theoretical point of view, the dynamics of mean field models corresponding to full replica symmetry breaking has been shown to exhibit rejuvenation and memory effects , and are actually very closely related to models of diffusion in self-similar landscapes such as the Sinai model , although the precise role of activation in these models is still rather obscure. However, these models are in principle incompatible with simple $`t/t_w`$ aging, as experimentally  and numerically  observed. It would be gratifying to understand whether or not these mean field models can correctly be interpreted in finite dimensions in terms of the simple trap models . This question of course has a far more general scope and relevant for other glassy systems as well.

Acknowledgements We are indebted to H. Yoshino for sharing with us many interesting ideas. We want to thank L. Berthier, E. Bertin, P. Doussineau, I. Giardina, J. Kurchan, J. Lamarcq, A. Levelut, O. Martin, M. Mézard and H. Takayama for interesting discussions. J.-P. B. wants to thank Harvard University for hospitality where part of this work was performed. Very stimulating discussions with D.S. Fisher, J.P. Sethna and R. da Silveira have to be acknowledged. This work owes a lot to the inspiring ideas of Daniel Fisher.


text:
# Collapsing Shells of Radiation in Higher Dimensional Spacetime and Cosmic Censorship Conjecture

## 1 Introduction

General relativity was formulated in spacetime with four dimensions ($`4D`$), of course. However, there are theoretical hints that we might live in a world with more dimensions. A generalization of general relativity to higher dimensions has been of considerable interest in recent times. It is believed that the underlying spacetime in the large energy limit of the Planck energy may have higher dimensions. At this level, all the basic forces of nature are supposed to unify and hence it would be pertinent in this context to consider solutions of the gravitational field equations in higher dimensions. Lately, there has been significant attention to studying gravitational collapse, in higher dimensions . llha et al  have generalized the Oppenhiemer-Snyder collapse model to higher dimensions.

Gravitational collapse continues to be a very important topic in gravitational research ever since Penrose  articulated the cosmic censorship conjecture some three decades ago. The cosmic censorship conjecture forbids the existence of naked singularities. In its strong version, the spacetime must be globally hyperbolic, whereas according to the weaker version, an event horizon must form during collapse and the singularity could be invisible to an asymptotic observer. The collapse of spherical matter in the form of dust or radiation forms shell focusing strong curvature singularities violating the cosmic censorship conjecture -.

The purpose of this letter is to study how these features get modified with extra dimensions. We study impolding radiation in a higher dimensional spacetime. The Vaidya metric , describing a spherically symmetric spacetime with radiation, is an exact solution of Einstein’s field equations. It has been generalized to higher dimensions  and we shall call it the higher dimensional Vaidya metric. We find that non self-similar higher dimensional Vaidya spacetime does admit strong curvature singularities in the sense of Tipler  for a sufficiently inhomogeneous collapse, providing an explicit counter example to cosmic censorship conjecture.

## 2 Higher Dimensional Vaidya Spacetime

The idea that spacetime should be extended from four to higher dimensions was introduced by Kaluza and Klein  to unify gravity and electromagnetism. Five dimensional $`(5D)`$ spacetime is particularly more relevant because both $`10D`$ and $`11D`$ super-gravity theories yield solutions where a $`5D`$ spacetime results after dimensional reduction . Hence, we shall confine ourselves to the $`5D`$ case.

The higher dimensional Vaidya spacetime which describes an implosion of radiation shells is

$$ds^2=(1\frac{m(v)}{r^2})dv^2+2dvdr+r^2d\mathrm{\Omega }^2$$
(1)
where $`d\mathrm{\Omega }^2=d\theta _1^2+sin^2\theta _1(d\theta _2^2+sin^2\theta _2^2d\theta _3^2)`$ is the metric of the 3-sphere, where $`v`$ is a null coordinate with $`\mathrm{}<v<\mathrm{}`$, $`r`$ is a radial coordinate with $`0r<\mathrm{}`$, and the arbitrary function $`m(v)`$ (which is restricted only by the energy conditions), represents the mass at advanced time $`v`$. The energy momentum tensor associated with eq. (1) can be written as

$$T_{ab}=\frac{3}{2r^3}\frac{dm}{dv}k_ak_b$$
(2)
with the null vector $`k_a`$ satisfying $`k_a=\delta _a^v`$ and $`k_ak^a=0`$. We have used units in which $`8\pi G=c=\mathrm{\hspace{0.33em}1}.`$ Clearly, for the weak energy condition to be satisfied, we require that $`dm/dv`$ be non negative. Thus the mass function is a non-negative increasing function of $`v`$. The Kretschmann scalar ($`K=R_{abcd}R^{abcd}`$, $`R_{abcd}`$ is the Riemann tensor) for the metric (1) takes the form

$$K=\frac{72m^2(v)}{r^8}$$
(3)
which diverges along $`r=0`$ establishing a scalar polynomial singularity. The Weyl scalar ($`C=C_{abcd}C^{abcd}`$, $`C_{abcd}`$ is the Weyl tensor) has the same expression as the Kretschmann scalar and thus the Weyl scalar also diverges whenever the Kretschmann scalar diverges and so the singularity is physically significant .

The physical situation here is that of a radial influx of a null fluid in an initially flat and empty region of the higher dimensional spacetime. For $`v<0`$ we have $`m(v)=\mathrm{\hspace{0.33em}0}`$, i.e., higher dimensional flat spacetime, and for $`v>T`$, $`dm/dv=\mathrm{\hspace{0.33em}0}`$, $`m(v)`$ is positive definite. The metric for $`v=0`$ to $`v=T`$ is higher dimensional Vaidya, and for $`v>T`$ we have the higher dimensional Schwarzschild solution. The first shell arrives at $`r=0`$ at time $`v=0`$ and the final at $`v=T`$. A central singularity of growing mass is developed at $`r=0`$. We shall now test whether future directed null geodesics terminate at the singularity in the past. If they do, the singularity is naked.

## 3 The Existence and Nature of Naked Singularities

In this section we investigate the existence of a naked singularity for the higher dimensional Vaidya spacetime. Let $`K^a=dx^a/dk`$ be the tangent vector to the null geodesic, where $`k`$ is an affine parameter. The geodesic equations, on using the null condition $`K^aK_a=0`$, take the simple form

$$\frac{dK^v}{dk}+\frac{m(v)}{r^3}(K^v)^2=0$$
(4)
$$\frac{dK^r}{dk}+\frac{1}{2r^2}\frac{dm}{dv}(K^v)^2=0$$
(5)
Following , we introduce

$$K^v=\frac{P}{r}$$
(6)
and, from the null condition, we obtain

$$K^r=\left(1\frac{m(v)}{r^2}\right)\frac{P}{2r}$$
(7)
The function $`P(v,r)`$ obeys the differential equation

$$\frac{dP}{dk}\left(1\frac{3m(v)}{r^2}\right)\frac{P^2}{2r^2}=0$$
(8)
In general, Eq. (8) may not yield an analytical solution. However, for our purpose the explicit solution of eq. (8) is not necessary. Radial null geodesics of the metric (1), by virtue of Eqs. (6) and (7), satisfy

$$\frac{dr}{dv}=\frac{1}{2}\left[1\frac{m(v)}{r^2}\right]$$
(9)
Clearly, the above differential equation has a singularity at $`r=0`$, $`v=0`$. The nature (a naked singularity or a black hole) of the collapsing solutions can be characterized by the existence of radial null geodesics coming out from the singularity. The character of the singularity depends on the exact form of $`m(v)`$. For example, with $`m(v)\lambda v^2`$ the spacetime is self-similar , admitting a homothetic Killing vector and singularities can be analyzed with ease. However, self-similarity is a strong geometric condition on the spacetime. It is therefore, of interest to us to examine more general forms of the function $`m(v)`$. Here, we construct specific examples which satisfy the weak energy condition but develop strong curvature naked singularities.

### Example-I

To analyze the nature of the singular point first we consider a higher dimensional analogue of Lake’s  solution, which would require

$$m(v)=\lambda v^2+f(v)$$
(10)
where $`\lambda `$ is a constant and $`f(v)=O(v^2)`$ as $`v0`$. It follows therefore that for $`0<\lambda 1/27`$, the singular point becomes an unstable node and the family of null geodesics meets the singularity with definite tangents. The possible values of the tangents are given by the roots of the characteristic for eq. (9)

$$\lambda \gamma _0^3\gamma _0^2+2=0$$
(11)
where

$$\gamma _0=\underset{r0v0}{lim}\gamma =\underset{r0v0}{lim}\frac{v}{r}$$
(12)
Thus gravitational collapse of null fluid in higher dimensions leads to a naked singularity if $`\lambda 1/27`$, and to formation of a black hole otherwise. The two positive roots of eq. (11), $`\gamma _0=2.21833`$ and $`5.69593`$, correspond to $`\lambda =1/50`$. For all such values, the singularity is naked. The degree of inhomogeneity of the collapse is defined as $`\mu 1/\lambda `$ (see ). We see that for a collapse sufficiently inhomogeneous, naked singularities develop. Comparison with the analogous $`4D`$ case shows that a naked singularity occurs for a slightly larger value of the inhomogeneity factor in higher dimensions. The global nakedness of the singularity can then be seen by making a junction onto the higher dimensional Schwarzschild spacetime. The Kretschmann scalar is $`K=a\lambda ^2/r^4`$, for some constant $`a`$. Thus as $`r0`$ the collapse forms a scalar polynomial curvature singularity.

The critical direction associated with the node is given by

$$r=\mu v+g(v)$$
(13)
where $`\mu =1/\gamma _0`$ and

$$f=(12\mu )g(2\mu v+g)g^{}(\mu v+g)^2$$
(14)
For a spacetime to be self-similar we require that $`g=0`$. The radial null geodesics given by Eq. (13) is the Cauchy horizon associated with the node, which is a strong curvature singularity (see below). One can see that the weak energy condition is satisfied for

$$2\lambda v(2\mu 1)2\mu g+2g^{}(3\mu +g^{}1)(\mu v+g)+g^{\prime \prime }(\mu v+g)^2$$
(15)
and the Cauchy horizon expanding for

$$g^{}>\mu $$
(16)
the function $`g`$ being restricted by (15) and (16).

The strength of a singularity is an important issue because there have been attempts to relate it to stability . A singularity is termed gravitationally strong or simply strong, if it destroys by crushing or stretching any object which falls into it. Along a null geodesic affinely parameterized by $`k`$, let $`\psi =R_{ab}K^aK^b`$ where $`R_{ab}`$ is the Ricci tensor and where the geodesic terminates at $`\lambda =0`$. We consider the following condition:

$$\underset{k0}{lim}k^2\psi >0$$
(17)
which is equivalent to the termination of a geodesic in a strong-curvature singularity in the sense of Tipler (cf. ). Eq. (17), with the help of eqs. (2), (10) and (6), can be expressed as

$$\underset{k0}{lim}k^2\psi =\underset{k0}{lim}\frac{3}{2r}\frac{dm}{dv}(\frac{kP}{r^2})^2$$
(18)
Our purpose here is to investigate the above condition along future directed null geodesics coming out from the singularity. Using the fact that, as the singularity is approached, $`k0`$ and $`r0`$ and using L’Hôpital’s rule, we first observe that $`lim_{k0}kP/r^2=2/1+\lambda \gamma _0^2`$ for $`P_0=\mathrm{}`$ and $`lim_{k0}kP/r^2=1/1\lambda \gamma _0^2`$ for $`P_0\mathrm{}`$ ($`P_0=lim_{k0}P`$) and hence eq. (18) gives

$$\underset{k0}{lim}k^2\psi =\frac{12\lambda \gamma _0}{(1+\lambda \gamma _0^2)^2}>0forP_0=\mathrm{}$$
(19)
and

$$\underset{k0}{lim}k^2\psi =\frac{3\lambda \gamma _0}{(1\lambda \gamma _0^2)^2}>0forP_0\mathrm{}$$
(20)
where $`\gamma _01/\sqrt{\lambda }`$. Thus along radial null geodesics, the strong curvature condition is satisfied.

Recently, Nolan  gave an alternative approach to check the nature of singularities without having to integrate the geodesics equations. It was shown in  that a radial null geodesic which runs into $`r=0`$ terminates in a gravitationally weak singularity if and only if $`\dot{r}`$ is finite in the limit as the singularity is approached (this occurs at $`k=0`$), the over-dot here indicates differentiation along the geodesics. So assuming a weak singularity, we have

$$\dot{r}d_0rd_0k$$
(21)
Using the asymptotic relationship above and that $`m(v)\lambda v^2`$ (as $`k0`$), the geodesic equations yield

$$\frac{d^2v}{dk^2}\delta k^1$$
(22)
where $`\delta =\lambda \gamma _0^4d_0=`$ a non-zero const., which is inconsistent with $`\dot{v}d_0\gamma _0`$, which is finite. Since the coefficient $`\delta `$ of $`k^1`$ is non-zero, the singularity is gravitationally strong .

### Example-II

We now construct a Joshi and Dwivedi  type solution for our higher dimensional spacetime. We choose the mass function as:

$$m(v)=\beta ^2v^{2\alpha }(12\beta \alpha v^{\alpha 1})$$
(23)
where $`\alpha >1`$ and $`\beta >0`$ are constants. This is a representative class of a more general problem $`m(v)v^n`$ ($`n>2`$). Clearly, $`\alpha =1`$ corresponds to a self-similar model. The null radiation shells start imploding at $`v=0`$ and the final shell arrives at $`v=T`$. To ensure that the weak energy condition is satisfied, i.e., $`dm/dv0`$, $`T`$ must satisfy,

$$T^{\alpha 1}<\frac{1}{\beta (3\alpha 1)}$$
(24)
which also guarantees that $`m(v)>0`$. Using (9), an outgoing radial null geodesic for the mass function (23) meeting the singularity $`v=0`$, $`r=0`$ in the past is given by

$$r=\beta v^\alpha $$
(25)
This integral curve meets the singularity with tangent $`r=0`$ and the singularity is naked. Since $`dr/dv>0`$ with increasing $`v`$, the null geodesics (25) escape to infinity and the singularity is globally naked. It is seen that the condition $`r^2>m(v)`$ is satisfied along the trajectories.

It is easy to see that the Kretschmann scalar diverges along singular null geodesics meeting the singularity in the past in the approach to the singualrity, establishing the presence of a scalar polynomial singularity. However, unless $`m(v)v^2`$ in the approach to the singularity, the strong curvature condition is not satisfied along radial null geodesics. In the analogous situation in $`4D`$, it was found that the singularities are strong-curvature only for mass functions which are initially linear functions of $`v`$ .

## 4 Conclusions

A rigorous formulation and proof for either version of the cosmic censorship conjecture is not available. Hence, examples showing the occurrence of naked singularities remain important and may be valuable if one attempts to formulate the notion of the conjecture in precise mathematical form. The Vaidya metric in the $`4D`$ case has been extensively used to study the formation of naked singularities in spherical gravitational collapse . We have extended this study to a higher dimensional Vaidya metric, and found that strong curvature naked singularities do arise for slightly higher values of the inhomogeneity parameter and only for functions which are initially quadratic functions of the advanced time, i.e., $`m(v)v^2`$. We have checked for naked singularities to be gravitationally strong by the method in  and by an alternative approach proposed by Nolan  as well, and both are in agreement. The models constructed here are not self-similar in general, and as a special case self-similar models arise. Also, it is straight forward to extend the above analysis for non radial causal curves. Whereas we confined our analysis to $`5D`$, there is no reason why it cannot be extended to spacetime of any dimensions ($`n4`$).

In conclusion, this offers a counter example to the cosmic censorship conjecture.

Acknowledgment: Two of the authors (SGG and RSV) would like to thank IUCAA, Pune (INDIA) for hospitality. Thanks are also due to Brien C. Nolan for helpful correspondence.