text:
# Scars of the Wigner function.

##

Sixteen years have passed since Heller  detected scars of periodic orbits in individual eigenfunctions of chaotic systems. Explanations in terms of wave packets  or the semiclassical Green function  do predict an enhancement of intensity near the projection of a Bohr-quantized periodic orbit. However, such theories apply to a collective superposition of states near this quantization condition. In spite of some quite striking visual evidence, the issue of scarring for individual states, is confused by the fact that different branches of the same periodic orbit may overlap so that their contributions interfere in the position space and this must be superposed on the expected random wave background for the chaotic state. This has led to the need to make quantitative assessments of scarring strength . We here show that the phase space picture can provide sharper qualitative evidence of the influence of a periodic orbit on its Wigner functions.

The old conjecture of Berry and Voros  that the Wigner functions for eigenstates of chaotic Hamiltonians are concentrated on the corresponding classical energy shells is widely disseminated , as it is compatible with Schnirelman’s theorem . Nonetheless, no such restriction arises in the more recent semiclassical theory for the mixture of these states over a narrow energy window, i.e the spectral Wigner function . Indeed, the computations presented in this letter show that individual Wigner functions, as well as their energy average can oscillate with large amplitudes deep inside the energy shell.

At points $`𝐱=(𝐪,𝐩)`$ inside the energy shell, the semiclassical spectral Wigner function is :

$$W(𝐱,E,\epsilon )=\underset{j}{}A_j(𝐱,E)e^{\epsilon t_j/\mathrm{}}\mathrm{cos}\left[\frac{S_j(𝐱,E)}{\mathrm{}}+\gamma _j\right].$$
(1)
The sum is over all the trajectory segments on the $`E`$-shell with endpoints, $`𝐱_{j\pm }`$, centered on $`𝐱`$, as sketched in Fig.1. The action is merely the symplectic area, $`S_j=𝐩𝑑𝐪`$, for the circuit taken along the orbit from $`𝐱_j`$ to $`𝐱_{j+}`$ and closed by the chord $`𝝃_j(𝐱)`$. The time of traversal for the stretch along the trajectory is $`t_j`$ and $`\epsilon `$ is the width of the energy window over which we average individual Wigner functions $`W_n(𝐱)`$, i.e.,

$$W(𝐱,E,\epsilon )=(2\pi \mathrm{})^{D/2}\underset{n}{}\frac{\epsilon /\pi }{(EE_n)^2+\epsilon ^2}W_n(𝐱),$$
(2)
where $`D`$ is the dimension of the phase space $`(𝐪,𝐩)`$. We shall not be concerned with the Maslov phase $`\gamma _j`$, nor with the amplitude $`A_j(𝐱,E)`$, except to note that these are purely classical quantities that vary smoothly with $`𝐱`$ inside the shell, as compared with the high frequency oscillations of the cosine factor.

The highly oscillatory nature of the spectral Wigner function inside the energy shell is mostly washed out on projection, i.e. for true probability densities. However, the oscillations are unavoidable in the study of interference effects.

The energy shell itself is a caustic, because the chord $`𝝃0`$ as $`𝐱`$ approaches the shell. In this limit, the contributing trajectories either shrink to a point, or they are very close to a periodic orbit. The modification of (1) leads to the Berry theory of scars , futher refined in , for evaluation points $`𝐱`$ close to the energy shell.

The point here is that large open segments of periodic orbits can be constructed by adding multiple windings to a primitive, small segment of a periodic orbit. If conditions for phase coherence, to be discussed, are satisfied, we can then obtain a scar deep inside the energy shell. This scar is located at the two-dimensional central surface constructed by the centers of all the chords with endpoints on the periodic orbit. For the case where the phase space dimension $`D=2`$, the central surface is the phase space region enclosed by the convex hull of the orbit that here is the energy shell ***Note that if the orbit is convex the central surface is just the phase space region inside it.. For $`D4`$, the two-dimensional central surfaces are still bounded by one-dimensional periodic orbits that lie within the higher dimensional energy shell. In the simplest case where the orbit is plane and convex, the central surface coincides with the part of the plane within the orbit just as in the case where $`D=2`$. If the orbit is not plane, then the central surface will be more complicated and it may exhibit singularities and self-intersections. In general it will not coincide everywhere with the invariant surface formed by changing continuosly the energy of this periodic orbit.

We emphasize that these scars of the Wigner function that correspond to a very recognizable oscillatory pattern have special localization properties in phase space of dimension $`D4`$. Previous explorations of phase space representations of eigenstates have been too blunt to discern these fine features. The work of Feingold et. al.  involves a projection from a three-dimensional section, which is fine for the Husimi function, but washes away the details of an oscillatory Wigner function. On the other hand, Agam and Fishman  restrict their investigation to the neighbourhood of the orbit and, hence, to the energy shell. This is just the edge of the central surface.

In this letter we tested our prediction for the simplest case of an unstable periodic orbit lying in a two dimensional plane. Accidentally this plane is a symmetry plane, however, we have checked that the enhanced amplitude of the spectral Wigner function accompanies the distortion of the orbit due to a nonlinear canonical transformation. Since the Wigner function is not invariant in this case, we thus obtain an essentially new system, in which the new central surface is not a symmetry surface. Thus, our scarring effect cannot be attributed to a special property of the symmetry plane. This case is essentially different from recent studies of higher dimensional systems where the scars arise over marginally stable invariant planes that are indeed special .

The chord structure for segments of a periodic orbit is sketched in Fig.2. This is similar to a system with $`D=2`$, for which Berry  showed that there is only one chord for most points $`𝐱`$. However, we must now distinguish between the chords $`𝝃_{in}=𝝃_{out}`$ and the sum in (1) includes each different winding of the periodic orbit for primitive orbit segments, “in” and “out”, associated with the two respective chords. Evidently all these “in” and “out” contributions build up if the energy is close to one of the values for which this periodic orbit is Bohr-quantized. We show in  that the condition for these two contributions to be in phase is the same as that for a maximal contribution to the Gutzwiller trace formula . Therefore, near the quantization energy, the sum of all the contributions in (1), owing to chords in the periodic orbit, can be approximated by an expression analogous to Berry’s simplest semiclassical aproximation  for the Wigner function in $`D=2`$ systems, but now evaluated over points $`𝐱`$ in the two-dimensional central surface. Hence, in this case, we have an enhanced contribution to the spectral Wigner function with the phase $`S_{in}(𝐱)/\mathrm{}`$ and the Maslov phase $`\gamma =0`$ . The successive contours $`S_{in}(𝐱)=`$constant thus determine rings of constant phase along the central surface. Therefore, we predict rings of positive and negative amplitude superimposed on the background of contributions from uncorrelated trajectory segments that also contribute to the spectral Wigner function. The edge of this system of rings along the orbit itself is the object of the Berry theory , though it does not deal with multiple windings.

We have tested this prediction for the coupled nonlinear oscillator ($`D=4`$)

$$\text{H}(𝐪,𝐩)=\frac{(p_1^2+p_2^2)}{2}+0.05q_1^2+\left(q_2\frac{q_1^2}{2}\right)^2,$$
(3)
sometimes referred to as the “Nelson Hamiltonian”. The classical periodic orbits together with the topology of their families in the energy.vs.period plot has been studied in  and the wave functions and their scars were discussed in . Evidently the plane $`q_1=p_1=0`$ is classically invariant, so we obtain an isolated periodic orbit on this plane for each energy (the vertical family). This is then, the simplest case where the central surface is flat - the central plane -. The restriction of the Hamiltonian to this plane defines a harmonic oscillator, so there is a single chord with endpoints in the periodic orbit (except for sign) for each point $`(p_2,q_2)`$ inside the orbit, except at $`q_2=p_2=0`$.

We computed the wave functions $`𝐪|n`$ for the eigenstates of (3) within the range of energies $`0.821E0.836`$ taking $`\mathrm{}=0.05`$ and using a basis of eigenfunctions $`\varphi _{n_1}^{(q_1)}(q_1)\varphi _{n_2}^{(q_2)}(q_2q_1^2/2)`$ , in the same way as in . The Wigner function for each state was then calculated directly by the double integral over $`𝐐=(Q_1,Q_2)`$:

$$W_n(𝐪,𝐩)=\frac{1}{(2\pi \mathrm{})^2}𝑑𝐐𝐪+𝐐/2|nn|𝐪𝐐/2e^{i\frac{𝐩𝐐}{\mathrm{}}}.$$
(4)
In Fig.3 we display the average over Wigner functions for the energy window chosen, surrounding the 11’th Bohr energy for the vertical orbit that has Maslov index $`3`$. This is more convenient than smoothing with the Lorentzian window in (2) and should produce essentially the same results.

As an example, Fig.4 shows some of the individual Wigner functions in this window for $`𝐱`$ in the $`(q_2,p_2)`$ plane. Immediatly, we notice that all these Wigner functions are remarkably regular, roughly in the region $`q_2<0`$. Indeed, Fig.3 coincides in this region with the Wigner function of the restricted Hamiltonian, also with respect to the precise phase of the oscillations. The individual Wigner functions follow the same phase contours, but there is a varying phase shift. If we extrapolate the semiclassical theory to energy smoothings of the order of the energy spacing, the conclusion is that there are no orbit segments, other than the periodic orbit itself, contributing to the Wigner functions on this surface up to the Heisenberg time. This is more dramatic for the state $`n=305`$ whose Wigner function on all the central surface is almost the same of that of the restricted Hamiltonian (see Fig.4). It is important to note that only this state seems to have a visual scar of the vertical orbit in position space, although superimposed to a random wave background.

For $`q_2>0`$, we have found other orbit segments that do not lie on the central plane, but which do contribute to the spectral Wigner function on it. Indeed, this is the case for segments of both the symmetric periodic orbits shown in Fig.6 projected onto position space. Needless to say, the multiple windings of these periodic orbits could produce total contributions of the same order as the plane orbit, if they are nearly Bohr-quantized. This supplies a qualitative explanation for the existence of interference patterns to the left of the Wigner functions in Fig.3 and 4. Notice that the second of the orbits in Fig.6 reaches into a region with $`q_2`$ greater than is attained by the plane periodic orbit. So we account for Wigner functions reaching outside the energy shell (which coincides with the periodic orbit on the central plane). Note that these are only particular examples of the many symmetric periodic orbits whose central surface intersects the central plane along lines with $`q_2>0`$ . Therefore, it will be hard to make quantitative predictions in this region in contrast to that where $`q_2<0`$, for which no such orbits of short period were found. The enhanced scar pattern of the Wigner function of the state $`n=305`$ over all the central plane evidently washes out other possible contributions.

The central surface of the harmonic orbit considered in this work coincides with an invariant plane of the Hamiltonian (3) which is a reflection symmetry plane of the system. This might induce to the conclusion that the symmetry is the actual responsible for the scarring observed. We show that this is not the case by applying a nonlinear canonical transformation, described in , deforming the central plane $`q_1=p_1=0`$ into a new invariant surface (Fig.5 a). Now, the new central surface of the distorted periodic orbit no longer coincides everywhere with the invariant surface (Fig.5 b) losing in this way any symmetry. We calculated the Wigner functions for this new system, over the central surface of the periodic orbit, for eigenstates within the same range of energies surrounding the 11’th Bohr level. Both the average and the individual Wigner functions, projected over the $`(q_2,p_2)`$ plane, display the same features as those shown in Fig.3 and Fig.4 respectively. Reference  also considers the thickness of scar surfaces.

Our main point is that individual trajectory segments do play a role in the spectral Wigner function. Indeed, phase coherent contributions from a periodic orbit can add up to scars in the form of concentric rings deep inside the energy shell. Eventhough the Berry-Voros hypothesis  ties in more intuitively with the concept of quantum ergodicity, it should be noted that the chord picture of Wigner functions does not contradict Shnirelman’s theorem and subsequent exact results , because the oscillations make negligible contributions to integrals with smooth functions.

It would seem that the structure of interfering chords, from which the spectral Wigner function is built up, could only generate a messy picture where scars would be even harder to recognize than in wave functions for which quantitative methods are needed. We have now shown that the converse can be true in the simplest case, i.e. that the relatively low dimension of the central surface can allow it to miss, over an appreciable region, the contribution of chords from other orbits. In this region we obtain a regular ring structure, not only for the energy average, but even for the pure Wigner functions shown in Fig.4. Though we cannot yet predict why the relative weight of this particular periodic orbit should vary so markedly among the states close to the Bohr-quantized energy, we are now able to recognize in Fig.4 a remarkable example of an individual scar.

###### Acknowledgements.

This work was financied by CNPq-CLAF, FAPERJ and Pronex-MCT.


format: text

text:
# Resistivity and 1/𝑓 Noise in Non-Metallic Phase Separated Manganites

## I Introduction

Recent experimental  and theoretical  papers provide a strong evidence for the existence of the phase separated state in the perovskite manganites with the colossal magnetoresistance (CMR). Both experiment and theory demonstrate that under certain conditions the material in the phase separated state consists of small ferromagnetic metallic droplets, or magnetic polarons, embedded into the insulating antiferromagnetic (AFM) matrix. The droplet in the ground state contains one charge carrier (electron or hole) confined in a potential well of ferromagnetically ordered local spins. In Refs.  the model allowing the estimation of the droplet radius and droplet concentration was proposed. According to this model, the radius of the droplet $`a`$ is defined from the minimization of the energy: $`Et(\pi d/a)^2+JS^2(4\pi /3)(a/d)^3`$, where the first term is related to the kinetic energy of electron delocalization in the bubble with the radius $`a`$, and the second term corresponds to the loss in the energy of the Heisenberg AFM exchange due to ferromagnetic ordering of local spins $`S`$ inside the bubble. The minimization of the energy with respect to droplet radius yields the following estimate for $`a`$ in $`3D`$ case: $`(a/d)(t/JS^2)^{1/5}`$, where $`J`$ is the AFM Heisenberg exchange, $`t`$ is the bandwidth, and $`d`$ is the intersite distance. The number of charge carriers is proportional to the electron (hole) doping $`x`$. Volume concentration of metallic droplets increases with $`x`$ and with the decrease of the temperature since the droplet radius $`a`$ decreases with temperature growth. As a result, at some critical concentration of holes $`x_c(d/a)^3`$ the droplets start to overlap and the percolation metal-insulator transition occurs in the system . In Ref.  it was observed that heating and cooling of the perovskite manganites in the phase separated state is accompanied by strong hysteresis in the magnetization and the resistance. In addition, the giant $`1/f`$ noise was measured in these experiments. The noise power is very high even far from the percolation threshold and drastically increases in its vicinity. The noise spectrum is close to $`1/f`$ form in the $`1`$$`1000`$ Hz frequency range.

In this paper, we calculate conductivity, magnetoresistance, and noise spectral power of the system in the phase separated state. The calculations are based on the results of Refs.  and on a simple model for the tunneling conductivity of the material accounting for electron jumps from one polaron to another. The concentration range not too close to the percolation transition is considered.

We use the terms (magnetic) polarons and droplets interchangeably throughout the paper.

## II Conductivity

Let us consider an insulating antiferromagnetic sample of volume $`V_s`$ in electric field $`𝑬`$. The total number of magnetic polarons in the volume is $`N`$, and thus their spatial density is $`n=N/V_s`$. As it was mentioned before, the number of polarons is assumed to be equal to the number of charge carriers introduced by doping. Neglecting the conductivity of the insulating phase, we assume that charge carriers are only located within the droplets. The charge transfer can thus occur either due to the motion of the droplets or due to the electron tunneling. The former mechanism is less effective: Indeed, the motion of a droplet is accompanied by a considerable rearrangement of the local magnetic structure, which results in the big effective mass of magnetic polarons. In addition, the droplets are expected to be easily pinned by crystal lattice defects. Thus, it is realistic to assume that the charge transport is essentially due to electron transitions between the droplets.

A magnetic polaron in the ground state contains one electron. As a result of a tunneling process, droplets with more than one electron are created, and some droplets become empty (the lifetime of such excitations is discussed in the end of this Section). If the energy of an empty droplet $`E(0)`$ is taken to be zero, the energy of a droplet with one electron can be estimated as $`E(1)t(d/a)^2`$. This is essentially the kinetic energy of an electron localized in the sphere of the radius $`a`$. In the same way, the energy of two-electron magnetic polaron $`E(2)2E(1)+U`$, with $`U`$ being the interaction energy of the two electrons. In all these estimates, we have disregarded the surface energy, which is expected to be small . Thus, $`E(2)+E(0)>2E(1)`$, and the creation of two-electron droplets is associated with the energy barrier of the order of $`AE(2)2E(1)U`$. It is clear that the interaction energy $`U`$ of two electrons in one droplet is determined mainly by the Coulomb repulsion of these electrons, hence $`Ae^2/ϵa`$, where $`ϵ`$ is the static dielectric constant, which in real manganites can be rather large ($`ϵ20`$). We assume below that the mean distance between the droplets is $`n^{1/3}a`$ (the droplets do not overlap). Then, $`A`$ is larger than the average Coulomb energy $`e^2n^{1/3}/ϵ`$. Since the characteristic value of the droplet radius $`a`$ is of the order of 10 Å, we have $`A/k_B1000K`$ and $`A>k_BT`$ in the case under study. In the following, we assume that the temperature is low, $`Ak_BT`$, and we do not consider a possibility of the formation of the droplets with three or more electrons. Even in the case when these excitations are stable, it can be shown that far from the percolation threshold the strong Coulomb interaction suppresses their contribution to the conductivity (giving rise only to the next order terms with respect to $`\mathrm{exp}(A/k_BT)`$).

Let us denote the numbers of single-electron, two-electron, and empty droplets as $`N_1`$, $`N_2`$, and $`N_3`$, respectively. According to our model, $`N_2=N_3`$, $`N_1+2N_2=N`$, and $`N`$ is constant. Before turning to conductivity, we evaluate the thermal averages of $`N_1`$ and $`N_2`$. To this end, we note that the number $`P_N^m`$ of possible states with $`m`$ two-electron droplets and $`m`$ empty droplets equals $`C_N^mC_{Nm}^m`$, with $`C_N^m`$ being the binomial coefficients. Since the created pairs of droplets are independent, we write the partition function in the form

$$Z=\underset{m=0}{\overset{N/2}{}}P_N^m\mathrm{exp}(m\beta ),\beta =A/k_BT.$$
(1)
Though the sum can be evaluated exactly and expressed in terms of the Legendre polynomials for arbitrary $`N`$, it is more convenient to use the Stirling formula for the factorials and the condition that the sample is macroscopic, $`N1`$. Approximating the sum by an integral,

$`Z`$ $`=`$ $`{\displaystyle _0^{N/2}}𝑑m`$
$`\times `$ $`\mathrm{exp}\left[m\beta N\mathrm{ln}\left(1{\displaystyle \frac{2m}{N}}\right)+2m\mathrm{ln}\left({\displaystyle \frac{N}{m}}2\right)\right],`$
calculating it in the saddle-point approximation, and subsequently evaluating in the same way the statistical average of $`N_2`$,

$$\overline{N}_2=Z^1\underset{m=0}{\overset{N/2}{}}mP_N^m\mathrm{exp}(m\beta )=\frac{}{\beta }\mathrm{ln}Z,$$
(2)
we easily obtain

$`\overline{N}_2`$ $`=`$ $`N\mathrm{exp}(A/2k_BT),`$ (3)
$`\overline{N}_1`$ $`=`$ $`N2\overline{N}_2=N\left[12\mathrm{exp}(A/2k_BT)\right].`$ (4)
Now we calculate the conductivity. Within the framework of the proposed model the electron tunneling occurs via one of the four following processes illustrated in Fig. 1,

(i) In the initial state we have two droplets in the ground state, and after tunneling in the final state we have an empty droplet and a droplet with two electrons;

(ii) An empty droplet and a two-electron droplet in the initial state transform into two droplets in the ground state (two droplets with one electron);

(iii) A two-electron droplet and a single-electron droplet exchange their positions by transferring an electron from one droplet to the other;

(iv) An empty droplet and a single-electron droplet exchange their positions by transferring an electron from one droplet to the other.

In the linear regime, all these processes contribute to the current density $`j`$ independently, $`j=j_1+j_2+j_3+j_4`$. The contributions of the first two processes read

$$j_{1,2}=en_{1,2}\underset{i}{}v_{1,2}^i,$$
(5)
where $`n_{1,2}=N_{1,2}/V_s`$ are the densities of the single- and two-electron droplets, and $`\mathrm{}`$ stands for statistical and time averages. The appearance of the factors $`n_{1,2}`$ reflects the fact that the electron tunnels from a single-electron droplet (process (i)) or two-electron (ii) droplet. The summation in Eq. (5) is performed over all magnetic polarons the electron can tunnel to — one-electron droplets for the process (i) and empty droplets for the process (ii). Finally, the components of average electron velocity $`v_{1,2}^i`$ along the direction of the electric field are obviously found as

$$\underset{i}{}v_{1,2}^i=\underset{i}{}\frac{r^i\mathrm{cos}\theta ^i}{\tau _{1,2}(r^i,\theta ^i)},$$
(6)
where $`r^i`$ and $`\theta ^i`$ are the electron tunneling length (the distance between the droplets) and the angle between the electric field and the direction of motion, respectively, and $`\tau _{1,2}(r^i,\theta ^i)`$ are characteristic times associated with the tunneling processes. The relation between $`\tau _1(r,\theta )`$ and $`\tau _2(r,\theta )`$ can be found from the following considerations. Near the equilibrium, the number of two-electron droplets, excited per unit time, equals to the number of the decaying two-electron droplets. We thus have the detail balance relation,

$$\frac{\overline{N}_1^2}{\tau _1(r,\theta )}=\frac{\overline{N}_2^2}{\tau _2(r,\theta )},$$
(7)
where we have taken into account that the probability of the formation of a two-electron droplet is proportional to the total number $`N_1`$ of the single-electron states multiplied by the number of available hopping destinations, which also equals $`N_1`$. Similarly, the probability of decay of a two-electron droplet is proportional to $`N_2N_3=N_2^2`$. Eq. (7) implies $`\tau _2(r,\theta )=\tau _1(r,\theta )\mathrm{exp}(A/k_BT)`$. We write then the conventional expression for the tunneling times  in the following form,

$$\tau _{1,2}(r,\theta )=\omega _0^1\mathrm{exp}\left(\frac{r}{l}\pm \frac{A}{2k_BT}\frac{eEr\mathrm{cos}\theta }{k_BT}\right),$$
(8)
where $`l`$ and $`\omega _0`$ are the characteristic tunneling length and magnon frequency, and we have taken into account the contribution of the external electric field to the tunneling probability.

To perform the averaging, we assume that the centers of the magnetic polarons are randomly positioned in space and the average distance $`n^{1/3}`$ between them is much larger than the droplet radius $`a`$. Both assumptions seem to be perfectly justified far below the percolation threshold. Then the averaged sum in Eq. (5) is essentially the space average of $`v^i`$, multiplied by the number of droplets available for hopping ($`N_1`$ for the process (i) and $`N_2`$ for the process (ii)). Expanding in $`eEl/k_BT1`$, we obtain

$`{\displaystyle \underset{i}{}}v_{1,2}^i`$ $`=`$ $`{\displaystyle \frac{eE\omega _0}{k_BT}}N_{1,2}e^{A/2k_BT}r^2\mathrm{cos}^2\theta e^{r/l}_V,`$ (9)
$`\mathrm{}_V`$ $`=`$ $`V_S^1{\displaystyle \mathrm{}d^3𝒓}.`$ (10)
In Eq. (9), electric field is outside the averaging. Rigorously speaking, this means that the characteristic hopping length $`l`$ is larger than the interdroplet distance $`n^{1/3}`$ and our approach is valid only when the droplet concentration is not too small.

Substituting Eq. (9) into Eq. (5) and performing the integration, we find

$$j_{1,2}=\frac{32\pi e^2E\omega _0l^5n_{1,2}^2}{k_BT}\mathrm{exp}(A/2k_BT).$$
(11)
In the processes (iii) and (iv) the free energy of the system is not changed after the tunneling, and we write the characteristic times as

$$\tau _{3,4}(r,\theta )=\omega _0^1\mathrm{exp}(r/leEr\mathrm{cos}\theta /k_BT).$$
(12)
The contribution of these two processes to the current is calculated similarly to that of (i) and (ii). For the process (iii) the number of magnetic polarons from which the electron may tunnel, is $`N_2`$, whereas the number of accepting droplets is $`N_1`$. In the same way, for the process (iv) these numbers are $`N_1`$ and $`N_3=N_2`$, respectively. Consequently, the factors $`n_{1,2}^2`$ in Eq. (11) are replaced by $`n_1n_2`$,

$$j_{3,4}=\frac{32\pi e^2E\omega _0l^5n_1n_2}{k_BT}.$$
(13)
From Eqs. (11) and (13) we now obtain the dc conductivity $`\sigma =j/E`$,

$$\sigma =\frac{32\pi e^2\omega _0l^5}{kT}\left[2n_1n_2+n_1^2e^{A/2k_BT}+n_2^2e^{A/2k_BT}\right].$$
(14)
In this Section we are only interested in the average conductivity; fluctuations lead to the appearance of noise and are considered in Section IV. Using Eq. (3), we find that all the four processes illustrated in Fig. 1 give identical contributions to the conductivity; for $`Ak_BT`$ the average conductivity (for which we retain the notation $`\sigma `$) reads

$$\sigma =\frac{128\pi e^2n^2\omega _0l^5}{k_BT}\mathrm{exp}(A/2k_BT).$$
(15)
We see that the conductivity increases with temperature as $`\sigma (T)T^1\mathrm{exp}(A/2k_BT)`$, which is typical for tunneling systems (see e.g. Ref. ).

At this point, let us discuss the applicability range of our model. The essence of our picture is the existence of different types of droplets. Of course, only single-electron droplets are stable. Obviously, an empty droplet decays during the time of the order of $`1/\omega _0`$. On the other hand, following the above discussion, the empty droplet should acquire an electron from neighboring one-electron or two-electron droplets during the characteristic time $`\tau _0`$, which can be easily calculated based on the following considerations. The probability $`P`$ per unit time for an empty droplet to acquire one electron can be written as

$$P=4\pi \omega _0_0^{\mathrm{}}e^{r/l}(n_1+n_2\mathrm{exp}(A/k_BT))r^2𝑑r,$$
(16)
where the terms with $`n_1`$ and $`n_2`$ correspond to the electron transfer from single- and two- electron droplets, respectively. Performing integration in Eq. (16) and using Eq. (3), we find

$$\tau _0=\frac{\mathrm{exp}(A/2k_BT)}{8\pi \omega _0l^3n}$$
(17)
Just the same estimate can be obtained for the characteristic time of electrons leaving two-electron droplets. For our picture with empty and two-electron droplets to be valid, the following condition must be met, $`\tau _0\omega _0^1`$. Thus, our approach is valid at sufficiently low temperatures, $`k_BTA`$, and for not too small droplet density $`n`$.

The applicability of our approach also implies that $`l>a,n^{1/3}`$. It is of interest to consider also the case of $`la`$ and/or low droplet concentrations. In this situation, in usual hopping systems, the conductivity strongly depends on the geometry of current paths . This causes exponential dependence of conductivity on the carrier concentration. However, our system turns out to be more complicated that those commonly invoked for hopping conductivity. It involves different types of hopping centers giving rise to unusual geometry of current paths. Therefore, the conventional approaches used for hopping can not be applied straightforwardly to the analysis of our model at low droplet concentration or at $`la`$. Despite these complications, we believe that the expression for the conductivity in the case $`ln^{1/3}`$ includes the percolation-related factor $`\mathrm{exp}(\beta /n^{1/3}a)`$, with $`\beta `$ of the order one , though currently we have no rigorous proof of this statement. The results below for magnetoresistance and noise are insensitive to this factor, and therefore we expect them to be valid in a general case.

## III Magnetoresistance

As we already discussed, below the percolation threshold when the volume fraction of droplets $`n<n_c`$, a typical value of $`A/k_B`$ is mainly determined by Coulomb interactions between two electrons inside the droplet $`Ae^2/ϵa`$ and has a typical value of $`1000K`$. Now we can use this estimate to analyze the magnetoresistance in non-metallic phase-separated manganites. To do that, we use the expression for the radius of the magnetic polaron, obtained in the introduction $`ad(t/JS^2)^{1/5}`$. Recall once more that here $`J`$100 K is an AFM Heisenberg exchange between the local spins $`S=3/2`$. It is natural to conclude that in the magnetic field $`H`$ the Heisenberg exchange integral $`J`$ decreases according to the formula $`J(H)S^2=J(0)S^2g\mu _BHS`$, where $`\mu _B`$ and $`g`$ are the Bohr magneton and the gyromagnetic ratio, respectively. Consequently, the value of $`A`$ is decreasing linearly in the experimentally accessible range of magnetic fields, and for the excitation energy we obtain

$$A(H)=A(0)[1bH],b=\frac{1}{5}\frac{g\mu _B}{J(0)S}.$$
(18)
It follows now from Eq. (15) that the magnetoresistance is negative and for temperatures $`T<A/k_B`$ reads

$`|MR|{\displaystyle \frac{\rho (0)\rho (H)}{\rho (H)}}`$ $`=`$ $`\mathrm{exp}\left[{\displaystyle \frac{A(0)A(H)}{2k_BT}}\right]1`$ (19)
$`=`$ $`\mathrm{exp}\left({\displaystyle \frac{bHA}{2k_BT}}\right)1.`$ (20)
For low magnetic fields and not very small temperatures the absolute value of the magnetoresistance is small, $`|MR|=bHA/2k_BT1`$. In higher fields (but still $`bH1`$) the absolute value of magnetoresistance eventually exceeds one and behaves in the exponential fashion, $`|MR|=\mathrm{exp}(bHA/2k_BT)`$. Note that for temperatures $`TA/k_B`$ and for typical gyromagnetic ratios $`g10`$ the magnetoresistance in our region of doping becomes larger than one by absolute value only in relatively high magnetic fields $`H10`$T.

## IV $`1/f`$ noise power

Recently, Podzorov et al  reported the observation of giant $`1/f`$ noise in perovskite manganites in the phase separated regime. Generally, systems with distributed hopping lengths are standard objects which exhibit $`1/f`$ noise (for review, see Refs. ). The purpose of this Section is to study low-frequency noise within the framework of the model used to calculate the conductivity in Section II, and show that it has, indeed, $`1/f`$ form.

Starting from the Ohm’s law $`U=IL/\sigma S`$ (where $`L`$ and $`S`$ are the sample length and the cross-section, respectively) and assuming that the measuring circuit is stabilized ($`I=const`$), we can present the voltage noise at the frequency $`\omega `$, $`\delta U^2_\omega `$, in the following way,

$$\delta U^2_\omega =U_{dc}^2\frac{\delta \sigma ^2_\omega }{\sigma ^2},$$
(21)
where $`U_{dc}`$ is the time-averaged voltage, and $`\delta \sigma ^2_\omega `$ is the noise spectrum of the fluctuations of the conductivity.

If we disregard possible fluctuations of temperature in the system, the only source of the fluctuations in our model is those of the occupation numbers $`n_1`$ and $`n_2`$. Using the conservation law $`n_1+2n_2=n`$, we find from Eq. (14)

$$\delta \sigma =\sigma \frac{\delta n_2}{\overline{n}_2}\left[12\mathrm{exp}(A/2k_BT)\right].$$
(22)
We thus need to find the fluctuation spectrum $`\delta n_2^2_\omega `$. Following the general prescription , we recollect that the two-electron droplets decay via the process (ii), and the relaxation equation has the form

$$\delta \dot{n}_2=\frac{\delta n_2}{\tau (r)},\tau (r)=\omega _0^1\mathrm{exp}(r/lA/2k_BT),$$
(23)
where we have neglected the effect of the electric field. The fluctuation spectrum then reads

$`\delta n_2^2_\omega `$ $`=`$ $`\delta n_2^2_T{\displaystyle \underset{i}{}}{\displaystyle \frac{2\tau (r^i)}{1+\omega ^2\tau ^2(r^i)}}`$ (24)
$`=`$ $`8\pi \overline{n}_2\delta n_2^2_T{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\tau (r)}{1+\omega ^2\tau ^2(r)}}r^2𝑑r,`$ (25)
where $`\delta n_2^2_T`$ is the thermal average of the variation of $`n_2`$, and the summation is performed over the pairs ”empty droplet – two-electron droplet”, with $`r_i`$ being the distance between the sites in a pair. The average in Eq. (24) is essentially a spatial integral, with the main contribution coming from short distances,

$$\delta n_2^2_\omega =8\pi \overline{n}_2\delta n_2^2_T_0^{\mathrm{}}\frac{\tau (r)}{1+\omega ^2\tau ^2(r)}r^2𝑑r.$$
(26)
Note that Eq.(26) is valid for an arbitrary relation between $`l`$ and $`a`$, not necessarily for $`al`$.

We are interested below in the frequency range

$$\stackrel{~}{\omega }_0\mathrm{exp}(L_s/l)\omega \stackrel{~}{\omega }_0,\stackrel{~}{\omega }_0\omega _0\mathrm{exp}(A/2k_BT),$$
(27)
where $`L_s`$ is the smallest of the sample sizes. In this case, with the logarithmic accuracy we obtain for $`Ak_BT`$,

$`\delta U^2_\omega `$ $`=`$ $`U_{dc}^2{\displaystyle \frac{\delta n_2^2_T}{\overline{n}_2}}{\displaystyle \frac{4\pi ^2l^3}{\omega }}\mathrm{ln}^2\left({\displaystyle \frac{\stackrel{~}{\omega }_0}{\omega }}\right).`$ (28)
Thus, in the wide range of sufficiently low frequencies (27) the noise power spectrum for our system has almost $`1/f`$ form.

The variation $`\delta n_2^2_T=V_S^2(\overline{N_2^2}\overline{N}_2^2)`$ is easily found in the same way as Eq. (3),

$$\delta n_2^2_T=\frac{\overline{n}_2}{2V_s}.$$
(29)
Combining this with Eq. (28), we write the final expression for the spectral density of noise for $`Ak_BT`$ in the form

$`\delta U^2_\omega `$ $`=`$ $`U_{dc}^2{\displaystyle \frac{2\pi ^2l^3}{V_s\omega }}\mathrm{ln}^2({\displaystyle \frac{\omega _0e^{A/2k_BT}}{\omega }}).`$ (30)
## V Discussion

For the further discussion, it is convenient to rewrite Eq. (30) in the form

$$\alpha =\frac{\delta U^2_\omega V_s\omega }{U_{dc}^2}=2\pi ^2l^3\mathrm{ln}^2(\frac{\stackrel{~}{\omega }_0}{\omega }).$$
(31)
It is remarkable that the noise spectrum in our model has $`1/f`$ form up to very low frequencies. This is due to fluctuations in occupation numbers of droplets, associated with creation and annihilation of extra electron – hole pairs. This mechanism of $`1/f`$ noise is specific for our model and is not present in standard hopping conduction.

Let us estimate the numerical value of the parameter $`\alpha `$, which is the standard measure of the strength of $`1/f`$ noise. This parameter is proportional to the third power of $`l`$. Simple estimates (analogous to that presented in the Introduction for $`a`$) reveal that, in general, $`l`$ of the order or higher than $`a`$. Assuming again that the excitation energy is of the order of the Coulomb energy $`Ae^2/aϵ`$, taking $`\omega _0`$ to be of the order of the Fermi energy inside droplets (which means $`\mathrm{}\omega _0300K`$ for $`n<n_c`$), and estimating the tunneling length $`l`$ as being $`l2a`$20 Å, we arrive to a conclusion that the parameter $`\alpha `$ is of the order $`\alpha 10^{17}`$$`10^{16}`$cm<sup>3</sup> for $`T<A/k_B`$ and $`\omega `$1 Hz – 1 MHz. This value of $`\alpha `$ is by several orders of magnitude higher than in usual semiconducting materials (see Ref. ). Such a large magnitude of the noise can be attributed to the relatively low height of the potential barrier $`A`$ and to the relatively large tunneling length $`l`$. Formally, it is also related to the large value of the logarithm squared in Eq. (31).

According to Eqs. (30) and (31), the noise power and the noise parameter $`\alpha `$ are independent of the volume fraction occupied by the droplets. This result is valid in the intermediate range of $`n`$, when the droplet density is not too high and not too low. First, we assumed that the droplets are isolated point objects and the tunneling between the two droplets is not affected by a third polaron. This is only valid provided the droplet density is far from the percolation threshold, $`nn_c`$. On the other hand, the droplet density must not be too low since the conditions $`N,N_1,N_21`$ are assumed to be met. Moreover, we neglected the possibilities of the disappearance of a droplet without an electron, the formation of a new droplet due to the electron tunneling, and the decay of two-electron droplets. Thus, the characteristic times of these processes should be longer than the characteristic tunneling time, and the average tunneling distance can not be too high (see Eq. (17) and the discussion below it).

The above speculations imply that the following set of inequalities should be met, $`an^{1/3}l`$, for formula (15) for the conductivity to be valid. In general, the tunneling length should not be much larger than the droplet radius since just the same physical parameters determine these two characteristic distances. So, these inequalities could not be valid for real physical systems and it is of interest to consider the situation where $`a,ln^{1/3}`$, which is beyond the scope of our model. However, some definite conclusions concerning the magnetoresistance and the noise power can be made at present.

First, the factor $`\mathrm{exp}(A/2k_BT)`$ in the temperature dependence of the conductivity is related to the number of carriers and appears due to strong Coulomb repulsion of electrons in the droplet. It seems rather obvious that such a factor appears in the formula for the conductivity below the percolation threshold for an arbitrary relation between $`a`$ and $`l`$. On the other hand, in contrast to common hopping systems, strong $`1/f`$ noise in our model results from fluctuations of state occupation numbers. Actually, our result for parameter $`\alpha `$ (31) only relies on the fact that $`\delta \sigma /\sigma \delta n_2/\overline{n}_2`$, and thus the noise power is expected to have $`1/f`$ form below the percolation threshold for any relation between $`a`$ and $`l`$. Then, we can also conclude that in our model strong $`1/f`$ noise will be found under the (experimentally relevant) conditions $`al`$ .

Another important point is that we disregard direct Coulomb interaction between the droplets in comparison with the energy $`A`$. This is justified if the gas of the droplets is diluted, $`n^{1/3}a`$. In this respect, we recollect that in standard hopping conduction systems (doped semiconductors) the main mechanism of low-frequency noise is an exchange of electrons between the infinite cluster and nearby finite clusters. In the absence of interactions it leads to the noise power proportional to $`\omega ^\alpha `$, with the exponent $`\alpha `$ being considerably below one . To explain $`1/f`$ noise in these systems, models involving Coulomb interactions have been proposed . These sources of low-frequency noise are thus beyond our discussion. We also did not consider sources of noise different from resistance fluctuations. At least two other types of noise are inevitably present in the system, Nyquist-Johnson (thermal) noise, which is a consequence of the fluctuation-dissipation theorem, and shot noise due to the discrete nature of electron charge (see Refs.  for review). Both these noises are frequency independent (white) at low frequencies. The magnitudes of Nyquist-Johnson, shot, and $`1/f`$ noises are governed by absolutely different parameters, and we do not attempt to compare them here, noting only that at low frequencies $`1/f`$ noise must dominate.

In our model, we assumed that the number of droplets $`N`$ is fixed and strictly equal to the number of extra electrons. In actual systems, $`N`$ can also fluctuate, and this can be an additional source of noise, and of $`1/f`$ noise, in particular. However, this contribution depends critically on the heights of corresponding energy barriers and can vary for different systems.

As we have already mentioned, the main motivation of our work was the experimental study , which observed high $`1/f`$ noise power at high temperatures far from the metal-insulator transition. In the same experiment, noise dropped to much lower levels at low temperatures in the metallic phase. This behavior of the noise power is consistent with the present model since in the metallic phase the electron tunneling contribution to the total conductivity is negligible. In the vicinity of the percolation transition the noise power increases drastically . In this paper we do not attempt to describe the system of magnetic polarons close to the percolation threshold. However, we argue that the amplitude of $`1/f`$ noise is already large in the phase-separated regime even far from the percolation threshold.

## VI Conclusions

We emphasize that even in our oversimplified model we get a reasonable behavior of resistivity and magnetoresistance for underdoped manganites. Moreover, we have shown that in the framework of our model $`1/f`$ noise appears in the natural way. The phase-separation ensures a large magnitude of the noise power as compared with homogeneous materials.

Of course, a more sophisticated theory should include both the ferromagnetic structure of the droplet and the antiferromagnetic structure of the insulating matrix. This can lead us to the physics resembling that observed in the process of spin-assisted tunneling, which attracts a considerable interest nowadays (see e.g. Ref. ). The work in this direction is in progress.

## Acknowledgments

This work was supported by Swiss National Science Foundation, INTAS (grants 97-0963 and 97-11954), Russian Foundation for Basic Research (grants 00-02-16255 and 00-15-96570), and the Russian President Program (grant 96-15-9694). We acknowledge stimulating discussions with P. Fulde, D. Vollhardt, K. Held, P. Wölfle, A. F. Andreev, Yu. Kagan, A. S. Ioselevich, M. S. Mar’enko and D. V. Efremov. M. Yu. K. would like also to acknowledge hospitality of the Max-Planck-Institut in Dresden where a part of this work was done.


text:
# Dynamics of Lennard-Jones clusters: A characterization of the activation-relaxation technique

## I Introduction

In many materials and systems, microscopic structural relaxation takes place on time scales much longer than those of the atomistic oscillations set by the phonon vibrations (10 ps). This is the case, for example, for glasses and other complex materials. Such a time spread in the dynamics can be understood from the configurational energy landscape picture. Indeed, the system finds itself in a deep minimum surrounded by energy barriers much higher than its thermal energy and only rare fluctuations can allow the system to jump over a barrier and move to a new minimum. These long time scales are especially prohibitive for numerical studies using traditional methods such as molecular dynamics (MD) and real space Monte Carlo, which are tied to the phonon time scale.

One way to reach this long time scale is through activated dynamics . In this case, the algorithm focuses directly on the appropriate mechanisms and describes the dynamics as a sequence of metastable states separated by energy barriers. These metastable configurations can be well identified by their atomic positions at zero K, which correspond to a local minimum in the configurational energy landscape. Knowledge of the distribution and properties of these local minima is sufficient to determine the thermodynamical properties of the system. A proper description of the dynamics, however, also requires knowledge of the rates controlling jumps from one local minimum to another.

Activated dynamics has been successfully applied to a range of discrete problems. It has turned out to be especially useful in the study of metallic surfaces, where it is possible to identify and compute a priori the whole set of barrier to be visited during the dynamics . These methods have also be applied, with further approximations, to other systems such as the hetero-epitaxial growth of semiconductor compounds . Such simulations can reach simulated times 10 or 12 orders of magnitude longer than what can be done with molecular dynamics.

For more complex systems, however, identifying the barriers and insuring proper statistical sampling remain a challenge. The activation-relaxation technique (ART)  was proposed recently to address this challenge. A number of questions remain regarding the biases of the method, how it samples the potential energy landscape and whether or not it is reversible.

Before achieving the long term objective of developing an algorithm for simulating the atomistic time evolution of complex materials, it is necessary to address these issues. Here, we answer some of these questions in a study of Lennard-Jones (LJ) clusters, comparing ART of Barkema and Mousseau with a similar algorithm introduced by Doye and Wales(DW) . In particular, we look at the sampling of barriers, the reversibility of the paths from the saddle-point and from the new minimum. The dynamical and thermodynamical properties of LJ clusters have been thoroughly studied  and they provide an ideal model for the development of global optimization techniques  as well as activated methods.

This Paper is constructed as follows. In the next section, we describe the details of the activation-relaxation technique. We then present the results of our simulation on 3 different LJ clusters (13, 38 and 55 atoms) and give a short discussion.

## II Technique

The activation-relaxation technique (ART),  is a generic method to explore the surface energy landscape and search for saddle-points. It has been applied with success to amorphous semiconductors, silica and metallic glasses .

ART defines events in the configurational energy landscape according to a two-step process: (a) starting from a local minimum, a configuration is pushed up to a local saddle point, representing the activation process; (b) from this saddle point, the configuration is relaxed into a new minimum; this whole process is called an event. For the activation, we use a modified version of the algorithm, ART nouveau , which was introduced recently . It is now possible to follow the direction corresponding to the negative eigenvalue exactly, ensuring a full convergence onto the saddle point.

The configuration is first pushed along a random direction until a negative eigenvalue appears. At each step along this trajectory, its total energy is relaxed in the hyperplane perpendicular to it. This ensure that the total energy and forces remain under control as the configuration leaves the harmonic well.

Once the lowest eigenvalue passes a threshold (set here at 10<sup>-3</sup>), we start the convergence to the saddle point by pushing the configuration along the eigenvector corresponding to this lowest eigenvalue, while minimizing the forces in all other directions. Unless the lowest eigenvalue turns positive, this procedure is guaranteed to converge onto a first-order saddle point, where forces on the configuration are zero. If the lowest eigenvalue changes sign, the iteration procedure is stopped and a new event is started.

Because ART is designed to work for systems with thousands of degrees of freedom, it is not appropriate to perform a direct diagonalization of the Hessian to extract eigenvalues and eigenvectors. We use instead the Lanczos algorithm . This algorithm works by iteratively projecting a vector on the Hamiltonian, extracting preferably the lowest eigenvalues and corresponding eigenvectors. In our case, 15 to 30 force evaluations are sufficient to extract the very lowest eigenvalue and its eigenvector, requiring the diagonalization of a trigonal matrix of the same dimensions. As a bonus, this iterative scheme can use the direction of the previous step as a seed, ensuring a better convergence.

The second step of the algorithm, the relaxation, is straightforward and can be achieved with any standard minimization technique. We use the conjugate gradient method (CG) .

Since moves are defined directly in the 3N-dimensional configurational space, ART is not sensitive to the constraints of real space algorithms: a complex collective motion, requiring the displacement of hundreds of atoms, is as easily produced as a one-atom jump; and a high energy barrier does not require more efforts to cross than a thermal one. Such versatility is particularly important in disordered and complex materials where events can involve collective rearrangements that are hard to foresee.

In this work, the energy landscape is described by the Lennard-Jones potential:

$$E=4ϵ\underset{i<j}{}\left[\left(\frac{\sigma }{r_{ij}}\right)^{12}\left(\frac{\sigma }{r_{ij}}\right)^6\right],$$
(1)
where $`ϵ`$ is the pair well depth and 2$`{}_{}{}^{1/6}\sigma `$ is the equilibrium pair separation. The energy and distance are described below in units of $`ϵ`$ and $`\sigma `$, respectively.

As mentioned above, we compare our results to those obtained by the Doye and Wales (DW) version of ART. DW propose a systematic technique for exploring the surface energy landscape. Transition states are found by using the eigenvector-following method , in which the energy is maximized along one direction and simultaneously minimized in all the others. In this approach, the Hessian matrix is diagonalized at the local minimum and each of its eigenvectors are followed in turn in both directions away from this minimum. Although there is no information regarding the position of saddle points in the Hessian at the local minimum, the eigenvalue of the eigenvectors followed from this point often moves down and, a some point, might becomes negative. From then on, the procedure is similar to that described above.

The main advantage of this algorithm is that it provides a systematic way of exploring the local energy landscape, moreover, we can expect that its biases will be different from those of ART nouveau described above. It suffers some limitations, however. First, the number of trial direction is finite, leading to a maximum of 6N saddle points. As we will see below, even in small clusters, this is not enough to sample all saddle points. Second, the method requires the full diagonalization of the Hessian matrix repeatedly, at least in its first stage, making it an order $`N^3`$ technique, too costly for problems of more than a few hundred degrees of freedom.

## III Simulations

Simulations are done for the 13-atoms, 38-atoms and 55-atoms Lennard-Jones clusters, using both ART and DW. In all cases, we start from a relatively well-relaxed generic configuration, i.e., one which does not have special symmetries, and explore the energy landscape around this minimum. The goal here is not to recreate the full connectivity tree –this was done already by Doye et al — but to study the biases of the methods in finding events.

### A A comparison between ART and DW

For DW, the number of search directions is limited to $`6N`$. A number of these directions do not converge to a saddle point or lead to degenerate activated points, producing of order $`N`$ structurally different saddles. Contrary to DW, ART can generate an infinite number of initial search directions. In the first part of the simulation, we limit ourselves to sets of 3000 trial events, starting from the same initial minimum, for each cluster size.

After eliminating degenerate saddle configurations, ART is found to have generated three to four times more events than DW for these clusters. This ratio is, of course, related to the number of trial directions used in the ART simulation. As discussed below, for example, a 20 000 trial run on the 13-atom cluster can find all first-order paths leaving a given local minimum.

Results for these runs, including statistics on the saddle points and new configurations are reported in Table I. In all cases, permutational isomers are eliminated and only structurally different configurations are counted.

The sets of saddle points and minima obtained by ART and DW are obviously not independent and we find, as one would expect, a significant overlap between them: about 80 % of DW events are also found by the 3000-attempt run of ART (while about 20 % of ART events are also found by DW).

We can do a similar analysis for the new minima found by the two methods. Table I shows that in many cases, a number of different saddle points lead to the same final minimum. Moreover, this degeneracy seems to increase with the size of the cluster. (The existence of multiple paths connecting two minima is much more common for open systems, such as clusters, than for bulk materials. In the later case, the constraints of volume and continuity make it much more difficult to find many paths connecting two events.) For all cluster sizes, the 3000-trial event ART simulations finds about 3 times more minima than the DW method, a ratio similar to the saddle configurations.

### B Ergodicity and reversibility

It is possible to examine the question of the ergodicity of ART by extending the simulations described above. In Figure 1, we trace the number of different saddle points and minima as a function of trial events in a 20 000-attempt run, using a single initial minimum as in the previous simulation. The sampling of minima seems to be complete after about 9000 events and we identify a total of 79 different minima, including all those found using the DW method. It takes about 50 % more events to generate all 195 saddle points that can be found with ART. Here again, all saddle points generated with DW are found by ART. Although this does not show formally the completeness of the sets found, the previous simulation, as well as its comparison with a different technique, provides a fairly solid base for claiming that ART is ergodic.

Dynamically, events should also be fully reversible. We check this in a two-step simulation. First, we check the reversibility of the paths from the saddle point. After converging to a first order saddle point, we pull the configuration back very slightly and apply the relaxation routine from there. For the 3000 ART-generated events on the 13-atom LJ cluster, all but 10 (0.33%) events relax to the original minimum and these last 10 events converge to a different but very close local minimum. Saddle points found by ART are therefore really delineating the boundary of the energy basin around the initial minimum.

Clearly, paths are fully reversible from their activated point.

The full trajectory must also be reversible; the initial minimum, saddle, final minimum sequence should also be found in the opposite order. For each LJ cluster, we select 20 new local minima, reached in a one-step event from the initial minimum and check that, indeed, the configuration can alway come back to that state, and that the same saddle points are found in both directions. This important result allows us to conclude that the whole pathway between local minima passing through first order saddle-points is reversible.

### C Stability under change of parameters

In the search for saddle points, a number of factors can influence the selection of the activation paths. In particular, it is important to verify that the step size along the direction of the lowest eigenvalue, in the convergence to the saddle point, does not result in some events missing. To examine this issue, we perform the same simulation as that described above for the 13-atom LJ cluster with two different step values: 0.01 and 0.03 $`\sigma `$. The distribution of energy barriers for all converged saddle points is given in Fig. 2. All saddle points found by the 0.01 $`\sigma `$ step are also found by the larger step. Moreover, the 0.03 $`\sigma `$ seems to better converge to higher energy barriers. The total number of saddle-points reached by using 0.03 $`\sigma `$ step value is also four times greater than that found by 0.01 $`\sigma `$. This sheds some light on the workings of the algorithm. As we first start to follow the eigendirection corresponding to the lowest eigenvalue, the configuration is in a very shallow valley and too much relaxation perpendicular to this valley can easily make it vanish. With a larger step moving away from the minimum, the configuration reaches a deep valley faster, increasing strongly the rate of convergence. It is clear, here, that this parameter can be adjusted to optimize the rate of convergence without fear of loosing particular saddle points.

### D Biases in searching for saddle points

We now address the question of biases in the search for saddle points and new minima. The range of activation and asymmetry energies is quite wide for the events generated on the three clusters, as can be seen in Table I. As expected, this range increases with the size of the system. This behaviour is typical of clusters. For bulk system, the distribution of relevant activation events is generally independent of the size and is usually bounded by a small multiple of the binding energy between two atoms. The higher bound for the distribution of activation energy reflects the overall properties of the systems studied and should be rather independent of the starting configuration; the lower bound, however, is not; it represents a direct measure of the stability of the initial metastable state.

Figure 3 gives two representations of the distribution of activation and asymmetry energies for the 55-atom clusters. In Fig. 3(a), we give the distribution of energies for the 151 different saddles and 89 different minima identified in the 3000-trial event run. The distribution of saddle energies is rather smooth over its whole range while that of the asymmetry energies presents a few peaks.

To study the biases of ART, we also plot in Fig. 3(b) and (c) the normalized distribution of activation and asymmetry energies for all 1856 completed ART events and the corresponding 109 successful DW events (although there are a total of 330 trial directions for DW, the method shows a poor success rate when eigendirections corresponding to very high eigenvalues are selected.) The distribution for DW is therefore almost discretized and peaks correspond to single saddle points being visited many times. The ART distribution is almost continuous and we have binned the events as a function of energy.

The structure of these distributions is quite different from that of Fig. 3(a). In particular, there are strongly peaked at low energies indicating that each method seems to enhance significantly a few paths over the other ones. The biases are not exactly the same, however, although there is considerable overlap for the two methods. This suggests that the topography of the energy landscape around the minimum is reflected, at least in part, in the choice of events.

One approach to identify the overall biases of ART is to plot the ratio of the distribution of barriers for all saddle points generated (frequency of convergence) over that of the different saddle points existing around our minimum. This ratio, plotted in Fig. 4, provides a first indication on how ART selects saddle points. In Fig. 4(a), the distribution is dominated by a strong peak at low activation energy, followied by a fairly extended tail. Surprisingly, most of the small structure that can be found in 3(a) has been eliminated in the ratio, indicating that there is fairly little bias towards a few events in particular but that the selection of an event over another one is mostly a matter of energy. This can be seen more clearly in the log-normal distribution plotted in Fig. 4 (b). The result is remarkable: the distribution is well fitted by an exponential function $`p(E)\mathrm{exp}(0.57E)`$ and the sampling follows a Boltzmann distribution! The same behavior can be observed for the 13-atom and the 38-atom clusters, as shown in Fig. 5.

## IV CONCLUSION

The configurational energy landscape of Lennard-Jones clusters was explored by using the activation relaxation technique (ART). Comparing ART with a different sampling method proposed by Doye and Wales, it does not seem that ART missed any class of first-order activated paths. Based on this and extensive simulations, we conclude that ART can find all first-order saddle points and new minima around a given minimum, indicating that it is ergodic. The trajectories it defines, initial minimum–saddle point–final minimum, are also reversible, indicating that the trajectories found are real activated paths.

We find, finally, that ART samples the surface energy landscape according to a Boltzmann distribution, i.e., that the probability of finding a given saddle point is proportional to $`\mathrm{exp}(E_{barrier})/E_0`$. The original of this bias is not yet understood but this opens the door to real-time activated dynamics in complex system, a very exciting prospect.

## V Acknowledgment

We acknowledge useful discussions with G. Barkema as well as partial support from NSF under grant number DMR-9805848.


text:
# Off equilibrium properties of vortex creep in superconductors

## Abstract

We study a model for the dynamics of vortices in type II superconductors. In particular, we discuss glassy “off equilibrium” properties and “aging” in magnetic creep. At low temperatures a crossover point is found, $`T_g`$, where relaxation times seem to diverge á la Vogel-Tamman-Fulcher. Magnetic creep changes by crossing $`T_g`$: above $`T_g`$ power law creep is found asymptotically followed by stretched exponential saturation; below $`T_g`$ the creep is logarithmic and vortex motion strongly subdiffusive. In this region violation of time translation invariance is found along with important dynamical scaling properties.

Despite apparent structural differences, the dynamics of vortices in superconductors  has many common features with glass formers such as supercooled liquids, spin glasses or polymers . In superconductors even simple quantities such as the sample magnetisation, $`M`$, exhibit strong “glassy” behaviours and a variety of “history dependent” effects has been discussed . Common examples are hysteretic magnetisation loops which at low temperature exhibit “memory” effects: the value of $`M`$ for a given applied field, $`H`$, depends on the sweep rate of $`H`$, i.e., on the history of the sample . “Aging” has been also recently observed in magnetic relaxation in type II superconductors . Interestingly, these phenomena are very similar to those found in other glassy systems such as supercooled liquids. For instance, after a quench, “cooling rate” dependences and “aging” effects are usually recorded in particle density measurements .

Those similarities are striking when one considers the differences which exist at a microscopic level between these systems. One is apparent: the interaction ranges of vortices in superconductors are typically much larger than those of molecules in glass formers. In this perspective our aim is to discuss the above connections and relate off equilibrium dynamics of vortex matter to the dynamics of glass former liquids, where relevant progresses have been recently accomplished . In all these systems, the glassy features are due to some general mechanisms and unrelated to specific material parameters. Thus, a statistical mechanics approach, disregarding sample details, seems appropriate to describe them. For instance, many of the features of glassy behaviours in fragile glasses, such as colloidal suspensions, are well described in the context of spin glass models .

Extensive work has been done on the equilibrium or close to equilibrium properties of the so-called vortex glass phase proposed by M.P.A. Fisher (see ref.s in ). Here we discuss, in particular, the dynamics of the system when far from equilibrium, a typical situation due to the enormous relaxation times found in the low temperatures region. We consider a simple tractable model for vortices  which can be well understood in the context of standard statistical mechanics of disordered systems . Such a model was already shown to successfully depict a unified picture of creep and transport phenomena in vortex physics, ranging from magnetisation loops with “anomalous” second peak, logarithmic relaxation at vanishing $`T`$, Bean profiles, to history dependent behaviours in vortex flow and I-V characteristics, to the reentrant nature of the equilibrium phase diagram .

As mentioned, here we study in details the system off-equilibrium creep dynamics. The model gives the following scenario for vortex dynamics. At low temperatures a crossover point is found, $`T_g`$, where relaxation times become longer than the typical observation scale . Actually, they seem to diverge at a lower temperature, $`T_c`$, where an “ideal” thermodynamic glass transition can be located (in 2d $`T_c`$ is numerically zero). Magnetic creep changes when $`T_g`$ is crossed: above $`T_g`$ power law creep is found asymptotically followed by stretched exponentials; below $`T_g`$ the creep is logarithmic and shows strong form of “aging”, such as violation of time translation invariance with apparent dynamical scaling properties. Microscopically, $`T_g`$ corresponds to a drastic change in vortex motion which becomes strongly subdiffusive. We also outlines the connections of vortex off equilibrium dynamics with other glassy systems.

We study a Restricted Occupancy Model (ROM) , a statistical mechanics model which in the limit of zero temperature and infinite upper critical field reduces to a cellular automaton introduced in  to study vortex avalanches. The model is a coarse grained version  of a system of straight vortex lines: a set of repulsive particles diffusing in a disordered pinning environment. The coarse graining length, $`l_0`$, is taken to be of the order of the natural screening length of the problem (which in our case is the London penetration length $`\lambda `$). Thus, the coarse graining technique  reduces the original interaction potential to an effective short ranged one. The price to pay is a loss of information on scales smaller than $`l_0`$ and the necessity to allow multiple occupancy on the lattice sites of the coarse grained system. The ROM model is defined by the Hamiltonian : $`=\frac{1}{2}_{ij}n_iA_{ij}n_j\frac{1}{2}_iA_{ii}n_i_iA_i^pn_i`$, together with the constraint $`0n_iN_{c2}`$ on the integer occupancy variable $`n_i`$ representing the number of vortices of the coarse grained cell corresponding to site $`i`$ (bounded by the upper critical field).

The first term in $``$ represents the repulsion between particles . For the above considerations on vortex interactions, in the present coarse grained representation (where $`l_0\lambda `$) a finite range potential must be considered. For simplicity we choose $`A_{ii}=A_0`$; $`A_{ij}=A_1`$ if $`i`$ and $`j`$ are nearest neighbours; $`A_{ij}=0`$ for all others couples of sites. The second term in $``$ concerns the particle self-interaction energy and the third one is a pinning potential, with a given distribution $`P(A^p)`$, acting on a fraction $`p`$ of lattice sites (below $`p=1/2`$). For simplicity we choose: $`P(A^p)=(1p)\delta (A^p)+p\delta (A^pA_0^p)`$. $`A_0`$ sets the energy scale and we choose $`A_0=1.0`$; $`A_0^p=0.3`$; $`N_{c2}=27`$; $`\kappa ^{}A_1/A_0[0,0.3]`$ (below we typically discuss data for $`\kappa ^{}=0.28`$). The relation of the parameters of the model to material parameters is shown elsewhere , here we only recall that $`\kappa ^{}`$ is an increasing function of the Ginzburg-Landau $`\kappa `$ and $`N_{c2}H_{c2}`$.

Through its surfaces, the system is in contact with an external reservoir of “particles”, which schematically corresponds to the applied field present in magnetic experiments on superconductors. Particles are introduced and escape the system only through the reservoir, which, by definition, has a given density $`N_{ext}`$. We performed Monte Carlo Kawasaki (MC) simulations  on a square lattices of linear size $`L`$ ($`L\{8,\mathrm{},128\}`$) in presence of a thermal bath of temperature $`T`$. Our numerical statistical averages run, according to system size, from $`128`$ to $`512`$ pinning realizations.

The magnetisation in the system at time $`t`$  is defined as: $`M(t)=N_{in}(t)N_{ext}(t)`$, where $`N_{in}=_in_i/L^2`$ is the overall density inside the system and $`N_{ext}`$ is the applied field. As in usual zero field cooled experiments, we record the isothermal relaxation of $`M(t)`$ after ramping at a given “sweep rate”, $`\gamma `$, the field from zero up to a given working value $`N_{ext}`$.

The presence, at low temperature, of sweep rate dependent hysteretic cycles, slowly relaxing magnetisation, and similar effects , indicate that our system, on the observed time scales, can be far from equilibrium. Off-equilibrium behaviour is appropriately described by two times correlation functions. Thus, along with $`M(t)`$, we also record the magnetic correlation function ($`t>t_w`$):

$$C(t,t_w)=[M(t)M(t_w)]^2.$$
(1)
Since the general dynamical features of $`C`$ and $`M`$ are very similar, for simplicity we mainly discuss $`C`$ which gives richer information.

At not too low temperatures, for instance at $`T=1.0`$, the system creep is characterised by finite relaxation times, but the dynamics is already highly non trivial. The two times correlator $`C(t,t_w)`$ is plotted in Fig.1 after a ramp with $`\gamma =10^3`$ at $`T=1.0`$. At long times, $`C(t,t_w)`$ is well fitted by the so called Kohlrausch-Williams-Watts (KWW) law (see Fig.1):

$$C(t,t_w)C_{\mathrm{}}\left\{1e^{[(tt_w)/\tau ]^\beta }\right\}$$
(2)
The KWW decay is observed in the asymptotic relaxation of superconductors as well as in glass formers (their so called $`\alpha `$-relaxation ) above the glass transition . The time scale $`\tau `$ and the Kohlrausch-exponent $`\beta `$ depend on the temperature $`T`$ (see Fig.2) and on the overall field $`N_{ext}`$. The pre-asymptotic dynamics (i.e., $`t,t_w<<\tau `$) is also interesting and characterised by various regimes. In particular, for not too short times, a power law is observed over several decades (see Fig.1):

$$C(t,t_w)C_0\left(\frac{tt_w}{\tau }\right)^a$$
(3)
The exponent $`a`$ is almost independent of $`N_{ext}`$, $`a1.7`$, except at very small or high fields. Notice that the $`\tau `$ in eq.(3) is the same as in eq.(2), but the exponents $`a`$ and $`\beta `$ are numerically different (see Fig.2).

Interestingly, the above finding in the model of power laws followed by KWW relaxations is typically observed, for not too low $`T`$, in superconductors . This behaviour is also typical of supercooled liquids, where the power law regime is called the $`\beta `$-relaxation .

At $`T=1.0`$ or, generally, at not too low temperatures, no “aging” is seen: $`C(t,t_w)`$ is a function of $`tt_w`$. This is clearly shown in Fig.1, where we plot $`C(t,t_w)/C_{\mathrm{}}`$ for several different values of $`N_{ext}`$ as a function of the scaling variable $`(tt_w)/\tau `$: all the data (for all $`t_w`$ and $`N_{ext}`$) fall on the same master function (which is more general than the above KWW fit).

The scenario described for $`T=1.0`$ is found in a broad region at low temperatures. However, around $`T=0.5`$, $`\beta `$ drastically decreases and, at the same time, a steep increase of $`\tau `$ is found (see Fig.2). For instance, at $`N_{ext}=10`$ for temperatures below $`T_g0.25`$, the characteristic time gets larger than our recording window. Thus, below $`T_g(N_{ext})`$ the system definitely loses contact with equilibrium during our observation and, as shown in details below, typical glassy phenomena, such as “aging”, are observed. The crossover temperature $`T_g`$ is itself a function of $`\gamma `$. It has a physical meaning similar to the so called phenomenological definition of the glass transition point in supercooled liquids. Exploiting this analogy we will call this temperature the glass temperature despite the fact that in glassy systems it is loosely defined . The presence of an underlying “ideal” glass transition point, $`T_c(N_{ext})`$, is a non trivial possibility which, in many cases (as supercooled liquids) is still under debate . $`T_c`$ is often located by some fit of raw $`\tau `$ data collected in the high $`T`$ regime .

In our model, in the region where $`\tau `$ has a steep increase a Vogel-Tamman-Fulcher law (VTF) fits the data  (see inset of Fig.2):

$$\tau =\tau _0\mathrm{exp}\left(\frac{E_0}{TT_c}\right)$$
(4)
For example, at $`N_{ext}=10`$, the characteristic time $`\tau _0`$ is very large, $`\tau _0=1.610^3`$, and the characteristic activation energy, $`E_0`$, is ten times larger than $`T_c`$: $`E_0=1.1`$ and $`T_c=0.1`$, thus the above fit is consistent with an Arrhenius low (i.e., with $`T_c=0`$). The presence of a strong increase of $`\tau `$ close to a power law or a VTF law is again a mark of the apparent similitude with glassy features of supercooled liquids and glasses .

Since below $`T_g`$ relaxation times are huge, one might expect that the motion of the particles essentially stops, apart from their vibration inside cages of other vortices. Instead, as shown below, the off equilibrium dynamics has remarkably rich properties. In particular, the properties of the system slowly evolve with time. These phenomena are usually summarised by claiming that the systems is “aging”. They occur in many different systems ranging from polymers, to supercooled liquids , spin glasses  or granular media , and their origin and apparent universality are still an open problem .

In the inset of Fig.3 we show that $`C(t,t_w)`$, at $`T=0.1`$, exhibits strong “aging”. Notice that $`C`$ depends on both times $`t`$ and $`t_w`$: the system evolution is slower the older is its “age” $`t_w`$ (here, as in eq.(1), $`t_w`$ is the time elapsed from the sample preparation at the working field and temperature). Note the contrast with the case at $`T=1.0`$ of Fig.1. Such a dynamical “stiffening” is typical of glass formers . It is important to stress that the presence of slowly relaxing quantities, such as $`M(t)`$, must not lead to conclude that the system is close to equilibrium . In the entire low $`T`$ region ($`T<T_g`$), after a short initial power law behaviour, $`C(t,t_w)`$ can be well fitted by a generalisation of a known interpolation formula, often experimentally used , which now depends on the waiting time, $`t_w`$:

$$C(t,t_w)C_{\mathrm{}}\left\{1\left[1+\frac{\mu T}{U_c}\mathrm{ln}\left(\frac{t+t_0}{t_w+t_0}\right)\right]^{1/\mu }\right\}$$
(5)
We found that to take $`\mu 1`$ is consistent with our data. In the above fit $`U_c/\mu T`$ only depends on $`N_{ext}`$ and, interestingly, $`t_0`$ is approximately a linear function of $`t_w`$: $`t_0t_w+t_0^{}`$, where $`t_0^{}`$ is a constant. Very interesting is the presence of scaling properties of purely dynamical origin in the off-equilibrium relaxation. This is shown in Fig.3, where data for different fields, $`N_{ext}`$, and different waiting times, $`t_w`$, are rescaled on a single master function. The above results imply that for times large enough (but smaller than the equilibration time), $`C(t,t_w)`$ is a universal function of the ratio $`t/t_w`$: $`C(t,t_w)𝒮(t/t_w)`$. Such a behaviour (called “simple aging” ) is in agreement with a general scenario of off-equilibrium dynamics (see Ref.) and has strong analogies with other glassy systems . Experimental data on $`C(t,t_w)`$ do not exist in vortex matter, but would be extremely important to shed light on the real nature of dynamical phenomena in superconductors.

Finally, we stress that in vortex physics transitions from low temperature logarithmic to higher temperature power law creep are usually experimentally found (see references in ) and give a way to approximately locate the position of $`T_g`$ in real samples.

In the above scenario, microscopic quantities concerning the internal rearrangement of the system, such as the vortex mean square displacement, $`R^2(t)`$, can be insightful. Since vortices can enter and exit the sample, $`R^2(t)`$ must be properly defined. To this aim, we also made a different kind of computer runs where after ramping the field from zero to $`N_{ext}`$ we close the system (i.e., remove the reservoir). Since now the number of particles is fixed , their positions, $`\stackrel{}{r}_i(t)`$, are well defined at each time step and we record: $`R^2(t)=\frac{1}{N}_i(\stackrel{}{r}_i(t)\stackrel{}{r}_i(0))^2`$, were $`N`$ is the total number of present particles.

$`R^2(t)`$ is plotted in Fig.4 for $`N_{ext}=10`$. At high enough $`T`$, $`R^2(t)`$ is linear in $`t`$, but at lower temperatures it shows a pronounced bending. Finally, below $`T_g`$, the process at long times becomes subdiffusive:

$$R^2(t)t^\nu $$
(6)
with $`\nu <<1`$. From this point of view, $`T_g`$ is the location of a sort of structural arrest of the system, where particle displacement is dramatically suppressed. Interestingly, a very similar scenario has been recorded in real superconducting samples: for instance in Ref. it was clearly shown that vortices are definitely mobile in the low temperature phase and only “freeze” below a characteristic field dependent temperature.

Summarising, the schematic coarse grained model we considered shows magnetic properties very close to those experimentally found in superconductors . Here, in particular, we have focused on the off equilibrium creep dynamics and an interesting scenario emerges. At low temperatures the system relaxation time $`\tau (N_{ext},T)`$ grows enormously with an “ideal” divergence á la VTF at some $`T_c(N_{ext})`$ (interestingly, molecular dynamics simulations of London-Langevin models seem to confirm these results ). In fact, below a crossover temperature, $`T_g>T_c`$, vortex displacement becomes sub-diffusive, signalling the presence of a structural arrest. It is impossible to equilibrate the system on the observation time scale. At low temperatures creep is logarithmic and becomes a power law above the crossover $`T_g(N_{ext})`$.

In typical experiments or computer simulations, magnetic properties are measured after ramping the external field at a given rate $`\gamma `$. Whenever $`\gamma `$ is much larger than the inverse of the characteristic relaxation time, $`\tau (N_{ext},T)`$, the system is naturally driven off equilibrium, simply because it is unable to follow the drive, and “hysteresis”, “memory” effects, along with dependences on the sweep rate, occur. In the present framework, the strict correspondence between the dynamics of vortices and other glass formers can be also rationalised.

Work partially supported by EPSRC, INFM-PRA(HOP)/PCI.


text:
# The Jahn-Teller theorem for the 3d magnetic ions*

## Abstract

We argue that the Jahn-Teller theorem, if applied to 3d-ion compounds, has to be considered in the orbital+spin space. It is despite of the weakness of the intraatomic spin-orbit coupling.

Keywords: Jahn-Teller theorem, spin-orbit coupling, 3d ions

PACS: 71.70.E. 75.10.D 75.30.Gw

Receipt date by Phys.Rev.Lett.: 16.02.1999

(16 February 1999)

There is ongoing strong discussiona in Phys.Rev.Lett. and Phys.Rev.B about a relationship between the electronic structure, magnetism, the local 3d-ion surrounding and a (off-cubic) lattice distortion in 3d ionic compounds like La1-xCaxMnO3, LaCoO3. In this discussion, terms like the Jahn-Teller (J-T) effect and the J-T distortions appear that obviously are related with the J-T theorem. This long-time discussion reveals that the physical picture for 3d-ion compounds is far from a consensus. On other side, the x-ray-absorption fine structure (XAFS) experiments provide more and more accurate data on the local surrounding in the atomic scale.

The aim of this short Letter is to put attention that the degeneracy space for the consideration of the J-T theorem for the 3d ions is the orbital+spin space. It is in contrast to the current-literature that considers the J-T theorem, if applied to 3d ions, in the orbital space only.

In 1937 Jahn and Teller had pointed out1,2 that no non-linear molecule could be stable in a degenerate state. It means that one would expect a spontaneous distortion to undergo in order to remove this degeneracy. In the current-literature the J-T theorem and related with it the degeneracy, if applied to 3d ions, is considered in the orbital space. It likely repeats the text-book knowledge (Ref.2 p.193 line 4 bottom, Ref. 3 p. 212. line 12 top, Ref. 4 p.40 line 8 bottom). This way of thinking is commonly used for the description of the 3d ions and 3d-ion compounds basing on the weakness of the spin-orbit coupling for these ions.

We argue that the J-T theorem has to be considered in the orbital+spin space. It is the consequence of the spin-orbit coupling that, though weak in the 3d ions, always exists.

The demanded orbital+spin space for the consideration of the J-T theorem is consistent with the standard treatment of rare-earth-ion compounds.

In conclusion, we argue that the J-T theorem, if applied to 3d-ion compounds, has to be considered in the orbital+spin space.

\*This paper has been submitted 16.02.1999 to Phys.Rev.Lett. \[LB7427\] as a consequence of a erroneous report of the referee of Phys.Rev.Lett. on our paper  that ”the Jahn-Teller theorem only applies to orbitally degenerate levels and therefore does not apply to the system considered in the paper” \[LL6530, report A from 4.12.1997\], i.e. the d8 system in the octahedral crystal field. This report was taken by The Editor of Phys.Rev.Lett. Dr G.Wells as the base for the rejection of the a/m paper. In the paper LL6530 we have shown  that the d8 system (the Ni2+ ion) in the octahedral surrounding is the Jahn-Teller ion, despite of the orbital singlet ground state. The long discussion on this paper and others (e.g. Ref. 6, submitted 30.05.1997 as LE6925) has resulted in the undertaking by the Managing Editor (Dr G.Wells) and the Editor in Chief (Dr M.Blume, 30.03.2000) of a special discriminating policy with respect to my papers. Such the policy is the manipulation of Science by the Editors of Phys.Rev.Lett. and violates the fundamental scientific rules. Our request to publish our paper with negative referee reports has been ignored.

aIn the submitted letter of 15.02.1999 some examplanary papers discussing the J-T effect have been mentioned: Phys.Rev.Lett. 80 (1998) 853, Phys.Rev.Lett. 76 (1996) 4825, Phys.Rev.Lett.77 (1996) 5296.


text:
# On the origin of biquadratic exchange in spin 1 chains

\[

## Abstract

One dimensional spin 1 systems may have a rich phase diagram including Haldane gap and dimerized phases if the usually very small biquadratic exchange becomes significant. We show that this unlikely condition may be fulfilled in electron systems with quasi-degenerate orbitals. This mechanism may have been experimentally realized in the spin 1 chain LiVGe2O6. The implications for the exploration of the physics and quantum chemistry of spin 1 chains are discussed.

PACS: 75.30.Et Exchange and superexchange interactions - 75.10.Jm Quantized spin models - 75.50.Ee Antiferromagnetics

\]

There has been a lot of activity recently on the physics of one-dimensional (1D) spin 1 systems, especially after the prediction by Haldane that the 1D Heisenberg model has a spin gap for integer spins . That prediction has been confirmed since then by the observation of a gap in many spin 1 chains . There has also been a lot of progress in the study of the most general Hamiltonian describing an isotropic coupling between neighboring spins 1, namely

$`H`$ $`=`$ $`{\displaystyle \underset{ij}{}}H_{ij},`$ (1)
$`H_{ij}`$ $`=`$ $`J_1\stackrel{}{S}_i\stackrel{}{S}_j+J_2(\stackrel{}{S}_i\stackrel{}{S}_j)^2,`$ (2)
where the sum $`ij`$ is over nearest–neighbor pairs. The phase diagram of this model is extremely rich . It is most easily described using the parametrization $`J_1=J\mathrm{cos}\theta `$, $`J_2=J\mathrm{sin}\theta `$. As illustrated in Fig. 1, the system has two gapped phases: The Haldane phase for $`\pi /4<\theta <\pi /4`$ and a dimerized phase for $`3\pi /4<\theta <\pi /4`$. The two phases are connected by a critical point at $`\theta =\pi /4`$ for which the model is exactly solvable and the spectrum is gapless. In addition, there is a Lifshitz point and a disordered point in the Haldane phase, as well as a Valence Bond Solid point , and the properties for $`\pi /4<\theta <\pi /2`$ and $`\theta \pi /4`$ are not totally agreed upon yet.

Hopelessly, it has not been possible until very recently to explore this phase diagram experimentally due to the lack of systems with a sizable biquadratic exchange, which is conveniently measured by the ratio $`\beta =J_2/J_1=\mathrm{tan}\theta `$. It has been generally assumed that $`\beta `$ is always small so that the system is always in the Haldane gap phase. The experiments on many spin 1 chain compounds have indeed supported this view until a very recent experiment on the 1D, spin 1 vanadium oxide LiVGe2O6. In that experiment, an abrupt drop in magnetic susceptibility typical of a spin-Peierls transition has been observed at a temperature of 22 K . This property is consistent with a gapless spectrum above this temperature, and it has been argued that this behavior is most likely due to the presence of a significant biquadratic exchange interaction .

Although biquadratic or higher-order spin exchange interactions have been discussed in the past and shown to induce anomalous magnetic properties in spin $`S1`$ systems, they are usually small . It is therefore of interest to study the underlying mechanism for predominant biquadratic interaction.

In this paper we propose a new microscopic mechanism leading to a significant biquadratic interaction in spin 1 systems. We use a microscopic model to derive the interaction couplings $`J_1`$ and $`J_2`$ in Eq. (2) based on a perturbation theory, which is compared with exact numerical calculations. The key in our model is to include a third atomic orbital for an electron with slightly higher energy than the two singly occupied lowest orbital states, i.e. to consider the situation where we have a quasi-degeneracy of the orbitals on one site. The virtual electron transition via the third orbital level favors ferromagnetic spin interaction, which may compensate largely the antiferromagnetic superexchange interaction, leading to the predominance of the biquadratic interaction with $`\beta >1`$. We believe that this mechanism may be relevant to LiVGe2O6. In that case the three orbitals come from the $`t_{2g}`$ orbitals of the vanadium ions, and the degeneracy is lifted due to a small distortion of the octahedra around the vanadium atoms. More generally, we believe that the exchange Hamiltonian is likely to have higher order contributions whenever the degeneracy between the last occupied and first empty orbitals is only slightly lifted. This situation could be realized in some other transition metal oxides.

Let us consider a lattice of atoms with two outer electrons per atom. The atomic orbitals of the lowest energy level is two-fold degenerate, labelled by indices 1 and 2. There is a nearby level of orbital 3, whose energy is higher by an amount of $`\mathrm{\Delta }`$. We focus on the electron interactions at the same atom, and denote by $`U`$ ($`U^{}<U`$) the direct Coulomb repulsion between two electrons in the same (different) orbitals, and by $`J_H>0`$ the exchange interaction. The exchange interaction favors a total spin 1 state for two electrons on each atom, where orbitals 1 and 2 are both singly occupied with parallel spins, and orbital 3 is empty. This is in accordance with Hund’s rule. The effective spin-spin coupling arises when the electron has a virtual transition to the neighboring atoms. Let $`t_{lm}`$ be the electron hopping integral from orbital $`l`$ at site $`i`$ to orbital $`m`$ at site $`j`$ with the electron spin conserved. If the hopping integral is small compared with $`U^{}`$, $`J_H`$ and $`\mathrm{\Delta }`$, only virtual transitions are possible, and their net effect is to induce an intersite spin-spin coupling.

The spin Hamiltonian $`H_{ij}`$ is uniquely determined by the total spin of two neighboring sites $`i`$ and $`j`$. Let $`E_S`$ the energy of the two-site system with total spin $`S`$, with $`S`$ = 0, 1 or 2. Then $`J_1=\frac{1}{2}(E_2E_1)`$, and $`J_2=\frac{E_0}{3}+\frac{E_2}{6}\frac{E_1}{2}`$. Treating the hopping integrals as small parameters, the energy $`E_S`$ can be calculated within perturbation theory. To illustrate the essential physics, we consider the simple case where $`t_{11}=t_{12}=t_{13}=0`$, and $`t_{22}=t_{33}`$. To second order in the hopping integrals, and in units of $`t_{22}^2/U^{}`$, the spin couplings are given by

$`J_1^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(a_2+a_3){\displaystyle \frac{2}{3}}\alpha ^2(a_1a_2),`$ (3)
$`J_2^{(2)}`$ $`=`$ $`0,`$ (4)
where $`\alpha =t_{23}/t_{22}`$, and the $`a^{}s`$ are related to the Hund’s couplings: $`a_1=(1j_H)^1`$, $`a_2=(1+2j_H)^1`$, $`a_3=(1+4j_H)^1`$, with $`j_H=J_H/U^{}`$. From Eqn. (4), we see that the biquadratic interaction vanishes in the second order of perturbation. This explains the smallness of the biquadratic coupling in most systems. The first term in $`J_1^{(2)}`$ arises from hopping to an occupied orbital of the neighbouring site (see Fig. 2a) and corresponds to superexchange interaction. It is positive and favors antiferromagnetic alignment of the two spins. The second term arises from the virtual transition via the orbital state $`3`$ (see Fig. 2b). Since $`a_1>a_2`$, the contribution to $`J_1`$ is negative and it favors ferromagnetic alignment of the spins. Typical virtual transitions for the first and the second terms in $`J_1`$ are illustrated in Fig. 2. These two exchange mechanisms may compensate each other, leading to a very small net value of $`J_1`$. In this situation, it is necessary to extend the perturbation to include the contributions from the fourth order terms. Since $`\mathrm{\Delta }<<U^{},J_H`$, the most important contribution arises from processes involving a virtual state of spin 1 on both sites and with one electron in the excited state of orbital $`3`$. Such a process gives an energy correction to the second order results by a ratio of order of $`\delta =t^2/U^{}\mathrm{\Delta }`$. In Fig. 3 we show a typical example of such processes. All the other fourth order contributions are smaller by a factor of $`\mathrm{\Delta }/J_H`$ or $`\mathrm{\Delta }/U^{}`$, and will be neglected in the present consideration.

Summing over these fourth order terms, we find

$`J_1^{(4)}`$ $`=`$ $`\delta [b_2^2b_1^2+(b_3+{\displaystyle \frac{b_1}{3}})^2],`$ (5)
$`J_2^{(4)}`$ $`=`$ $`{\displaystyle \frac{\delta }{2}}[b_2^2+(b_3{\displaystyle \frac{2}{3}}b_1)^2],`$ (6)
where $`b_1=\alpha ^2a_1`$, $`b_2=\alpha (a_2+a_3)`$, and $`b_3=(a_2a_3)/2\alpha ^2a_2/3`$. The correction to $`J_1`$ is simply an energy shift, which modifies slightly the location of the transition between ferromagnetic and antiferromagnetic bilinear interaction. The biquadratic interaction first appears in the fourth order in perturbation theory. While the ratio $`\beta `$ is usually small, of the order of $`\delta `$, it can be significant if the two contributions in the second order balance out and if $`\mathrm{\Delta }`$ is relatively small. From equation (6), we also note that $`J_2<0`$, which suggests such a system is always in the lower half plane of the phase diagram in Fig. 1. This also implies $`\beta >0`$ if $`J_1>0`$, and $`\beta <0`$ if $`J_1<0`$. In Fig. 4b, we plot the value of $`\beta `$ as a function of $`\alpha `$ for given $`\mathrm{\Delta }`$ and $`J_H/U^{}`$. At $`\alpha <0.5`$, $`\beta `$ is vanishingly small. As $`\alpha `$ increases, $`\beta `$ increases sharply, and the divergent point corresponds to the transition from antiferromagnetic to ferromagnetic bilinear interaction.

As stated above, the spin couplings are related to the energy levels of the two spin 1 system as a function of the total spin. These energies can be numerically calculated exactly. For a two-site system with four electrons and three orbitals at each site, we have calculated exactly the 495 energy levels for various parameters. As long as $`\mathrm{\Delta }`$ is not too small as compared to the largest hopping integral, the low-energy part of the spectrum consists of three levels with total spin $`S=0`$, $`S=1`$ and $`S=2`$ respectively (Fig. 4a). So the low energy effective Hamiltonian is indeed a spin 1 model. From the lowest energies $`E_S`$ for total spin $`S=0,1`$ and $`2`$, we have extracted $`J_1`$ and $`J_2`$. Our numerical results are in good agreement with the perturbation theory. In particular, we find a transition from antiferromagnetic to ferromagnetic $`J_1`$ as a function of $`\alpha =t_{23}/t_{22}`$, and a sharp increase in $`\beta `$ in the same region. This transition corresponds to the crossings in Fig. 4a. The exact numerical results for $`\beta `$ as function of $`\alpha `$ for typical parameters are also plotted in Fig. 4b. Except for the precise value of $`\alpha `$ for which $`\beta `$ diverges, the essential features found in the exact numerical calculations are the same as in the perturbation theory. Let us note that the numerical results do not depend qualitatively on the details of the parameters as long as the fundamental processes are present: There is always a value of $`t_{23}/t_{22}`$ around 1 at which $`\beta `$ diverges. So although the simplified model presented here in order to have compact expressions for the fourth order perturbation results ($`t_{11}=t_{12}=t_{13}=0`$, $`t_{33}=t_{22}`$) is not exactly realized in LiVGe2O6, we have checked numerically that the appropriate extension gives the same behaviour.

The relative value of the biquadratic interaction $`\beta `$ depends sensitively on the ratio of the hopping integrals, but not on their magnitude. Application of a pressure on the system, if it is uniformly added, may not change the ratio of these hopping integrals, hence the biquadratic interaction sensitively. A uniaxial pressure, on the other hand, may change the hopping integrals in different proportions, and could dramatically affect the magnetic properties of the system. This possibility would be worth investigating in the case of LiVGe2O6.

Finally, let us put these results in perspective. Although this mechanism would apply in any dimension, there is something very special about 1D. In higher dimension, fourth–order perturbation theory produces other types of four-spin interactions involving spins at four different sites. However, in 1D these processes do not appear as long as the hopping to further neighbours can be neglected, which usually is the case. So this mechanism is somehow specific to 1D. This might be the reason why it was overlooked in early attempts at producing strong biquadratic interactions because they were done before Haldane’s conjecture and were more concerned with 3D magnets.

This mechanism is also somehow related to the problem of orbital degeneracy. It is well known that orbital degeneracy leads to a Hamiltonian which is not purely Heisenberg, but involves a pseudo-spin for the orbital degrees of freedom. For spin 1/2 systems, for which most of the work has been done, the bilinear interaction is the only candidate - higher order terms do not appear because they can be rewritten in terms of the bilinear interaction - and the situation is a clear cut: Either the orbitals do not play a role, and the Hamiltonian is purely Heisenberg, or they do, and the Hamiltonian is a spin–orbital model. For spin 1 systems, we showed that, on going from a Heisenberg to a spin–orbital model by pulling down an empty orbital, there is an intermediate region where the Hamiltonian is still a pure spin Hamiltonian, but with a more general interaction than just bilinear due to a non-trivial level crossing (Fig. 4a). This simple idea should be useful in other contexts as well in looking for general spin Hamiltonians. In that respect, we hope that the present results will encourage quantum chemists to look more closely at exchange in situations where orbitals are quasi-degenerate. Quantum chemistry has proven to be able to reproduce very accurately on the basis of ab initio calculations the values of the exchange integrals of several systems. It should thus be possible along the same lines to specify more precisely the conditions under which the mechanism proposed in this paper will apply.

In conclusion, we have provided a simple but efficient mechanism to produce significant biquadratic interactions in spin 1 chains. It is our hope that these ideas will ultimately lead to an experimental investigation of the fascinating phase diagram of the general spin 1 model in one dimension.

NOTE: After completion of this project, we became aware of a paper by Bhatt and Yang \[R. N. Bhatt and K. Yang, J. Appl. Phys. 83, 7231 (1998)\] who considered a similar problem in the context of random antiferromagnetic spin chains. Their model is slightly different from ours (they consider the case of $`N`$ degenerate orbitals with a simplified interaction), and they restrict themselves to a perturbation calculation, but our results are qualitatively consistent with theirs.

ACKNOWLEDGEMENTS: We acknowledge useful discussions with T. Jolicoeur, K. Le Hur and O. Golinelli. This work was supported in part by DOE Grant No. DE/FG03-98ER45687.