text:
# 1 Introduction

## 1 Introduction

Since the discovery of quasicrystals (QC’s) , it has been an important issue in condensed-matter physics to find out which properties distinguish this type of material, which has long-range order but no lattice periodicity, from crystals and amorphous materials. Theoretically this problem was studied by means of so-called approximants. An approximant is a lattice periodic crystal with a large unit cell with a structure close to that of a QC. For icosahedral QC’s approximants are obtained if in a higher dimensional description of a QC , the golden mean $`\tau =\frac{1}{2}(\sqrt{5}+1)`$ is approximated by a fraction. By approaching $`\tau `$ in a systematic way, a series of approximants can be obtained with ever larger unit cells and ever better structural convergence with respect to the ideal QC. By studying such a series of approximants and meanwhile monitoring convergence with respect to unit-cell size, results may be obtained that are valid for an ideal QC. Due to computational limitations most first principles calculations  were restricted to so-called $`1/1`$ approximants, where the fraction indicates the approximation made for $`\tau `$, with, depending on the particular QC, approximately 150 atoms per unit cell. These calculations were done by means of the linear muffin-tin orbital (LMTO) method. They led to the surprising prediction that the electronic density of states (dos) of QC’s is spiky, a property that would be specific for QC’s. It started in 1989, when Fujiwara  predicted on the basis of a calculation for 1/1 AlMn that the dos of crystalline and quasicrystalline AlMn consists of a set of spikes. In 1993, he found that the electronic dos of $`1/1`$ AlMn  and $`1/1`$ AlCuLi  has dense spiky peaks of width 10-20 meV or less, which in his view suggested that spikes are generic in QC’s and their crystalline approximants. In the same year, Hafner and Krajčí  conjectured on the basis of results for $`1/1`$ AlZnMg that the dos of an ideal QC is spiky. In 1994, De Laissardière and Fujiwara  found spikes of width 14 meV in $`1/1`$ ALCuFe. They stressed that the spikes are in the first place a consequence of quasiperiodicity, accentuated by the presence of transition metal atoms. In 1995, Krajčí and co-workers  found a spiky dos for $`1/1`$ AlPdMn and $`2/1`$ AlPdMn. The latter dos was obtained by folding the eigenvalues of the Hamiltonian at the $`\mathrm{\Gamma }`$-point with a Gaussian of width 25 meV. Apart from being interesting, it was indicated that spikes may also be relevant to electronic transport in QC’s . Disappointingly spikes have not been observed directly in experiments, although a sufficiently high resolution was obtained by amongst others Stadnik and co-workers  by means of ultrahigh energy resolution photoemission spectroscopy (5 meV), Dadydov and co-workers  in tunnelling experiments (1 meV) and Escudero and co-workers  in point contact and tunnelling experiments (0.35 meV). By convoluting the dos of $`1/1`$ AlCuFe  with a Gaussian of full width at half maximum of 31.6 meV, Stadnik and co-workers  showed that even at a relatively low resolution, spiky features should according to the prediction be observed near the Fermi level. In ref. References we gave an overview of explanations that have been proposed for the discrepancy between theory and experiment, which we will not repeat here. Instead, we will focus on two explanations. The first was proposed by Hafner and Krajčí . They indicated the possibility that spikes are specific for small periodic approximants of QC’s and do not persist in ideal QC’s. The second explanation is that spikes are an artefact of the calculation . We have studied these explanations by calculating the dos for the two-dimensional (2D) and the three-dimensional (3D) Penrose tiling, which are quasicrystalline model systems. The 3D Penrose tiling has icosahedral symmetry and is the basis framework for many realistic models of icosahedral QC’s . It is important to realise that these models do not take into account different kinds of electrons, nor do they represent realistic local atomic configurations in QC’s, but they are quasiperiodic, and as a consequence of their simplicity it is possible to calculate properties of high approximants. Whereas in most LMTO studies results extrapolated for the quasiperiodic limit are based on calculations for a $`1/1`$ approximant, we could calculate the dos of the 3D Penrose tiling with a resolution of 10 meV, a result valid in the quasiperiodic limit. In this letter, we show that spikes present in a small approximant of the 3D Penrose tiling, do not survive in the ideal QC, in agreement with the original suggestion of Hafner and Krajčí . This result falsifies the prediction that all QC’s have a spiky dos.

## 2 Model

For the 2D Penrose tiling, we have considered the same model as in ref. References. The 3D Penrose tiling  is made with the cut-and-project method . We consider three types of approximants. First of all, if in the internal space $`\tau `$ is approximated by $`F_{n+1}/F_n`$ with $`F_i`$ a Fibonacci number, then a so-called a periodic $`F_{n+1}/F_n`$ approximant is obtained. If $`F_{n+2}`$ is even, there is a body centring . Second, if in the internal space an additional translation is made along a vector ($`1/F_n`$) ($`\frac{1}{2}`$, $`\frac{1}{2}`$, $`\frac{1}{2}`$), another $`F_{n+1}/F_n`$ approximant is formed, which has a non-symmorphic cubic space group, e.g., Pa$`\overline{3}`$. We will refer to these approximants as cubic approximants. Finally, taking different approximations for $`\tau `$ in different directions in the internal space, we consider $`\tau _x=F_n/F_{n1}`$, $`\tau _y=F_{n+1}/F_n`$ and $`\tau _z=F_{n+2}/F_{n+1}`$. This gives non-cubic ($`F_n/F_{n1}`$, $`F_{n+1}/F_n`$, $`F_{n+2}/F_{n+1}`$) approximants, of which the lattice is body-centred orthorhombic. An important parameter for each approximant is $`N`$, the number of vertices of the approximant tiling per unit cell. We studied the so-called vertex model, which has one electronic basis state on each vertex and nearest neighbour hopping via the edges of the tiling. We chose the on-site energy to be zero. The electronic basis states can be thought of as s-electrons. The Hamiltonian is:

$$H=\underset{nm}{}|nT_{nm}m\left|+\underset{n}{}\right|nV_nn|,$$
(1)
with $`V_n=0`$ and $`T_{nm}=1\mathrm{eV}`$ if the vertices $`n`$ and $`m`$ are connected by an edge of the tiling and zero otherwise. This model has been studied before , but so far it was not possible to draw definitive conclusions about the nature of the spectrum. Krajčí and Fujiwara  have shown that at $`E=0`$ there are infinitely many, strictly localised states. An important symmetry of this model is the bipartite property , as a consequence of which the density of states is symmetric with respect to $`E=0`$ . Also, for the non-body-centred cubic and for the non-cubic approximants, the bipartiteness of the lattice can be used to reduce the problem of finding eigenvalues of $`H`$ from an $`N\times N`$ problem to an $`N/2\times N/2`$ problem .

## 3 Method

In order to determine the dos, we perform a bandstructure calculation. Here $`N_k`$, the number of $`𝐤`$ points calculated in the asymmetric unit, is important. For the evaluation of the dos, we follow roughly the same scheme as in ref. References: First we calculate the exact average dos of a linearly interpolated  bandstructure of an approximant in intervals of width 0.2 meV, where the averaging ensures uniform convergence . Then we convolute with a Gaussian of full width at half maximum of typically 10 meV, which reduces the error in the dos by smearing away all features on a scale smaller than 10 meV including effects of the finite averaging interval of 0.2 meV, but which does not smear away spikes of width 10-20 meV. This last step is not general practice, but we will argue that it is useful. Maximal errors for the dos of approximants were extrapolated by making a systematic comparison of the results for various $`N_k`$ and simultaneously assuming asymptotic behaviour: the error in the dos at a certain energy is the integral between the true constant energy surface and the polyhedron approximating it , which in three dimensions converges as $`N_k^{2/3}`$  for sufficiently large $`N_k`$. For an ideal tiling, the maximal error was extrapolated assuming that the asymptotic behaviour of the error is proportional to $`1/N`$ in a given series of approximants , which holds provided that the error with respect to $`N_k`$ is negligible, which was checked explicitly.

## 4 Results

Figure 1

shows the dos of the cubic 3/2 approximant of the 3D Penrose tiling with $`N=576`$ and $`N_k=1376`$. It has space group Pa$`\overline{3}`$. A grid of so-called special $`𝐤`$ points  was used. Figure 1a) shows the result before convolution with a Gaussian: there are many densely distributed spikes, but the error is so large that it is not possible to say whether these spikes are real or an artefact connected to the small value of $`N_k`$. After smoothening, fig. 1b) is obtained. The error is considerably smaller, so that conclusions can be drawn about the spikes. Comparing the error and the amplitude of some of the spiky features, we see that certainly some of these features are real. So, the dos of small approximants of icosahedral QC’s may have sharp features. In fig. 3,

the dos of the ideal 2D Penrose tiling is shown. It was calculated with $`N=167761`$ and $`N_k=3`$, where $`𝐤`$ points were chosen that give real phase factors. The error indicates convergence with respect to both $`N`$ and $`N_k`$. Our conclusion of ref. References, that the 2D Penrose tiling has a spiky dos, is confirmed. In fig. 3, the dos of the 3D Penrose tiling is shown. Actually, it is the dos of the non-cubic (8/5, 13/8, 21/13) approximant with $`N=21892`$ and $`N_k=260`$, where the $`𝐤`$ points lie on a regular grid which includes the $`\mathrm{\Gamma }`$ point, but convergence with respect to $`N`$ was checked. As far as our data can tell, this result does not depend on the sequence of approximants that is considered. For example, the 8/5 approximant with $`N=10336`$ and $`N_k=8196`$ has a dos that deviates less than 0.008 from the graph of fig. 3 for all energies except $`E=0`$. Unfortunately, the 13/8 approximant with $`N=43784`$ is too large for us in order to be able to compute a sufficient number of $`𝐤`$ points. Since the dos of the ideal 3D Penrose tiling in fig. 3 appears to be quite smooth, we conclude that the spikes of fig. 1b) do not survive in the quasiperiodic limit. Figure 4a)

\[t\] gives the dos of the non-cubic (3/2, 5/3, 8/5) approximant of the 3D Penrose tiling with $`N=1220`$ and $`N_k=36`$. The number of electronic basis states per unit cell and $`N_k`$ were chosen as close as possible to those of the calculation for icosahedral AlCuFe , and we chose the same width for the Gaussian as Stadnik and co-workers  did when smoothening it. Figure 4b) shows the same dos with $`N_k=16388`$. Comparing the two graphs, we see that taking into account a too small number of $`𝐤`$ points may well lead to quite broad ($``$ 32 meV) additional features in the dos. These features then are an artefact. Because of the indicated error in fig. 4a) it is clear these features should not be taken seriously. The prediction for AlFeCu of ref. References does not indicate an error and from the published data it is not clear whether the predicted spikes are an artefact or not. A too small $`N_k`$ may well explain some of the spikes predicted in literature. Expressly, we want to remark that it is in principle legitimate to try to keep the number of $`𝐤`$ points that are calculated, small in order to save limited computer resources, but one has to realise the price that is paid: either the error is larger or the data should be presented with a worse resolution.

## 5 Conclusions

We have carried out an analysis of the influence of the size of the approximant unit cell and the number of $`𝒌`$ vectors in the Brillouin zone on the finer features of the electronic dos in an icosahedral tight-binding model. We have shown that the dos of the 2D Penrose tiling is spiky (fig. 3), whereas the dos of the icosahedral 3D Penrose tiling is smooth (fig. 3), at the given resolution of 10 meV. We have also shown how a spiky dos would be predicted for our 3D model when (i) only a small periodic approximant would have been studied (fig. 1b) or (ii) a too small number of $`𝐤`$ points would have been taken into account in the calculations, leading to large numerical errors (fig. 4). Possibly the published calculations predictions that spikes are generic for icosahedral QC’s, have not sufficiently taken these two points into account. Our results further suggest that although experimentally spikes could not be observed in QC’s, they might be found in approximant phases.

## Acknowledgments

This work has been supported by the Stichting voor Fundamenteel Onderzoek der Materie with financial support of the Nederlandse Organisatie voor Wetenschappelijk Onderzoek. We used free software from LAPACK . We thank Gilles de Wijs for discussions.


format: text

text:
# Collective oscillations in quantum rings: a broken symmetry case

## 1 Introduction

Quantum rings have attracted much interest in recent years. They constitute a clear example of the great technological advances in the fabrication of electronic nanostructures, which also include a rich variety of quantum dots. Theoretical studies of rings have addressed the appearance of the so called ‘persistent current’ as a function of an applied magnetic field, as well as the electronic structure and optical excitations within a fully microscopic approach for small rings Chak94 ; Chak96 ; Wend96 . On the other hand, semiclassical methods have been applied to larger rings Zar96 , for which there exist experimental results in the mesoscopic domain Dha93 . These semiclassical methods have indeed provided a good description of the $`B`$-dispersion in the measured far-infrared collective excitations for large rings. At an intermediate theoretical level, microcopic density-functional calculations have been applied to describe rings containing around $`N10`$ electrons Rei99 ; Emp99 . An exciting prediction of these latter calculations is the transition for increasingly narrow rings (approaching the quasi-1d limit) to a state with broken rotational symmetry in the intrinsic reference frame containing a static spin-density wave (SDW). The fabrication of such rings, with a small number of electrons, is a very demanding challenge that only recently has been addressed by Lorke et al., who reported measurements of the far infrared response of two electron rings Lor99 ; Lor98 .

All calculations of ring excitations for intermediate sizes ($`N10`$), using either semiclassical methods or density-functional theory (DFT), have imposed circular symmetry to the system and, therefore, have not been able to signal the influence of the broken symmetry ground state on the spectra. In this paper we relax the symmetry condition on the ground and excited states of rings with the purpose to explore these effects within DFT. To this end we will use real time simulations of the collective density oscillations developed by two of us in the context of quantum dots Pue99 ; Ser99 . The paper is organized as follows: in Sec. II we consider the ring ground-state within our approach, focussing on the transition to the symmetry-broken ground state; in Sec. III we address the real time simulation of the charge and spin-density oscillations. Finally, the conclusions are drawn in Sec. IV.

## 2 DFT description of rings

In order to model two-dimensional (2d) quantum rings of varying widths we have used the model proposed by Reimann et al. Rei99  and consider a ring confining potential on the $`xy`$ plane depending only on the radial distance $`r`$, such as

$$V^{(ext)}(r)=\frac{1}{2}\frac{\omega _0^2}{\delta }(rR_0)^2.$$
(1)
In this expression $`R_0`$ corresponds to the ring radius and we take $`\omega _0^2=1/r_s^3N_p^{1/2}`$, with $`N_p`$ and $`r_s`$ parameters analogous to those giving for a circular quantum dot the number of positive charges and the Wigner-Seitz radius, respectively. The $`\delta `$ parameter is introduced in order to model the radial thickness of the ring. As delta decreases, the ring becomes thinner and thinner, approaching the quasi-1d limit.

We describe the electronic structure within the local-spin-density approximation of DFT. The reader is referred to Refs. Kos97 ; Pi98  for details of the approach in the present context. The set of single-particle (sp) orbitals is obtained by selfconsistently solving the Kohn-Sham (KS) equations

$$\left[\frac{1}{2}^2+V^{(ext)}(r)+V^{(H)}(𝐫)+V_\eta ^{(xc)}(𝐫)\right]\phi _{i\eta }(𝐫)=\epsilon _{i\eta }\phi _{i\eta }(𝐫),$$
(2)
where $`V^{(H)}`$ and $`V_\eta ^{(xc)}`$ are respectively the Hartree and exchange-correlation potentials. The latter one is obtained from the local energy density as

$$V_\eta ^{(xc)}=\frac{_{xc}(\rho ,m)}{\rho _\eta },$$
(3)
In these expressions $`\eta =,`$ labels the spin components while the spin densities are $`\rho _\eta (𝐫)=_i|\phi _{i\eta }(𝐫)|^2`$. Total density and spin magnetization are given by $`\rho =\rho _{}+\rho _{}`$ and $`m=\rho _{}\rho _{}`$, respectively. Our exchange-correlation functional is based on the Tanatar-Ceperley results for the 2d electron gas Tan89  with the von Barth-Hedin interpolation for intermediate polarizations Bar72 .

We emphasize here that our solution to Eq. (2) is not assuming a priori any symmetry for the orbitals $`\phi (𝐫)`$. The 2d plane is discretized into a uniform grid of points in Cartesian coordinates and the Laplacian operator is approximated by the corresponding finite differences (we have used typically 7 point formulas). The KS solutions are then obtained by iteratively applying the imaginary-time-step method. This technique is very robust although convergence may be rather slow in some cases. The lowest energy solution in each case is obtained by performing calculations with different starting points in order to assure that the result does not correspond to a local minimum representing an excited state of the system.

### 2.1 Circular rings

Figure 1 shows the density and magnetization for a ring with units  $`R_0=8`$ , containing $`N=10`$ electrons and for two values of the width $`\delta `$. The wide one ($`\delta =2`$) has perfect circular symmetry and its magnetization vanishes. On the contrary, the narrow one ($`\delta =0.08`$) presents a modulation in charge density, a charge density wave (CDW), and spin polarization of antiferromagnetic type, with alternate orientations of spin up and down (a SDW). This drammatic change of behaviour is associated with the increase of the radial confinement and agrees with the findings of Ref. Rei99 . We compare in Fig. 2 the energy of the unrestricted solution to Eq. (2) with that obtained with the contraint of circular symmetry, i.e., by solving only the radial KS equation as in Ref. Emp99 . Also shown in Fig. 2 is the ratio of maximum magnetization to maximum density and the amplitude of the CDW. Panels (b) and (c) show that the transition is rather abrupt with the value of $`\delta `$, although the energy (panel (a)) is smoothly changing. Nevertheless, the energy trend is very clear and manifests the gain in binding energy with the formation of CDW’s and SDW’s.

Although some controversy has existed in the quantum dot community Kos97 ; Hi99 ; Yan99 , the formation of broken-symmetry solutions as those displayed by mean-field theories is a well known phenomenon in nuclear physics Rin90  and in atomic physics Ber89 . It has been recently discused in the context of 2d quantum dots by Yannouleas and Landman Yan99  and by Koskinen et al Kos97 ; Kos00 . It has been shown that the mean field solution corresponds to the intrinsic structure of the system although the full exact solution in the laboratory frame preserves the symmetry because of an underlying degeneracy. This is convincingly shown by comparing internal structure properties of the exact solution, such as the conditional pair probability, with the mean field resultYan99 . Another strong evidence of the internal symmetry breaking is given by the comparison of the low lying states with the roto-vibrational states of the structure formed by the localized electrons Yan99 ; Kos00 . Quantum rings constitute indeed a very good scenario for this behaviour since, as shown in Figs. 1 and 2, they exhibit strong breaking of the circular symmetry. A critical comparison of the mean field solution with the exact one for two-electron quantum rings will be presented in a future publication.

The broken-symmetry solution is accompanied by a strong change in the sp level scheme. While the circular structure is characterized by shells, with closings at electron numbers of 2, 6, 10, 14 and a small sensitivity on the width $`\delta `$, the formation of SDW’s and CDW’s leads to a bunching of the occupied levels and the formation of a relatively large energy gap at the Fermi level. This behaviour is shown in panel (d) of Fig. 2.

### 2.2 Elliptic rings

Deformed rings can be simulated by using an elliptic contour for $`R_0(\theta )`$ of the type

$$R_0(\theta )=a\left(\mathrm{cos}^2\theta +\frac{1}{\beta ^2}\mathrm{sin}^2\theta \right)^{\frac{1}{2}},$$
(4)
where $`a`$ is the major axis of the ellipse, $`\beta `$ gives the ratio of minor to major axis and $`\theta `$ is the polar angle. In order to compare with the circular rings of the preceding subsection we have taken elliptic rings with the same perimeter and have considered two deformations, namely $`\beta =0.75`$ and $`\beta =0.5`$, corresponding to medium and large ring deformation, respectively.

One of the most striking findings in elliptic rings is the appearance of the broken-symmetry solution along the ring contour for increasing deformation, at fixed width. This happens even for wide rings although the wave structure is more clearly marked in the narrow ones (Fig. 3). Another behaviour which is apparent in Fig. 3 is the formation of charge concentrations at the ends of the long axis. These are attributed to end states similar to those found in finite wires Rei99 . In broad rings, the formation of end states appears abruptly after a minimum deformation is reached. For instance, at $`\delta =2`$ and $`\beta =0.75`$ no evidence is yet found of their formation.

## 3 Time dependent spin and charge density oscillations

The analysis of excited states is often made by invoking perturbation theory for an external probing field. The feasibility of this approach usually relies on the symmetries of the system, which greatly simplify the treatment. For instance, in circular dots one can take density oscillations of the type $`\delta \rho _\eta (r)e^{im\theta }`$, i.e., a radial function times a multipolar field, and then reduce the equations to the radial parts. In the present context the possibility of having a broken symmetry ground state forbids this approach.

A suitable alternative to perturbation theory is to study the density oscillations by using real-time simulations. This is particularly well adapted to DFT and leads to its time-dependent generalization (TDDFT). The simplest version of this scheme is the adiabatic-local-density approximation, which uses the energy functional of the ground state to explore small amplitude oscillations around it. This is one of the simplest versions of TDDFT and yet it is quite robust and satisfies known exact properties as the energy weighted sum rule and the generalized harmonic-potential theorem Dob94 ; Lipp97 . We will use it to calculate the ring spin-density oscillations.

The time-dependent Schrödinger equations can be written

$$i\frac{}{t}\phi _{i\eta }(𝐫,t)=h_\eta [\rho ,m]\phi _{i\eta }(𝐫,t),$$
(5)
where the KS single-particle Hamiltonian $`h_\eta [\rho ,m]`$ is given by the square bracket in Eq. (2), whose stationary solution was found in the preceding Section. Here we will perform an initial (at $`t=0`$) pertubation of the wave functions

$$\phi _{i\eta }(𝐫,0)=𝒫\phi _{i\eta }(𝐫),$$
(6)
and track the time evolution through Eq. (5). Technical details of the integration method, as well as of the analysis of the dipole signals $`𝐫_\eta `$ which allow to extract the excitation spectrum can be found in Refs. Pue99 ; Ser99 .

We will distinguish three different types of excitation modes, two of them depending on the initial perturbation $`𝒫`$ and another one corresponding to a theoretical model in which the electrons are not interacting. Namely: density modes associated with an initial rigid displacement of the total density; spin modes for which at $`t=0`$ the spin densities are rigidly shifted in opposite directions; and free modes for which the sp Hamiltonian is kept fixed to the stationary densitites, i.e., $`h_\eta [\rho ,m]h_\eta [\rho _0,m_0]`$. The free oscillations, also known as sp modes, are equally excited by shifting either the total density or the spin densities.

### 3.1 circular rings

Figure 4 shows the free, density and spin spectra for the broad ring with $`N=10`$ and $`\delta =2`$. As known from other calculations Zar96 ; Emp99  the spectrum is roughly divided in two regions in the sp model which are then shifted by the interaction in the density and spin channels. The lower region contains a single peak while the higher one exhibits an important fragmentation. The density response is characterized by the blue shift from the sp peaks while the spin response exhibits a shift to lower energy. The residual interaction is also quite effective in enhancing one peak above the others in the high energy region. The peak that collects the strength is then associated to the more collective excitation. These features can be considered standard in finite Fermi systems with a repulsive interaction in the density channel and an attractive one for the spin channel. Indeed, a similar behaviour has been found in quantum dots. Figure 4 also shows the corresponding spectra obtained by using the circularly constrained solution in the perturbative random-phase approximation (RPA) Emp99 . It is obvious that both results are equivalent, as one would expect since the ground state has circular symmetry. We attribute the minor diferences to the quite different techniques that have been used in the two cases. It is worth to notice that the real time method is much more demanding computationally and that very fine details of the spectrum are harder to describe, such as very low intensity peaks, or the resolution of two closely lying excitations.

The results for the narrow ring ($`\delta =0.08`$) are shown in Fig 5. Qualitative differences appear with respect to the preceding ring and, most important, with respect to the circularly constrained solution. They must be attributed to the symmetry breaking ground state of this ring. The separation of high and low energy region is now clearer, with an important gap between both. The most striking result comes from the comparison of the low energy peaks with the circularly contrained solution. First we notice that the full solution in this energy region yields a density response which is very close to, but just below, the sp peaks and overlaps with some excitations of the spin response. In the circularly contrained calculation the density response is above the sp one and, furthermore, the spin peak lies at an essentially vanishing energy. The different behaviour is indicating that the symmetry breaking is giving stability to the ring, which would otherwise become unstable against spin oscillations. This is in agreement with the appearance of an energy gap in the sp energies associated to the symmetry breaking, mentioned in Sec. II.

The high energy region is associated to a locally radial oscillation in the ring, and thus similar to the oscillation in a one dimensional oscillator. This is confirmed by the similarity of the energy with the radial curvature at the minimum given by $`\omega _0/\sqrt{\delta }`$ that, by the well known generalized Kohn’s theorem, is the only allowed energy for rigid density oscillations of dipole type in parabolic confinement Mak90 . On the other hand, the low energy region is associated with oscillations locally tangential to the ring. Within this interpretation, the fact that the density and spin tangential modes are close in energy and essentially overlap with the sp excitations can be understood as a signature of a 1-dimensional Luttinger liquid Hal81 , for which the lowest excitations are not the sp but the collective ones.

In the symmetry broken ground state there is also a possibility to excite pure tangential modes by means of a twist of the two spin densities in opposite directions. These exotic modes were proposed recently in the context of quantum dots Pue99 . It is obvious that pure tangential modes can not exist in circularly symmetric rings since a twist of spin densities does not change the energy. However, in a symmetry broken state a restoring force appears because of the distortion of the spin-density wave. The associated oscillation energy is easily obtained by analyzing the time dependence of the orbital currents $`_i\mathrm{}_i^{(z)}\sigma _i^{(z)}`$ that appear after an initial spin twist.

Figure 6 shows the spin twist spectrum for the narrow ring with 10 electrons. We observe that the peaks are in the same region as the low energy dipole excitations, which is proving that the latter ones are indeed of tangential nature. Both free and interacting spectra are characterized by a rather regular energy spacing, after an initial gap. We attribute this to the parabolicity of the potential in the tangential direction for points close to the local minima.

### 3.2 Elliptic rings

The systematics with deformation is very similar to that discussed for the circular case, with the additional result that both dipole energy regions are now more fragmented. Figure 7 shows the corresponding spectra for a strongly deformed ring $`\beta =0.5`$ with thicknesses given by $`\delta =2`$ and $`\delta =0.08`$. This fragmentation is more regular in the wide rings, with a doubling of high and low energy peaks by the same energy separation. In the narrow ring the relative fragmentation is lower with respect to the separation of high and low energy regions which is in agreement with the previously mentioned transverse and longitudinal character of both regions.

## 4 Conclusions

The ground state and oscillations of quantum rings have been discussed within density-functional theory relaxing the constraint of circular symmetry of the electronic density. A transition to a broken symmetry ground state with static charge and spin density waves has been found for increasingly narrow rings, in agreement with previous results. The energy gain with respect to the circularly contraint solution, the amplitude of charge and spin-density wave as well as the sp level scheme have been discussed as a function of ring width. The transition to a wave structure along the ring contour has also been obtained in elliptic rings for increasing deformations. In this case, there also appear charge concentrations at the long axis ends indicating the formation of end states similar to those obtained in finite quantum wires.

Signatures of the broken symmetry ground state on the collective oscillation energies have been searched for in the dipole spin and density oscillations. We have shown that the broken symmetry is accompanied by an overlap of the density and spin response peaks in the low energy region corresponding to tangential oscillations. At the same time, the low lying sp excitations are at a slightly higher energy, which is indicating a quasi 1d behaviour similar to that of Luttinger liquids. The energy of these excitations is close to the gap associated to the broken symmetry. It has also been shown that the broken symmetry gives stability to the ring, which would otherwise become unstable against spin dipole oscillations. The oscillations in the radial direction are located at higher energy, close to the radial curvature of the parabola as indicated by the generalized Kohn’s theorem. The separation between high and low energy regions incresases as the ring width decreases.

Pure tangential modes in rings with broken symmetry, associated to a twist of the spin densities have been discussed and it has been shown that their energy is close to the low dipole peaks. This supports the interpretation of this dipole excitations. Besides, we have shown that the more exotic spin-twist modes are regularly spaced in energy because of a high degree of harmonicity of the potential. Finally, the dipole oscillations in elliptic rings have been obtained and their fragmentation depending of the ring width has been shown.

Further work is necessary in order to analyze in detail the convergence to the strict 1d limit, for very narrow rings, and how the peculiarities of a Luttinger liquid such as power-law (as opposite to long range) order in the SDW correlation function and spin-charge separation Hal81  may manifest in a finite size system.

## 5 Acknowledgements

M.V.-R. gratefully acknowledges support by the Consell de Mallorca-UIB. This work was performed under grant No. PB98-0124.


text:
# Parity Effect in Mesoscopic and Nanoscopic Superconducting Particles

## 1 Introduction

Superconductivity in metals is related to condensation of electrons into pairs which are extended objects with a size of order $`\xi \mathrm{}v_F/\mathrm{\Delta }_0`$. A fundamental question  is whether pairing can survive in small specimens with size $`d`$ less than $`\xi `$. In energy scale this corresponds to condition that $`\delta _1>\mathrm{\Delta }_0`$ where $`\delta _1=\mathrm{}v_F/d`$ is typical scale of dimensional energy quantization. The energy levels can remain highly degenerate in a perfect sample (see Fig. 1) but in a disordered sample, level spacing $`\delta _2\epsilon _F/N_s`$, where $`N_s`$ is the total number of sites (which is of the order of total number of electrons, $`N`$), is much smaller than $`\delta _1`$ and can put a limit to superconducting pairing at $`\delta _2\mathrm{\Delta }_0`$.

In recent experiments, signature of superconducting pairing was traced in ultrasmall specimens down to sizes of 10 nm. The superconductivity reveals itself through the number-parity effect, the dependence of superconducting energy on whether the number of electrons in sample is even or odd. Matveev and Larkin  introduced a parity gap $`\mathrm{\Delta }_P=\frac{1}{2}(E_{2n+1}+E_{2n1})E_{2n}`$ where $`E_N`$ is superconducting energy vs number of electrons $`N`$, to characterize that the paired state is preferable over the state with odd number of electrons. In a number of papers, in particular in the paper of Braun and Von Delft , $`E_N`$ was calculated assuming that in an odd-number sample, electron state near the Fermi energy is blocked from pairing, similar to the way in which finite-number systems are treated in nuclear physics . In present paper, we investigate the suppression of superconductivity in a finite size system by using an approximation similar to that of , and also in the opposite limit, by starting from atomic-scale system (small atomic cluster) for which an exact solution of pairing problem is provided, and trying to match small- and large-size approximations through an interpolation procedure.

## 2 Large-scale approximation

In the Bardeen-Cooper-Schrieffer (BCS) approximation  a variational density matrix $`\rho _{\alpha \beta }=\frac{1}{Z}\varphi _\alpha ^{}\mathrm{exp}((H\mu N)/kT)\varphi _\beta `$ is exploited to obtain the ground state energy with trial wave functions

$$\varphi _\alpha ^{(1)}=\underset{𝐤}{}(u_𝐤+v_𝐤a_𝐤^{}a_𝐤^{})|0\text{ and }\varphi _\alpha ^{(2)}=\underset{𝐤𝐤_F}{}(u_𝐤+v_𝐤a_𝐤^{}a_𝐤^{})a_{k_F}^{}|0$$
(1)
for even and odd number of electrons, respectively. In the BCS pairing scheme (case 1), assuming equally spaced levels with distance $`\delta `$, we receive $`\mathrm{\Delta }_{2m}=\mathrm{\Delta }_{2m+1}=\mathrm{\Delta }_0`$ whereas in the blocked pairing state (case 2) $`\mathrm{\Delta }_{2m}=\mathrm{\Delta }_0`$ and $`\mathrm{\Delta }_{2m+1}=\mathrm{\Delta }_0\delta /2`$. The energy of superconductor at $`T=0`$ is calculated as $`H\mu N=_n(\xi _n\sqrt{\mathrm{\Delta }^2+\xi _n^2})`$ resulting in a parity gap value $`\mathrm{\Delta }_P^{(1)}=\frac{\delta }{4}`$ (for $`\delta \mathrm{\Delta }_0`$) in the first case, and $`\mathrm{\Delta }_P^{(2)}\mathrm{\Delta }_0\mathrm{ln}\frac{\omega _D\delta }{\mathrm{\Delta }_0^2}+\frac{1}{2}\delta `$ in the second. Energy of the blocked state is larger than that of the BCS state and therefore there is no reason to believe that the former is a better approximation for pairing. There is no direct relation of the above calculation to the question that there is a large uncertainty in the number of particles, and assertion that number of particles in $`\varphi _\alpha ^{(2)}`$ is odd is merely an illusion.

## 3 Josephson effect in mesoscopic metal particles

Consider two identical small superconducting grains interacting through the Hamiltonian $`H_T=_{nm}T_{nm}a_n^{}b_m+H.c.`$ where $`a_n^{}(a_n)`$ creates (annihilates) an electron in the first grain in state $`n`$. Similarly $`b_n^{}`$ and $`b_n`$ are operators for the second grain. By introducing the phases $`\varphi _1`$ and $`\varphi _2`$ to the Bogoliubov $`u,v`$\- coefficients in Eq.(1), we calculate coupling energy due to $`H_T`$ as a function of $`\varphi _1`$ and $`\varphi _2`$ in a standard way  $`\mathrm{\Delta }E=\mathrm{\Delta }E_S^{}\mathrm{\Delta }E_S\mathrm{cos}(\varphi _1\varphi _2)`$ where

$$\mathrm{\Delta }E_S=|T|^2\mathrm{\Delta }^2\underset{|\xi _n|,|\xi _m|<\omega _D}{}\frac{1}{(\epsilon _n+\epsilon _m)\epsilon _n\epsilon _m}.$$
(2)
The stationary Josephson current at zero temperature is $`J_c\mathrm{sin}(\varphi _1\varphi _2)`$ where $`J_c=2e\mathrm{\Delta }E_S/\mathrm{}`$. Figure 2 shows dependence of Josephson critical current and energy gap as a function of gap suppression parameter $`\delta /\mathrm{\Delta }_0`$. The curves show that superconductivity is suppressed at level spacing $`\delta _c=3.56\mathrm{\Delta }_0`$. This however is a conclusion within the mean field approximation. It does not exclude the possibility that pairing may survive for $`\delta >\delta _c`$ as a gapless superconducting state or as a state of superconducting insulator with an effective insulator gap $`\delta `$ larger than the superconducting pseudogap.

## 4 Richardson-Sherman approximation

A finite size solution to the pairing problem was suggested by Richardson and Sherman  (and recently revived by von Delft ) assuming that states with energy within the interaction shell $`\omega _D<\xi _n<\omega _D`$ can be approximated as hard-core bosons whereas states away from this shell are treated exactly as non-interacting particles. The interaction Hamiltonian for hard core bosons is

$$H=2\delta \underset{n\alpha }{}nc_{n\alpha }^{}c_{n\alpha }\frac{V\delta }{d}\underset{|\xi _n|,|\xi _m|<\omega _D,\alpha \beta }{}c_{n\alpha }^{}c_{m\beta }$$
(3)
where boson operators are $`c_{n\alpha }=a_{n\overline{\alpha }}a_{n\alpha }`$, $`d`$ is level degeneracy, and $`\alpha `$ labels states within the degenerate energy level $`\xi _n=n\delta `$. Hamiltonian 3 is solved  by iteration scheme proposed by Richardson and Sherman . Dependence of $`\mathrm{\Delta }_P`$ on $`\delta /\mathrm{\Delta }_0`$ is shown in Fig. 3. In nondegenerate case ($`d=1`$) $`\mathrm{\Delta }_P(\delta )`$ shows a behavior similar to previous results with smaller number of levels. In case of degenerate energy levels, pairing gap is multi-valued depending upon which three successive states are used in calculation of $`\mathrm{\Delta }_P`$. The true pairing gap is possibly the largest one in Fig. 3(b), and others represent a new kind of mesoscopic effect which may or may not be related to superconductivity. Small $`\delta `$ of $`\mathrm{\Delta }_P(\delta )`$ does not match with any of pairing gaps calculated in Section 3, which means possibly that neither of the approximations (original BCS or BCS with one level blocked) are good for small systems.

## 5 Parity effect in atomic clusters

We calculate the exact eigenstates of the negative-$`U`$ Hubbard Hamiltonian

$$H=t\underset{<i,j>,\sigma }{}a_{i\sigma }^{}a_{j\sigma }U\underset{i=1}{\overset{N}{}}a_i^{}a_i^{}a_ia_i+\underset{i\sigma }{}V_ia_{i\sigma }^{}a_{i\sigma }$$
(4)
(with $`U>0`$) in small system of $`N`$ sites. $`V_i`$ is random potential at site $`i`$. The fact that Hamiltonian 4 commutes with the total number of electrons and also with the total spin leads to lowering of dimensionality $`2^{2N}\times 2^{2N}`$ of the matrix. Interaction is non-retarded (unlike BCS Hamiltonian) and momentum dependence of energy is nonparabolic. We do not believe however that these are crucial differences, and indeed the small-$`U`$ limit mean-field solution  corresponds to an energy gap $`\mathrm{\Delta }_012t\mathrm{exp}(12t/U)`$. Figure 4,a shows minimal energy $`E_N`$ of cubic cluster as a function of number of electrons $`N`$. Parity effect clearly manifests itself in lower $`E_N`$ values for even $`N`$ and, when plotting $`\mathrm{\Delta }_P`$ against $`\delta /\mathrm{\Delta }_0`$ (Fig.4,b), displays a behavior very similar to one received in the Richardson-Sherman model.

One of the authors (I.O.K.) acknowledges helpful discussions with the participants of International Symposium on Mesoscopic Superconductivity (NTT, March 2000, Atsugi, Japan), especially with A.F. Andreev, H. Takayanagi and R.I. Shekhter. This work was partially supported by the Scientific and Technical Research Council of Turkey (TUBITAK) under grant No. TBAG 1736 and by the National Research Council of Italy, under the Research and Training Program for the Third Mediterranean Countries.


text:
# Analysing Lyapunov spectra of chaotic dynamical systems

## Abstract

It is shown that the asymptotic spectra of finite-time Lyapunov exponents of a variety of fully chaotic dynamical systems can be understood in terms of a statistical analysis. Using random matrix theory we derive numerical and in particular analytical results which provide insights into the overall behaviour of the Lyapunov exponents particularly for strange attractors. The corresponding distributions for the unstable periodic orbits are investigated for comparison.

Finite-time Lyapunov exponents represent an important tool for the quantitative description of the geometrical and dynamical properties of chaotic dynamical systems  and apply therefore to a variety of different physical situations. Numerous recent works (see, for example, refs.) try to describe the spectrum of these exponents in terms of products of random matrices. The majority of these investigations focus on the maximum average Lyapunov exponent which can be estimated by taking the proper mean of an ensemble of chaotic trajectories for a sufficiently large time interval $`p`$. Details of the corresponding distribution of Lyapunov exponents (LE) are thereby not important. In order to analyze strange attractors in more detail it is however crucial to know the spectrum of the finite-time LE which provides valuable information on the structures and properties emerging in phase space . For the asymptotic (large $`p`$) Lyapunov spectra it is known  that, apart from a few exceptions , a Gaussian approximation fits the behaviour around the maximum very well. Very little is however known with respect to the overall behaviour of the distribution. This is in contrast to the fact that the non-Gaussian tails of the distribution have significant influence on physical processes .

In the present letter we analyze finite-time Lyapunov spectra for low-dimensional discrete dynamical systems with fully developed chaos for the case of large time intervals $`p1`$, i.e. their asymptotic form. Our main interest is twofold. First we explore common features of these spectra beyond their behaviour in the vicinity of their maxima and compare them to the corresponding spectra of the unstable periodic orbits (UPOs). It turns out that these features can be understood in terms of products of random matrices. Second we investigate the origin of spectral properties which depend on the dynamical system. System-dependent characteristics can e.g. be due to invariant structures in phase space like UPOs.

It is an open question how the random-like features of the chaotic dynamics determine the distribution of the finite-time LE. One way to investigate this is to replace certain dynamical quantities of the original system by random ones and to study the resulting changes and/or common properties in the spectra of LE. In several recent works  such random-like modifications have been suggested. Adopting this methodology here, we will use four different ensembles of trajectories leading to different distributions for the corresponding LE. Each ensemble $`E_i`$ consists of $`N`$ trajectories (typically $`N10^6`$) of a given length $`p`$ ($`p1`$) and is given by (i) the distribution of the starting points of the finite trajectories and by (ii) the rule applied to generate the trajectory itself. We will then be primarily interested in the distributions of the finite-time maximum LE $`\lambda =\frac{1}{p}\mathrm{ln}\mathrm{\Lambda }`$ belonging to these ensembles $`\{E_i\}`$. $`\mathrm{\Lambda }`$ is the largest of the absolute values of the eigenvalues belonging to the transfer matrix $`M_i^{(p)}={\displaystyle \underset{k=1}{\overset{p}{}}}M_i(\stackrel{}{x}_k)`$ with $`\{\stackrel{}{x}_k|k=1,\mathrm{},p\}`$ being a finite-time trajectory of the ensemble $`E_i`$. $`M_i(\stackrel{}{x}_k)`$ is the stability matrix belonging to the (chaotic) dynamical law $`\stackrel{}{x}_{n+1}=\stackrel{}{F}(\stackrel{}{x}_n)`$. The index $`i`$ of $`M_i`$ indicates the ensemble $`E_i`$ according to which the points $`\{\stackrel{}{x}_k\}`$ of the trajectories are determined. It is important to note that all ensembles $`\{E_i\}`$ (see below) use the specific functional form of the stability matrix belonging to $`\stackrel{}{F}`$ but involve trajectories with different degrees of randomness. In the following we specify the ensembles $`\{E_i\}`$.

$`E_1`$ consists of trajectories of length $`p`$ obtained via iteration of the chaotic dynamical law $`\stackrel{}{F}`$. The corresponding initial conditions are distributed according to the invariant density of the map. This gives the so-called finite-time Lyapunov exponent distribution (FTLED) (see ref. and refs. therein). The ensemble $`E_2`$ consists of trajectories which are generated by a random variable distributed according to the invariant density of the chaotic map $`\stackrel{}{F}`$. This yields the bootstrap Lyapunov exponent distribution (BLED) . Compared to the FTLED the BLED corresponds to a dynamics with enhanced random character. The successive points of the bootstrap trajectory are completely uncorrelated. The third ensemble $`E_3`$ uses a uniformly distributed random variable for the generation of the trajectories. The range of the uniform distribution is chosen according to the phase space of the dynamical system. This case corresponds to a random matrix simulation of the dynamical system however respecting the form of the stability matrix $`𝐌`$ which belongs to the map $`\stackrel{}{F}`$. The resulting distribution of the Lyapunov exponents is called the random matrix Lyapunov exponent distribution (RMLED). Within the present investigation the ensemble $`E_3`$ possesses the highest degree of randomness. The fourth ensemble $`E_4`$ consists of the UPOs of the dynamical system $`\stackrel{}{F}`$ with period $`p`$ and correspondingly the distribution of their maximal LE. For fixed $`p`$ the number of trajectories contained in $`E_4`$ is finite according to the topological entropy of the corresponding phase space.

Let us now explore the distributions defined above for a variety of low-dimensional fully chaotic systems. Our main goal is to analyse and understand the overall behaviour of the FTLED for these systems. We begin with a simple $`1d`$ example: the logistic map. Results on the FTLED for this system can be found in ref.. Noteworthy the FTLED has a non-Gaussian form with one dominating central cusp. In comparison to this our numerical calculations on both the BLED and RMLED show that they are smooth functions with a Gaussian-like maximum but with characteristic asymmetric tails. For the maximum Lyapunov exponent distribution of the UPOs the exact result is a $`\delta `$-function i.e. $`\rho _p(\lambda )=\delta (\lambda \mathrm{ln}2)`$ independent of the period $`p`$. Looking at the FTLED one observes that the tails of the distribution are rather similar to the tails of the BLED and the RMLED while the cusp (maximum) is located exactly at $`\lambda =\mathrm{ln}2`$. The BLED (or RMLED) reflects therefore the overall behaviour of the FTLED i.e. describes the envelope of the FTLED. The latter possesses an additional central peak at the position of the Lyapunov exponents of the UPO’s. These features will in the following turn out to be common for a broad class of dynamical systems.

We focus in the following on two-dimensional systems with a strange attractor. The Henon map  for $`a=1.4`$ and $`b=0.3`$ and the Ikeda map  for $`a=0.9`$ and $`b=6`$ are prototypes of such systems. The FTLED of both are presented in Fig. 1. The ensemble $`E_1`$ consists of $`10^6`$ trajectories of length $`p=27`$ for both the Henon and Ikeda map. The results are surprisingly similar to the $`1d`$ case: A characteristic envelope with asymmetric tails dominates the distribution while superimposed peaks indicate the presence of additional structures in phase space. As we shall see in the following three regions with different functional forms of the envelope can be distinguished in the FTLED as well as the RMLED and BLED: A fast asymptotic algebraic decay for large values of $`\lambda `$, a Gaussian-like behaviour around the maximum and a dominating exponential decay for sufficiently small $`\lambda `$. Figure 2 shows the BLED and RMLED for these maps for the same length $`p`$. The position of the maximum of the RMLED is sensitive with respect to the random number intervals chosen, i.e. its location carries the information of the position of the attractor in phase space. Based on the above results and observations we are naturally lead to the following conclusion: the basic possible features of the smooth envelope of the FTLED (asymmetric structure, asymptotic tail properties) of a chaotic dynamical system are of random origin and can be obtained and understood by a corresponding study of random matrices (see below). Additional superimposed structures are signatures of e.g. invariant sets in phase space and are therefore of exclusively deterministic dynamical origin.

To elucidate and quantify the above observations we perform in the following an analytical investigation of the RMLED. This will allow us to thoroughly understand the behaviour of the RMLED and consequently the corresponding aspects of the FTLED. We begin by introducing a fictitious dynamical system with a stability matrix $`𝐌_i`$ of strongly random character i.e.

$$𝑴_i=r_i𝑨=\left(\begin{array}{cc}\mathrm{a}r_i& \mathrm{b}r_i\\ \mathrm{c}r_i& \mathrm{d}r_i\end{array}\right)$$

where $`r_i`$ is a random variable uniformly distributed in $`[0,R]`$ and $`i`$ labels the fictitious trajectory of length $`p`$. The simple form of the matrix $`𝑴_i`$ allows us to factorize the random variables $`\{r_i|i=1,\mathrm{},p\}`$ of the transfer matrix $`𝑴^{(p)}={\displaystyle \underset{k=1}{\overset{p}{}}}𝑴_k`$ and to reduce the problem of the product of random matrices to that of a product of random numbers. The matrix structure is then retained in the constant matrix $`𝐀`$ which is assumed to be nonsingular.

The distribution of the maximum LE for trajectories of length $`p`$ of this system is determined as:

$$\rho _p(\lambda )=_0^R𝑑r_1_0^R𝑑r_2\mathrm{}_0^R𝑑r_p\delta (\lambda \frac{1}{p}\mathrm{ln}\underset{i=1}{\overset{p}{}}r_i|\mathrm{\Lambda }_{max}|)\underset{l=1}{\overset{p}{}}\stackrel{~}{\rho }(r_l)$$
(1)
where $`\stackrel{~}{\rho }(z)=\mathrm{\Theta }(z)\mathrm{\Theta }(Rz)\frac{1}{R}`$, $`\mathrm{\Theta }`$ being the step function and $`\mathrm{\Lambda }_{max}`$ is given by $`\mathrm{\Lambda }_{max}=\frac{1}{2}(\mathrm{Tr}𝑨+\mathrm{sign}(\mathrm{Tr}𝑨)\sqrt{(\mathrm{Tr}𝑨)^24\mathrm{d}\mathrm{e}\mathrm{t}𝑨})`$. Using the substitution $`t_i=\mathrm{ln}(r_i|\mathrm{\Lambda }_{max}|)`$ and performing the Fourier transform of the $`\delta `$-function involved in (1) we obtain:

$$\rho _p(\lambda )=\frac{p}{2\pi (R|\mathrm{\Lambda }_{max}|)^p}_{\mathrm{}}^{\mathrm{}}𝑑ke^{ikp\lambda }_{\mathrm{}}^{\mathrm{ln}(R|\mathrm{\Lambda }_{max}|)}𝑑t_1\mathrm{}._{\mathrm{}}^{\mathrm{ln}(R|\mathrm{\Lambda }_{max}|)}𝑑t_pexp\left((1+ik)\underset{i=1}{\overset{p}{}}t_i\right)$$
(2)
The integrations over $`t_i`$ in (2) can be easily performed leading to a $`p`$-th order pole in the complex $`k`$-plane at $`k=i`$. This pole structure, corresponding to values of $`k`$ for which the exponent of the second exponential term in (2) vanishes and for which the integration of $`t_i`$ leads to singularities, is responsible for the features of the LED described above. Complex integration finally yields:

$$\rho _p(\lambda )=\frac{p^p}{(p1)!}\left(\mathrm{ln}(R|\mathrm{\Lambda }_{max}|)\lambda \right)^{p1}e^{p(\mathrm{ln}(R|\mathrm{\Lambda }_{max}|)\lambda )}\mathrm{\Theta }(\mathrm{ln}(R|\mathrm{\Lambda }_{max}|)\lambda )$$
(3)
Eq.(3) demonstrates that the exponential behaviour dominates for sufficiently small values of $`\lambda `$. Around the maximum at $`\lambda _o=\mathrm{ln}(R|\mathrm{\Lambda }_{max}|)(1\frac{1}{p})`$ the saddle point approximation ($`p1`$) gives a Gaussian. For values of $`\lambda `$ close to the maximum value $`\lambda _{max}=\mathrm{ln}(R|\mathrm{\Lambda }_{max}|)`$ we arrive at a power law behaviour, i.e. an algebraic decay.

Although the fictitious dynamical model discussed above captures the main features of the statistical properties of the distributions of Lyapunov exponents it is clearly desirable to investigate chaotic dynamical systems for which the RMLED can be obtained analytically. To this end let us consider the dynamical system defined by the quadratic equations: $`x_{n+1}=ay_n^2+b;y_n=cx_n+d`$. This system possesses for $`a=c=d=1`$ and $`b=2.5`$ a strange attractor which contains repeating crosses of decreasing size. For that reason we call it in the following the cross map. The average maximum Lyapunov exponent of this attractor is $`\overline{\lambda }0.123`$ while its fractal dimension is $`d_F1.78`$. The RMLED for the cross map can be calculated analytically following the line described above for our fictitious model. One pecularity of the cross map is that one has to distinguish between the RMLED obtained through trajectories with even and odd length $`p`$. The reason is that for odd $`p`$ both eigenvalues of the transfer matrix $`M^{(p)}`$ have the same absolute value, while for even $`p`$ there are two eigenvalues different with respect to their absolute value and one has to select the maximal one. After a tedious calculation using complex contour integration techniques we find for the RMLED of the cross map for odd values of $`p`$ the result:

$$\rho _p(\lambda )=\frac{2p}{p!}\left[\frac{pR_2}{R_2R_1}\right]^pexp\left(p(\mathrm{\Delta }2\lambda )\right)\underset{j=0}{\overset{p}{}}\left(\begin{array}{c}p\\ j\end{array}\right)(\mathrm{\Delta }\frac{j}{p}\mathrm{ln}|\frac{R_2}{R_1}|2\lambda )^{p1}\mathrm{\Theta }(\mathrm{\Delta }\frac{j}{p}\mathrm{ln}|\frac{R_2}{R_1}|2\lambda )$$
(4)
Here we have taken the random variables appearing in the stability matrix of the cross map to be uniformly distributed in the interval $`[R_1,R_2]`$. The parameter $`\mathrm{\Delta }`$ is given as: $`\mathrm{\Delta }=\mathrm{ln}|2acR_2|`$. From eq.(4) we see that the RMLED of the cross map is essentially a product of a single exponential and a sum of power laws. The latter possess all the same power $`p1`$ and differ only with respect to the constants involved. The similarity to the RMLED result of our model system in eq.(3) is obvious which confirms the universality of certain features of the RMLED. The RMLED for the case of even $`p`$ is given by the integral:

$$\rho _p(\lambda )=2\rho _{\frac{p}{2}}(\lambda )_{\mathrm{}}^\lambda 𝑑z\rho _{\frac{p}{2}}(z)$$
(5)
with $`\rho _{\frac{p}{2}}(x)`$ according to eq.(4). The integration in (5) can be performed analytically leading to a lengthy expression which will be not given here. The main characteristics of the function $`\rho _p(\lambda )`$ are again the features stated previously: a dominating exponential behaviour for sufficiently small values of $`\lambda `$, a Gaussian maximum and a fast algebraic decay for $`\lambda `$ close to its maximum value. We have also studied the FTLED, BLED and the distribution of the Lyapunov exponents of the UPOs for the cross map for $`p=28`$. The results are shown in Fig. 3. The envelope of these distributions exhibit the same features as discussed above. The additional structures present in the FTLED of the cross map are, as can be seen from Fig. 3, due to the presence of the UPOs which provide signatures of deterministic dynamical origin.

Finally let us consider the distributions of the Lyapunov exponents of the UPO’s of the Henon ($`p=27`$) and the Ikeda ($`p=14`$) map  which are presented in Fig. 4. The signs of the characteristic properties of the envelopes discussed above are also visible here. An interesting feature appears for the distribution of the LE of the UPOs of the Ikeda map: it is shifted significantly compared to the FTLED. This shift is probably due to the fact that the UPOs fail to reproduce the invariant density of the attractor, at least up to the above-considered period. It is well-known that the Ikeda attractor needs a description going beyond the linear neighbourhood .

Summarizing our results we have demonstrated that the overall behaviour of the finite-time Lyapunov exponent distributions of fully chaotic dynamical systems show general characteristics i.e. they can be understood in terms of statistical random matrix simulations of the systems. Seemingly this holds also for the distributions of the Lyapunov exponents of the unstable periodic orbits embedded into the chaotic phase space. Since Lyapunov spectra are at the heart of our understanding of chaotic systems in general our results apply to a variety of different physical systems.

Financial support by the Deutsche Forschungsgemeinschaft and the Landesgraduiertenförderungsgesetz Baden-Württemberg (D.P.) are gratefully acknowledged. We appreciate valuable discussions with O. Biham.

Figure Captions


text:
# Time-temperature superposition in viscous liquids

## Abstract

Dielectric relaxation measurements on supercooled triphenyl phosphite show that at low temperatures time-temperature superposition (TTS) is accurately obeyed for the primary ($`\alpha `$) relaxation process. Measurements on 6 other molecular liquids close to the calorimetric glass transition indicate that TTS is linked to an $`\omega ^{1/2}`$ high-frequency decay of the $`\alpha `$ loss, while the loss peak width is nonuniversal.

Liquids approaching the calorimetric glass transition have very high viscosities and very long average relaxation times compared to, e.g., room-temperature water. The study of highly viscous liquids is fascinating because – independent of the chemical nature of the intermolecular bonds – these liquids share a number of common features . The most significant common features are: i) non-Arrhenius temperature-dependence of the average relaxation time, and ii) non-Debye linear response functions. This paper addresses an old empiricism relating to the latter, time-temperature superposition (TTS) .

The average relaxation time is basically the time it takes for the liquid structure to relax to equilibrium after a small perturbation. This quantity is easily determined from the frequency-dependence of quantities like dielectric constant , specific heat , shear modulus , or bulk modulus . Measurements are usually reported in terms of the frequency-dependence of the response function’s imaginary part, the “loss.” The inverse loss peak frequency is a measure of the average relaxation time .

When temperature is lowered any loss peak moves towards lower frequencies. This is a rather dramatic effect in the sense that a temperature-change of just a few percent of the calorimetric glass transition temperature $`T_g`$ usually changes the loss peak frequency by more than one decade.

Loss peaks are often plotted in log-log plots. If the shape of the loss peak in this plot is temperature-independent the liquid obeys TTS with respect to the response function in question $`\chi (\omega )`$. Mathematically, if $`T`$ is temperature and $`\tau (T)`$ the average relaxation time TTS is obeyed whenever functions $`N`$ and $`\varphi `$ exist such that $`\chi (\omega ,T)=N(T)\varphi \left[\omega \tau (T)\right]`$. The factor $`N(T)`$ is usually a weak function of temperature.

When TTS applies response functions are easily measured over many decades of frequency. This is done by combining measurements done at different temperatures. This procedure works even if only one or two decades of frequency are directly accessible (as in, e.g., mechanical relaxation probed by the torsion pendulum technique or some of the first ac specific heat measurements ). Clearly, TTS is extremely useful if it is correct. But is this the case?

TTS has a long history . For many years TTS was assumed uncritically and used extensively. As accuracy increased and frequency ranges widened, more and more systems were found not obeying TTS. It turned out that TTS-violations were not just exceptions but, in fact, quite common. These violations were initially explained as being due to interference from an additional minor relaxation processes, the so-called $`\beta `$ relaxation often found at frequencies higher than those of the $`\alpha `$ relaxation. Several problems remained, however. From a practical point of view, if TTS does not always work it cannot be assumed a priori. Thus the option of constructing wide-frequency master curves was no longer available. Also, additional relaxation processes may always be postulated to “explain” TTS-violations, but what is to be learned from this?

In this situation the nonpolymeric scientific community became skeptical of TTS (for polymers TTS is still often assumed ). For several years now there has been little interest in TTS for nonpolymeric viscous liquids close $`T_g`$. In the less viscous regime at higher temperatures where mode-coupling theory is believed to apply, this theory’s TTS-prediction has received some confirmation from both experiment  and computer simulations . From the mode-coupling perspective, breakdowns of TTS in the very high viscosity regime correlate with the well-known breakdown of ideal mode-coupling theory here.

We are mainly interested in the behavior of \[equilibrium\] highly viscous liquids just above $`T_g`$. If TTS is general this would give an important input to theory. Similarly, connecting TTS to other features of the relaxation would also be important. Our motivation for reinvestigating the validity of TTS close to $`T_g`$ is the following. Dielectric relaxation often exhibits a high-frequency tail a few decades above the $`\alpha `$ peak, known as the Nagel wing. The Nagel wing is temperature-dependent and plays an important role in TTS-violations . Although generally thought to be part of the $`\alpha `$ process , it has recently been suggested that the Nagel wing instead derives from an additional, partly hidden $`\beta `$ process . A few months ago this was confirmed by long-time annealing experiments on propylene carbonate and glycerol by Schneider, Brand, Lunkenheimer and Loidl  (for glycerol $`\beta `$ relaxation was reported at high pressures long time ago ). The obvious question now is: What happens at temperatures low enough that any influence of $`\beta `$ relaxation is eliminated – is it possible that $`\alpha `$ relaxation does obey TTS after all?

Our dielectric set-up, which is briefly described in Ref. , covers frequencies from 1 mHz to 1 MHz. The dielectric loss is determined with a precision better than $`10^4`$. The absolute temperature uncertainty is below 0.1 K with relative variations during a frequency sweep below 1 mK .

Figure 1 shows the dielectric loss for triphenyl phosphite \[Aldrich, 97%\]. Figure 1a gives the raw data, Fig. 1b shows the same data scaled to check TTS. Visually, TTS is obeyed. A closer look is provided by Fig. 1c which compares the difference \[noisy curve\] between two curves of Fig. 1b to the difference \[sinusoidal curve\] between one curve and its best Havriliak-Negami fit . We conclude that TTS is obeyed quite well. In fact, deviations are roughly as expected from the 1 mK temperature inaccuracy.

For triphenyl phosphite $`\alpha `$ and $`\beta `$ relaxations are unusually well separated. This is because the $`\beta `$ peak is at higher frequencies than for most other viscous liquids. Two things are to be noted from Fig. 1. First, TTS applies for the $`\alpha `$ relaxation whenever temperature is low enough that $`\beta `$ relaxation does not interfere. Secondly, the high-frequency slope of the $`\alpha `$ peak is close to $`1/2`$. How general are these findings? Is TTS universally obeyed whenever additional relaxations do not interfere? Is the high-frequency $`\alpha `$ slope close to $`1/2`$ whenever TTS is obeyed? These questions are looked into in Fig. 2.

Figure 2a shows how $`\alpha `$ peak widths vary with temperature for triphenyl phosphite ($``$) and 6 other molecular liquids. This is done  by giving $`w`$, the logarithm of the width at half maximum relative to a Debye process, as a function of loss peak frequency – the latter quantity providing a convenient measure of temperature. When TTS is obeyed $`w`$ is temperature-independent. Most liquids seem to converge to TTS at low temperatures (the same is seen in the data of Nagel and coworkers from 1990 ), but this is not without exceptions. In these cases (e.g., glycerol) TTS may be violated because of interference from additional relaxations, but it shall be hard to prove this.

How about the $`\alpha `$ loss peaks for the liquids that obey TTS, do they have anything in common? Clearly, the width is nonuniversal. If the $`\alpha `$ process follows a high-frequency power-law the exponent may be estimated by evaluating the minimum slope above the loss peak in the log-log plot. Figure 2b shows this minimum slope for the same 7 liquids. Comparing with Fig. 2a we find that the minimum slope is close to $`1/2`$ whenever $`w`$ is virtually temperature-independent.

Figures 2a and 2b suggest that for many liquids TTS is obeyed at sufficiently low temperatures and that, whenever this happens, the $`\alpha `$ loss decays as $`\omega ^{1/2}`$ at high frequencies . We have indications that this picture applies also for the frequency-dependent shear modulus (which we can measure between 1 mHz and 50 kHz ), but the shear data are more noisy than the dielectric data.

To conclude – generalizing the triphenyl phosphite observations – we find that

* TTS is often obeyed for the $`\alpha `$ relaxation at low temperatures .
* Whenever this happens the $`\alpha `$ loss decays as $`\omega ^{1/2}`$ at high frequencies.
* The $`\alpha `$ loss peak width is nonuniversal.

An extension of this is the following conjecture:

* There is always a $`\beta `$ relaxation , though it is sometimes hidden under the $`\alpha `$ process and only visible as a Nagel wing.
* Whenever temperature is so low that additional relaxations do not interfere, the $`\alpha `$ relaxation obeys TTS.

To prove this beyond any reasonable doubt, however, one would have to improve sensitivity at ultra-low frequencies and, for some liquids, measure over many years.

To some extent the above 3+2 points reflect the opinion prevailing in this research field 30 years ago, now generally regarded as obsolete. Thus, Johari and Goldstein originally conjectured that $`\beta `$ relaxation is “a characteristic property of the liquid in or near the glassy state”. In this situation, of course, any TTS-violation gets a logical explanation since $`\alpha `$ and $`\beta `$ relaxations have quite different temperature-dependencies . The $`\omega ^{1/2}`$ high-frequency $`\alpha `$ peak decay was also in focus at that time and several theories were proposed to explain this .


text:
# LA-UR-00-2839 Geometric Morphology of Granular Materials

## 1 INTRODUCTION

Mathematical morphology was born in 1964, when G. Matheron  was asked to investigate the relationships between the geometry of porous media and their permeabilities, and when at the same time J. Serra  was asked to quantify the petrography of iron ores, in order to predict their milling properties. This initial period (1964 – 1968) has resulted in a first body of theoretical notions. The notion of a geometrical structure, or texture, is not purely objective. It does not exist in the phenomenon itself, nor in the observer, but somewhere between the two. Mathematical morphology quantifies this intuition by introducing the concept of structuring elements. Chosen by the morphologist, they interact with the object under study, modifying its shape and reducing it to a sort of caricature which is more expressive than the actual initial phenomenon. The power of this approach, but also its difficulty, lies in this structural analysis.

In this paper, we are going to present a new method of structural analysis for images (micrographs) of granular materials to the immediate benefit of the field of micromechanics of materials. In this particular field of physics, hydrodynamic model calculations are performed to investigate the propagation of shockwaves through such heterogeneous media. To be specific, one needs to setup the calculations with very accurate initial conditions, e.g., distributions of grain sizes in case one is interested in the mechanical properties of granular materials. In our presentation we demonstrate how to prepare the micrographs for extraction of proper input statistics data by using low and high level image processing algorithms.

We perform image segmentation of a micrograph (original image and its smoothed version) by applying either pulse coupled neural networks (PCNN)  or a spectral segmentation. Contour extraction of the segmented blobs will be performed with a newly developed algorithm, which is highly parallel and therefore fast. The contours and its points form the input to the Chordal Axis Transform  (CAT), which provides a morphological decomposition of grain shapes into simplicial chain complexes of limbs and torsos. This decomposition is used to perform geometric filtering operations such as grain separation based on grain morphology. Subsequent to such filtering, we explain how the metrical and statistical properties of the granular material can be efficiently computed to obtain input for the hydrodynamic model calculations.

## 2 IMAGE PROCESSING

In this paper, we discuss the image analysis of a micrograph (cf. Fig. 1.a) of an explosive material, PBX 9501 (95% HMX and 5% polymeric binder, by weight). More images similar to the one shown in Fig. 1.a can be found in the paper by Skidmore  et al.

### 2.1 Image Segmentation and Image Smoothing

An image can in general be viewed as the spacial arrangement (positions) of pixels, which contain spectral information (such as grey-levels, etc.). In Fig 1.b we show the grey-level histogram of the explosive material, which is shown in Fig 1.a. We can identify the pixels, which have rather low grey-levels (dark) in the range of 0–100, with the polymeric binder, and the remaining pixels (light) with the HMX crystals (grains). If we assign to the binder pixels black pixels and to the grain pixels white pixels, respectively, then this spectral segmentation leads to the bi-level image shown in Fig. 3.a.

A more sophisticated way of performing image segmentation can be acchieved by applying a PCNN . PCNN’s use, simultaneously, the spacial and spectral information content of an image for its segmentation into blobs. We have applied such a PCNN to Fig 1.a. The result is shown in Fig. 3.b (for more detail, cf. the paper by Skourikhine  et al.). In addition to segmentation, PCNN’s can also help to smooth an image, in order to suppress noise generated by the imager or to remove undesired texture. We have smoothed Fig. 1.a while using a PCNN. The resulting image

and its corresponding pixel grey-level histogram are displayed in Fig. 2. We would like to emphasize, that the application of a PCNN leads to smoothing without blurring the image . In Fig. 4. we have applied spectral and PCNN segmentation, respectively, as described above to the smoothed image Fig. 2.a. In total we have generated four different segmented versions of the original image, where apparantly the level of pixel noise decreases when going from Fig. 3.b over Fig. 3.a and Fig. 4.a to Fig. 4.b. In the following, we shall use the four segmented images for their further morphological analysis. This section also demonstrates, that image analysis is not purely objective, as we indicated in the introduction.

### 2.2 Contour Extraction

The next step after the image segmentation is the contour extraction of the blobs. For the material science application under consideration we desire blob contours, which are nondegenerate, i.e., they always enclose an area larger than zero, and they never cross or overlap each other. Furthermore, the contours should be oriented with respect to the objects and their possible holes.

In Fig. 5.a we show for the sake of illustration a cropped subimage of Fig. 1.a with increased magnification. In Fig. 5.b we have applied a spectral segmentation of the pixels grey-levels as described in the previous section. Finally, Fig. 5.c shows the dilated contours as generated by the algorithm , which is described in the paper by Schlei et al.  and which fulfills all of the above stated requirements for the contours.

### 2.3 Constrained Delaunay Tesselation and Chordal Axis Transform

The Constrained Delaunay Tesselation (CDT) of a simple planar polygon (contour) is a decomposition of a polygon into triangles, such that the circumcircle of each triangle contains no vertex of the polygon inside it that is simultaneously visible to two vertices of the triangle . Note, that the CDT of the contours is the key step that allows for the shape feature extraction. In Fig. 6.a we transform through the CDT the spectral pixel information of our given micrograph into an affine geometric description. As a next step we use the pointset of the dilated oriented contours and the generated triangular mesh with its neighbor information of the triangles to compute the CAT skeleton  of the grains (cf. Fig. 6.b). Each arc of the skeleton represents a simplical chain complex of either a limb or a torso, respectively. A limb is a chain complex of pairwise adjacent triangles, which begins with a junction triangle and ends with a termination triangle. A torso is a chain complex of pairwise adjacent triangles, which begins and ends with a junction triangle.

### 2.4 Grain Separation

In the micrograph (cf. Fig. 1.a) many grains might be very close to each other spatially such that they are connected through a single CAT skeleton. This is not satisfactory, if one wants to use the CAT skeleton directly for the computation of material statistics, because the number of grains is unequal to the number of skeletons. If we compute for each CDT triangle the average grey level of the pixels which are covered by it, we obtain Fig. 6.b. For each given torso of the CAT we can investigate the fluctuation of grey levels for its simplical chain complex. If two grains are connected in the CAT skeleton through a torso, we would expect a larger grey level fluctuation compared to the case, where a torso would be located within a single grain. In the following, we have used a threshold of 30

(among 256 grey levels) for cutting a torso, thus leading to grain separation. In Fig. 6.c we show the resulting refined grain contours. We note, that it is neccessary for the shape manipulation to use the unpruned CAT skeleton.

### 2.5 Material Statistics

After the grain separation it is in general neccessary to generate a new CDT of the refined contours (cf. Fig. 7.a). In the subsequent CAT our computer code  generates a connectivity hierarchy of the triangular chain complexes, which allows the efficient pruning of morphologically insignificant shape features. The results are the pruned skeleton for grains (cf. Fig. 7.b) and the pruned skeleton for the polymeric binder (cf. Fig. 7.c). We note, that the CAT skeleton of a grain can be sometimes pruned to a single point. After pruning, each connected arc of the attributed CAT skeletons can be regarded as the axis of a single particle (grain or binder particle). This allows for the computation of area, length, width, location, orientation, etc. for each particle. In case of present holes in a single HMX grain we even are able to remove them through the knowledge of the morphology of a single grains CAT skeleton. Of course, this requires a third application of CDT and CAT, respectively.

Figs. 8. and 9. show the processed images for Figs. 3 and 4. E.g., we find the following grain numbers when going from Fig. 8.a over Fig. 8.b and Fig. 9.a to Fig. 9.b: 912, 919, 624, and 534. Correspondingly, the two-dimensional HMX crystal percentages are: 75%, 72%, 74%, and 71%. Obviously, although the ratio between grains and polymer binder are more or less insensitive to the choice of image segmentation, there are apparently significant differences in

the grain numbers. Further experimental analysis should help to determine, which image segmentation is the most favorable for the computer vision micrograph analysis of explosive materials such as the PBX 9501 considered here.

## 3 SUMMARY

In summary, we have presented a method of image decomposition for micrographs into meaningful parts such as single HMX and binder grains. We have demonstrated that our algorithms are capable of grain separation and the efficient computation of statistical features of granular materials.

## ACKNOWLEDGMENTS

This work has been supported by the Department of Energy.