text:
# Coupled Two-Way Clustering Analysis of Gene Microarray Data

## Introduction

In a typical DNA microarray experiment expression levels of thousands of genes are recorded over a few tens of different samples111 By “sample” we refer to any kind of living matter that is being tested, e.g. different tissues cell populations collected at different times etc.  . Hence this new technology gave rise to a new computational challenge: to make sense of such massive expression data . The sizes of the datasets and their complexity call for multi-variant clustering techniques , which are essential for extracting correlated patterns and the natural classes present in a set of $`N`$ data points, or objects, represented as points in the multidimensional space defined by $`D`$ measured features.

Gene microarray data are fairly special in that it makes good sense to perform clustering analysis in two ways . The first views the $`n_s`$ samples as the $`N=n_s`$ objects to be clustered, with the $`n_g`$ genes’ levels of expression in a particular sample playing the role of the features, representing that sample as a point in a $`D=n_g`$ dimensional space. The different phases of a cellular process emerge from grouping together samples with similar or related expression profiles. The other, not less natural way, looks for clusters of genes that act correlatively on the different samples. This view considers the $`N=n_g`$ genes as the objects to be clustered, each represented by its expression profile, as measured over all the samples, as a point in a $`D=n_s`$ dimensional space.

Whereas in previous work  the samples and genes were clustered completely independently, we introduce and perform here a coupled two-way clustering (CTWC) analysis.

Our philosophy is to narrow down both the features that we use and the data points that are clustered. We believe that only a small subset of the genes participate in any cellular process of interest, which takes place only in a subset of the samples; by focusing on small subsets, we lower the noise induced by the other samples and genes. We look for pairs of a relatively small subset $``$ of features (either genes or samples) and of objects $`𝒪`$, (samples or genes), such that when the set $`𝒪`$ is clustered using the features $``$, stable and significant partitions are obtained. Finding such pairs of subsets is a rather complex mathematical problem; the CTWC method produces such pairs in an iterative clustering process.

CTWC can be performed with any clustering algorithm. We tested it in conjunction with several clustering methods, but present here only results that were obtained using the super-paramagnetic clustering algorithm (SPC) , which is especially suitable for gene microarray data analysis due to its robustness against noise and its “natural” ability to identify stable clusters.

The CTWC clustering scheme was applied to two gene microarray data sets, one from a colon cancer experiment  and the other from a leukemia experiment . From both datasets we were able to “mine” new partitions and correlations that have not been obtained in an unsupervised fashion by previously used methods. Some of these new partitions have clear, well understood biological interpretation. We do not report here discoveries of biologically relevant, previously unknown results. The main point of our message is twofold: (a) we were able to identify biologically relevant partitions in an unsupervised way and (b) other, not less natural new partitions were also found, which may contain new, important information and for which one should seek biological interpretation.

## Coupled Two Way Clustering

### Motivation and Algorithm

The results of every gene microarray experiment can be summarized as a set of numbers, which we organize in an expression level matrix $`𝒜`$. A row of this matrix corresponds to a single gene, while each column represents a particular sample. Our normalization is described in detail later.

In a typical experiment simultaneous expression levels of thousands of genes are measured. Gene expression is influenced by the cell type, cell phase, external signals and more . The expression level matrix is therefore the result of all these processes mixed together. Our goal is to separate and identify these processes and to extract as much information as possible about them. The main point is that each biological process on which we wish to focus may involve a relatively small subset of the genes that are present on a microarray; the large majority of the genes constitute a noisy background which may mask the effect of the small subset. The same may happen with respect to samples.

The CTWC procedure which we now describe is designed to identify subsets of genes and samples, such that a single process is the main contributor to the expression of the gene subset over the sample subset. We start with clustering the samples and the genes of the full data set and identify all stable clusters of either samples or genes. We scan these clusters one by one. The expression levels of the genes of each cluster are used as the feature set $``$ to represent object sets. The different object sets $`𝒪`$ contain either all the samples or any sample cluster. Similarly, we scan all stable clusters of samples and use them as the feature set $``$ to identify stable clusters of genes. We keep track of all the stable clusters that are generated, of both genes, denoted as $`v^g`$, and samples $`v^s`$. The gene clusters are accumulated in a list $`V^g`$ and the sample clusters in $`V^s`$. Furthermore, we keep all the chain of clustering analyses that has been performed (which subset was used as objects, which subset was used as features, and which were the stable clusters that have been identified).

When new clusters are found, we use them in the next iteration. At each iteration step we cluster a subset of the objects (either samples or genes) using a subset of the features (genes or samples). The procedure stops when no new relevant information is generated. The outcome of the CTWC algorithm are the final sets $`V^g`$ and $`V^s`$ and the pointers that identify how all stable clusters of genes and samples were generated.

A precise, step by step definition of the algorithm is given in Fig. 1.

### Analyzing the clusters obtained by CTWC

The output of CTWC has two important components. First, it provides a broad list of gene and sample clusters. Second, for each cluster (of samples, say) we know which subset (of samples) was clustered to find it, and which were the features (genes) used to represent it. We also know for every cluster $`𝒞`$, which other clusters can be identified by using $`𝒞`$ as the feature set. We present here a brief selection of the possible ways one can utilize this kind of information. Implementations of the particular uses listed here are described in the Applications section.

Identifying genes that partition the samples according to a known classification.  This particular application is supervised. Denote by $`C`$ a known classification of the samples, say into two classes, $`c_1`$ and $`c_2`$. CTWC provides an easy way to rank the clusters of genes in $`V^g`$ by their ability to separate the samples according to $`C`$. It should be noted that CTWC not only provides a list of candidate gene clusters one should check, but also a unique method of testing them.

First we evaluate for each cluster of samples $`v^s`$ in $`V^s`$ two scores, purity and efficiency, which reflect the extent to which assignment of the samples to $`v^s`$ corresponds to the classification $`C`$. These figures of merit are defined (for $`c_1`$, say) as

$$\text{purity}(v^s|c_1)=\frac{|v^sc_1|}{|v^s|};\text{efficiency}(v^s|c_1)=\frac{|v^sc_1|}{|c_1|}.$$

Once a cluster $`v^s`$ with high purity and efficiency has been found, we can use the saved pointers to read off the cluster (or clusters) of genes that were used as the feature set to yield $`v^s`$ in our clustering procedure. Clustering, as opposed to classification, discovers only those partitions of the data which are, in some sense, “natural”. Hence by this method we identify the most natural group of genes that can be used to induce a desired classification.

Needless to say, one can also test a gene cluster $`v^g`$ that was provided by CTWC using more standard statistics, such as the t-test  or the Jensen-Shannon distance . Both compare the expression levels of the genes of $`v^g`$ on the two groups of samples, $`c_1,c_2`$, partitioned according to $`C`$. Alternatively, one can also use the genes of $`v^g`$ to train a classifier to separate the samples according to $`C`$ , and use the success of the classifier to measure whether the expression levels of the genes in $`v^g`$ do or do not correspond to the classification.

Discovering new partitions. Every cluster $`v^s`$ of $`V^s`$ is a subset of all the samples, the members of which have been linked to each other and separated from the other samples on the basis of the expression levels of some co-expressed subset of genes. It is reasonable therefore to argue that the cluster $`v^s`$ has been formed for some biological or experimental reason.

As a first step to understand the reason for the formation of a robust cluster $`v^s`$, one should try to relate it to some previously known classification (for example, in terms of purity and efficiency). Clusters which cannot be associated with any known classification, have to be inspected more carefully. Useful hints for the meaning of such a cluster of samples may come from the identity of the cluster of genes which was used to find it. Clearly, the CTWC clusters can be used in the same way to interpret clusters of genes which were not previously known to belong to the same process.

CTWC is a sensitive tool to identify sub-partitions. Some of the sample clusters in $`V^s`$ may have have emerged from clustering a subset of the samples, say $`v_0^s`$. These clusters reflect a sub-partition of the samples which belong to $`v_0^s`$. When trying to cluster the full sample set, this sub-partition may be missed, since other samples, unrelated to $`v_0^s`$, are masking it.

CTWC reveals conditional correlations among genes. The CTWC method collects stable gene clusters in $`V^g`$. In many cases the same groups of genes may be added to $`V^g`$ more than once. This is caused by the fact that some genes are co-regulated in all cells, and therefore are clustered together, no matter which subset of the samples is used as the feature set. For example, ribosomal proteins are expected to be clustered together for any set of samples which is not unreasonably small.

Some gene clusters, however, are different; they are co-regulated only in a specific subset of samples. We call this situation conditional correlation. The identity of the sample cluster which reveals the conditionally correlated gene cluster is clearly important to understand the biological process which makes these genes correlated.

## Clustering method and similarity measures

Any reasonable choice of clustering method and definition of stable clusters can be used within the framework of CTWC. We describe here the benefits of the particular clustering algorithm and similarity measure we used, which we found to be particularly suitable to handle the special properties of gene microarray data.

### SPC provides clear identification of stable clusters in a robust manner.

Super-paramagnetic clustering (SPC) is a hierarchical clustering method recently introduced by Blatt et al . The intuition that led to it is based on an analogy to the physics of inhomogeneous ferromagnets. Full details of the algorithm and the underlying philosophy are given elsewhere .

As for many hierarchical clustering algorithms, the input for SPC is a distance or similarity matrix $`d_{ij}`$ between the objects $`𝒪`$, calculated according to the feature set $``$. A tunable parameter $`T`$ (’temperature’) controls the resolution of the performed clustering. One starts at $`T=0`$, with a single cluster that contains all the objects. As $`T`$ increases, phase transitions take place, and this cluster breaks into several sub-clusters which reflect the structure of the data. Clusters keep breaking up as $`T`$ is further increased, until at high enough values of $`T`$ each object forms its own cluster.

Blatt et al showed that the SPC algorithm is robust against variation of its parameters, initialization and against noise in the data. The following advantages of SPC makes it especially suitable for gene microarray data analysis: (i) No prior knowledge of the structure of the data is assumed; (ii) SPC provides information about the different self organizing regimes of the data; (iii) The number of “macroscopic” clusters is an output of the algorithm; and (iv) Hierarchical organization of the data is reflected in the manner clusters merge or split when the control parameter (the ’temperature’ $`T`$) is varied.

Moreover, the control parameter can be used to provide a natural measure for the stability of any particular cluster by the range of temperatures $`\mathrm{\Delta }T`$ at which the cluster remains unchanged. A stable cluster is expected to ’survive’ throughout a large $`\mathrm{\Delta }T`$, one which constitutes a significant fraction of the range it takes the data to break into single point clusters. Inspection of the gene dendrograms of Fig. 4 reveals stable clusters and stable branches.

### Normalization of the gene expression array

The Pearson correlation is commonly used as the similarity measure between genes or samples . This measure conforms with the intuitive biological notion of what it means for two genes to be co-expressed; this statistic captures similarity of the “shapes” of two expression profiles, and ignores differences between the magnitudes of the two series of measurements . The correlation coefficient is high between two genes that are affected by the same process, even if each has a different gain due to the process, over different background expression levels (caused by other processes). One problem of using the correlation coefficient is that its reliability depends on the absolute expression level of the compared genes; a positive correlation between two highly expressed genes is much more significant than the same value between two poorly expressed genes. This information is ignored in the clustering process.

However, we find that correlations do not always capture similarity between samples. For example, consider two samples taken at different stages of some process, with the expression levels of a family of genes much below average in one sample and much higher in the other. Even if the expression levels of the two samples over these genes are correlated, one would like to assign them into different clusters. Furthermore, the distance between the two samples should be affected by the statistical significance of their expression differences.

We therefore used the following normalization scheme. Denote by $`𝒟`$ the matrix of the raw data. $`𝒟`$ is a $`n_g\times n_s`$ matrix, where $`n_g`$ is the number of genes and $`n_s`$ the number of samples.

We normalize our expression level matrix in two steps. First, divide each column by its mean: $`𝒟_{}^{}{}_{ij}{}^{}=𝒟_{ij}/\overline{𝒟}_j`$ ; $`\overline{𝒟}_j=\frac{1}{n_g}_{i=1}^{n_g}𝒟_{ij}`$. We then normalize each row, such that its mean vanishes and its norm is one:

$$𝒜_{ij}=\frac{𝒟_{}^{}{}_{ij}{}^{}\overline{𝒟^{}}_i}{𝒟_{}^{}{}_{i}{}^{}},$$

where $`\overline{𝒟^{}}_i=\frac{1}{n_s}_{j=1}^{n_s}𝒟_{}^{}{}_{ij}{}^{}`$ and $`𝒟_{}^{}{}_{i}{}^{}^2=_{j=1}^{n_s}\left(𝒟_{}^{}{}_{ij}{}^{}\overline{𝒟^{}}_i\right)^2`$.

For genes and samples we use the Euclidean distance as the dissimilarity measure. For two genes (rows of $`𝒜`$) the Euclidean distance is closely related to the Pearson correlation between them.

## Applications

In order to show the strength of the CTWC algorithm, we apply it to two gene microarray experiment data sets. Here we report only the results which were obtained by CTWC, and could not be found using a straightforward clustering analysis. We highlight a small subset of the partitions that our method was able to extract from the data; these are the results for which we were able to find satisfactory biological explanation. We do not report here new discoveries of biologically relevant, previously unknown results. Rather, we claim to have discovered a method that is capable to mine such information out of the available data. New, relevant information may be contained in the new partitions which were found, to which we were not yet able to assign biological meaning. These new, uninterpreted results are reviewed briefly below; full lists of the clusters associated with these results, as well as their constituent samples or genes can be found at http://www.weizmann.ac.il/physics/complex/compphys.

### Analysis of Leukemia samples

We analyzed data obtained by Golub et al  from 72 samples collected from acute leukemia patients at the time of diagnosis. 47 cases were diagnosed as ALL (acute lymphoblastic leukemia) and the other 25 as AML (acute myeloid leukemia). RNA prepared from the bone marrow mononuclear cells was hibridized to high-density oligonucleotide micorarrays, produced by Affymetrix, containing 6817 human genes.

After rescaling the data in the manner described by Golub et al , we selected only those genes whose minimal expression over all samples is greater than 20. As a result of this thresholding operation 1753 genes were left. The resulting array was then normalized as described previously, to give the $`1753\times 72`$ expression level matrix $`𝒜`$ (see Fig. 2).

We found that two iterations of the CTWC algorithm sufficed to converge to 49 stable gene clusters (LG1-49) and 35 stable sample clusters (LS1-35). We highlight here four of our findings, which demonstrate the power of the method to solve problems listed above.

Identifying genes that partition the samples according to a known classification.  First we use the known ALL/AML classification of the samples to determine which gene clusters can distinguish between the two classes. We found only a single gene cluster (LG1) which enables stable separation into AML/ALL clusters222A cluster is identified with a certain class if both its purity and efficiency exceeds 3/4.. This well demonstrates the strength of CTWC, since it turned out that SPC was not able to clearly identify the AML/ALL separation using the full set of genes.

Discovering new partitions. Next, we search the stable sample clusters for unknown partitions of the samples. We focus our attention on sample clusters which were repeatedly found to be stable. One such cluster, denoted LS1, may be of interest; it includes 37 samples and was found to be stable when either a cluster of 27 genes (LG2) or another unrelated cluster of 36 genes (LG3) was used to provide the features. LG3 includes many genes that participate in the glycolysis pathway. Due to lack of additional information about the patients we cannot determine the biological origin of the formation of this sample cluster.

Identifying sub-partitions Using a 28 gene cluster (LG4) as features, we tried to cluster only the samples that were identified as AML patients (leaving out ALL samples). A stable cluster, LS2, of 16 samples was found (see Fig. 2(B)); it contains most of the samples (14/15) that were taken from patients that underwent treatment and whose treatment results were known (either success or failure). For none of the other AML patients was any information about treatment available in the data. Some of the 16 genes of this cluster, LG4, are ribosomal proteins and some others are related to cell growth. Apparently these genes can partition the AML patients according to whether they did or did not undergo treatment.

This result demonstrates a possible diagnostic use of the CTWC approach; one can identify different responses to treatment, and the groups of genes to be used as the appropriate probe.

We repeated the same procedure, but discarding AML and keeping only the ALL samples. We discovered that when any one of 5 different gene clusters (LG4-8) are used to provide the features, the ALL samples break into two stable clusters; LS5, which consists mostly of T-Cell ALL patients and LS4, that contains mostly B-Cell ALL patients (see Fig. 2(A)). When all the genes were used to cluster all samples, no such clear separation into T-ALL vs B-ALL was observed. One of the gene clusters used, LG5, with T/B separating ability, contains 29 genes, many of which are T-cell related. Another gene cluster, LG6, which also gave rise to T/B differentiation contains many HLA histocompatability genes.

These results demonstrate how CTWC can be used to characterize different types of cancer. Imagine that the nature of the sub-classification of ALL had not been known. On the basis of our results we could predict that there are two distinct sub-classes of ALL; moreover, by the fact that many genes which induce separation into these sub-classes are either T-Cell related or HLA genes, one could suspect that these sub-classes were immunology related.

As a different possible use of our results, note that some of the genes in the T-Cell related gene cluster LG5 have no determined function, and may be candidates for new T-Cell genes. This assumption is supported both by the fact that these genes were found to be correlated with other T-Cell genes, and by the fact that they support the differentiation between T-ALL and B-ALL.

### Analysis of Colon cancer data

The data set we consider next contains 40 colon tumor samples and 22 normal colon samples, analyzed with an Affymetrix oligonucleotide array complementary to more than 6500 human genes and ESTs. Following Alon et al , we chose to work only with the 2000 genes of greatest minimal expression over the samples. We normalized the data to get a $`2000\times 62`$ expression level matrix $`𝒜`$.

The CTWC algorithm was applied to this data set. 97 stable gene clusters (CG1-97) and 76 stable sample clusters (CS1-76) were obtained in two iterations.

Identifying genes that partition the samples according to a known classification.  Again we search first for gene clusters which differentiate the samples according to the known normal/tumor classification. We found 4 gene clusters (CG1-4) that partition the samples this way. The genes of these clusters can be used if one wishes to construct a classifier for diagnosis purposes (see Fig. 3(A)).

Discovering new partitions. Five clusters of genes (CG2,CG4-CG7) generated very stable clusters of samples. Two of the five (CG2,CG4) differentiated tumor and normal; two other were less interesting since the clusters they generated contained most of the samples. The gene cluster CG5, however, gave rise to a clear partition of the samples into two clusters, of 39 and 23 tissues (see Fig. 3(B)). Checking with the experimentalists333U.Alon, K.Gish, D.Mack & A.Levine, Private communication. We discovered that this separation coincides almost precisely with a change of the experimental protocol; 22 RNA samples were extracted using a poly-A detector (’protocol-A’), and the other 40 samples were prepared by extracting total RNA from the cells (’protocol-B’). Cursory examination did not yield any obvious common features among the 29 genes of the cluster CG5 that gave rise to this separation of the tissues.

Identifying conditionally correlated genes and sub-partitions Finally, we turn to identify conditionally correlated genes by comparing stable gene clusters formed when using different sample sets as features. We found that most gene clusters form irrespectively of the samples that are used. We did find, however, 4 special groups of genes (CG8-11) that formed clear and stable clusters when using only the tumor samples as features, but were relatively uncorrelated, i.e. spread across the dendrogram of genes, when clustering was performed based on all the samples or only the normal ones.

One of these 4 clusters, (CG9), breaks up, at a higher resolution, into two sub-clusters, as shown in Fig. 4(B). One of these sub-clusters, (CG12), consists of 51 genes, all of which are related to cell growth (ribosomal proteins and elongation factors). The other sub-cluster, (CG13), contains 17 genes, many of which are related to intestinal epithelial cells (e.g. mucin, cathespin proteases). Interestingly, when clustering the genes on the basis of the either all samples or only the normal ones, both clusters (CG12 and CG13) appear as two uncorrelated distinct clusters, and their positions in the dendrogram are quite far from each other (Fig. 4).

The high correlation between growth genes and epithelial genes, observed in tumor tissue, suggests that it is the epithelial cells that are rapidly growing. In the normal samples there is smaller correlation, indicating that the expression of growth genes is not especially high in the normal epithelial cells. These results are consistent with the epithelial origin of colon tumor.

Two other groups of genes formed clusters only over the tumor cells. One (CG11, of 34 genes) is related to the immune system (HLA genes and immunoglobulin receptors). The second (CG10, of 62 genes) seems to be a concatenation of genes related to epithelial cells (endothelial growth factor and retinoic acid), and of muscle and nerve related genes. We could not find any common function for the genes in the fourth cluster (CG8).

Clustering the genes on the basis of their expression over only the normal samples revealed three gene clusters (CG14-16) which did not form when either the entire set of samples or the tumor tissues were used. Again, we could not find a clear common function for these genes. Each cluster contains genes that apparently take part in some process that takes place in normal cells, but is suppressed in tumor tissues.

## Summary and discussion

We proposed a new method for analysis of gene microarray data. The main underlying idea of our method is to zero in on small subsets of the massive expression patterns obtained from thousands of genes for a large number of samples. A cellular process of interest may involve a relatively small subset of the genes in the dataset, and the process may take place only in a small number of samples. Hence when the full data set is analyzed, the “signal” of this process may be completely overwhelmed by the “noise” generated by the vast majority of unrelated data.

We are looking for a relatively small group of genes, which can be used as the features used to cluster a subset of the samples. Alternatively, we try to identify a subset of the samples that can be used in a similar way to identify genes with correlated expression levels. Identifying pairs of subsets of genes and samples, which produce significant stable clusters in this way, is a computationally complex task. We demonstrated that the Coupled Two-Way Clustering technique provides an efficient method to produce such subgroups.

The CTWC algorithm provides a broad list of stable gene and sample clusters, together with various connections among them. This information can be used to perform the most important tasks in microarray data analysis, such as identification of cellular processes and the conditions for their activation; establishing connection between gene groups and biological processes; and finding partitions of known classes of samples into sub-groups.

We reemphasize that CTWC is applicable with any reasonable choice of clustering algorithm, as long as it is capable of identifying stable clusters. In this work we reported results obtained using the super-paramagnetic clustering algorithm (SPC), which is especially suitable for gene microarray data analysis due to its robustness against noise which is inherent in such experiments.

The power of the CTWC method was demonstrated on data obtained in two gene microarray experiments. In the first experiment the gene expression profile in bone marrow and peripheral blood cells of 72 leukemia patients was measured using gene microarray technology. Our main results for this data were the following: (i) The connection between T-Cell related genes and the sub-classification of the ALL samples, into T and B-ALL, was revealed in an unsupervised fashion. (ii) We found a stable partition of the AML patients into two groups: those who were treated (with known results), and all others. This partition was revealed by a cluster of cell growth related genes. This observation may serve as a clue for a possible use of the CTWC method in understanding the effects of treatment.

The second experiment used gene microarray technology to probe the gene expression profile of 40 colon tumor samples and 22 normal colon tissues. Using CTWC we find a different, less obvious stable partition of the samples into two clusters. To find this partition, we had to use a subset of the genes. The new partition turned out to reflect two different experimental protocols. We deduce that the genes which gave rise to this partition of the samples are the ones which were sensitive to the change of protocol.

Another result that was obtained in an unsupervised manner using CTWC, is the connection between epithelial cells and the growth of cancer. When we looked at the expression profiles over only the tumor tissues, a cluster of cell growth genes was found to be highly correlated with epithelial genes. This correlation was absent when the normal tissues were used.

These novel features, discovered in data sets which were previously investigated by conventional clustering analysis, demonstrate the strength of CTWC. We find CTWC to be especially useful for gene microarray data analysis, but it may be a useful tool for investigating other kinds of data as well.

#### Acknowledgments

We thank N. Barkai for helpful discussions. The help provided by U.Alon in all stages of this work has been invaluable; he discussed with us his results at an early stage, provided us his data files and shared generously his understanding and insights. This research was partially supported by the Germany - Israel Science Foundation (GIF).


format: text

text:
# New Class of Eigenstates in Generic Hamiltonian Systems

## Abstract

In mixed systems, besides regular and chaotic states, there are states supported by the chaotic region mainly living in the vicinity of the hierarchy of regular islands. We show that the fraction of these hierarchical states scales as $`\mathrm{}^\alpha `$ and relate the exponent $`\alpha =11/\gamma `$ to the decay of the classical staying probability $`P(t)t^\gamma `$. This is numerically confirmed for the kicked rotor by studying the influence of hierarchical states on eigenfunction and level statistics.

Typical Hamiltonian systems are neither integrable nor ergodic  but have a mixed phase space, where regular and chaotic regions coexist. The regular regions are organized in a hierarchical way  (see, e.g., Fig. 2a) and chaotic dynamics is clearly distinct from the dynamics of fully chaotic systems. In particular, chaotic trajectories are trapped in the vicinity of the hierarchy of regular islands. The most prominent quantity reflecting this, is the probability $`P(t)`$ to be trapped longer than a time $`t`$, which decays as

$$P(t)t^\gamma ,\gamma >1,$$
(1)
in contrast to the typically exponential decay in fully chaotic systems. While the power-law decay is universal, the exponent $`\gamma `$ is system and parameter dependent. The origin of the algebraic decay are partial transport barriers , e.g., Cantori, leading to a hierarchical structure of the chaotic region . Quantum mechanically the classical algebraic decay of $`P(t)`$ is mimicked at most until Heisenberg time .

Even after two decades of studying quantum chaos, the search for quantum signatures of this universal power-law trapping is still in its infancy: In fact, only conductance fluctuations of open systems have been investigated so far. It was semiclassically derived that these fluctuations should have a fractal dimension $`D=2\gamma /2`$ , which was confirmed in gold nanowires , semiconductor nanostructures , and numerics . Quite recently, a second type of conductance fluctuations in mixed systems has been discovered numerically , namely isolated resonances. There the classical exponent $`\gamma `$ seems to appear in the scaling of the variance of conductance increments, surprisingly, on scales below the mean level spacing, what is not understood so far. Thus, even for this subject there is a lack of basic understanding.

In this paper we present consequences of the classical $`P(t)t^\gamma `$ in closed quantum systems, namely the appearance of a new class of eigenstates. Different from the well studied regular and chaotic states, they are supported by the chaotic region, but predominantly live in the vicinity of the regular islands . As for decreasing $`\mathrm{}`$ they move deeper into the hierarchical structure of the chaotic region, we like to call them hierarchical states. We show that the fraction $`f_{\mathrm{hier}}`$ of such states scales as

$$f_{\mathrm{hier}}\mathrm{}^{11/\gamma },$$
(2)
which is confirmed numerically for the kicked rotor by studying their influence on eigenfunction and level statistics. We deduce relation (2) by combining the finite resolution of quantum mechanics for a given $`\mathrm{}`$ with the simplest model  describing classical transport in a mixed phase space. We are encouraged by the success of this approach, and think that it will help in understanding the quantum signatures of a mixed phase space.

We motivate and numerically verify our results using the well known kicked rotor, which is a paradigm for a generic Hamiltonian system  and has a time evolution described by the map:

$`q_{n+1}`$ $`=`$ $`q_n+p_n`$ (3)
$`p_{n+1}`$ $`=`$ $`p_n+K\mathrm{sin}q_{n+1}.`$ (4)
Its quantum properties are determined by the time evolution operator for one period

$$𝒰=\mathrm{exp}\left(\frac{i}{\mathrm{}}K\mathrm{cos}q\right)\mathrm{exp}\left(\frac{i}{\mathrm{}}\frac{p^2}{2}\right).$$
(5)
Fig. 2a shows a typical classical phase space, which is governed by one big chaotic invariant set with fractional phase space volume $`\rho _{\mathrm{ch}}`$ and islands of regular motion of total fraction $`\rho _{\mathrm{reg}}`$  arranged in a hierarchical way. Quantum mechanically, however, we find three types of eigenstates of $`𝒰`$ (Fig. 1) : There are ’regular’ states living on KAM-tori of the regular islands and there are ’chaotic’ states extending across most of the chaotic region as first described by Percival . In addition, there are eigenstates supported by the chaotic region, that are of a third type. They predominantly live in the vicinity of the regular islands with only a small contribution in the main part of the chaotic sea. These states are separated from the main chaotic sea by partial barriers of the classical phase space, e.g., Cantori or stable and unstable manifolds . For decreasing $`\mathrm{}`$ these states live deeper in the hierarchy of classical phase space and thus we like to call them hierarchical states. Denoting the fraction of the three types of states by $`f_{\mathrm{reg}},f_{\mathrm{ch}}`$ and $`f_{\mathrm{hier}}`$, respectively, one has

$$f_{\mathrm{reg}}+f_{\mathrm{ch}}+f_{\mathrm{hier}}=1=\rho _{\mathrm{reg}}+\rho _{\mathrm{ch}}.$$
(6)
Obviously for any given $`\mathrm{}`$ not all regular islands are resolved and therefore $`f_{\mathrm{reg}}(\mathrm{})<\rho _{\mathrm{reg}}`$ holds. We determined the island volumes, revealing that the $`\mathrm{}`$-dependence of $`f_{\mathrm{reg}}`$ is negligible compared to the $`\mathrm{}`$-dependence of $`f_{\mathrm{hier}}`$ and $`f_{\mathrm{ch}}`$. We thus use $`f_{\mathrm{reg}}=\rho _{\mathrm{reg}}`$ in the following and Eq. (6) reduces to

$$f_{\mathrm{ch}}+f_{\mathrm{hier}}=\rho _{\mathrm{ch}}.$$
(7)
We want to focus on the fraction $`f_{\mathrm{hier}}`$ of hierarchical states, which we now determine from their influence on eigenfunction and level statistics.

As a first quantity we determine $`\eta `$, the density of an eigenstate averaged over a large area in the chaotic region, far away from the islands. As long as this averaging area is large compared to $`\mathrm{}`$ its specific choice and the used representation (e.g., Husimi) does not affect the following analysis. Normalization for $`\eta `$ is chosen such that, an eigenstate uniformly covering the entire phase space would give $`\eta =1`$. The distribution of $`\eta `$ (Fig. 1) shows a peak for small $`\eta `$ stemming from the regular states, intermediate values from hierarchical states and a Gaussian-like peak that we ascribe to the chaotic states. Assuming that the fraction $`f_{\mathrm{ch}}`$ of chaotic states extends on average uniformly (neglecting scar effects ) across a fraction $`f_{\mathrm{ch}}\rho _{\mathrm{ch}}`$ of phase space, we fit a Gaussian with area $`f_{\mathrm{ch}}`$ and mean $`1/f_{\mathrm{ch}}`$ to the latter part of the $`\eta `$-distribution. We find that the mean $`1/f_{\mathrm{ch}}`$ is clearly distinct from $`1/\rho _{\mathrm{ch}}`$ (Fig. 1), the value one would obtain if there were eigenstates uniformly covering the entire chaotic region. The fit determines the fraction $`f_{\mathrm{hier}}=\rho _{\mathrm{ch}}f_{\mathrm{ch}}`$ of hierarchical eigenstates, which is shown as a function of $`\mathrm{}`$ in Fig. 2b.

Secondly, we want to quantify the influence of hierarchical states on level statistics for mixed systems. To this end we extend the approach of Berry and Robnik , which assumes a random superposition of a regular (Poissonian) spectrum with measure $`\rho _{\mathrm{BR}}`$ and a chaotic (GOE) spectrum with measure $`1\rho _{\mathrm{BR}}`$ : We take into account the hierarchical states and their spectrum. As most of them couple weakly to regular as well as chaotic states, their spectrum is superimposed independently to the rest of the spectrum. Typically, $`f_{\mathrm{hier}}`$ is smaller than $`f_{\mathrm{ch}}+f_{\mathrm{reg}}`$ and therefore to first order these levels are randomly placed in the spectrum. Thus, effectively one can use the original Berry-Robnik approach, but one has to interpret the parameter $`\rho _{\mathrm{BR}}`$ as the sum of regular and hierarchical fraction of states, $`\rho _{\mathrm{BR}}=f_{\mathrm{reg}}+f_{\mathrm{hier}}`$. We determine $`\rho _{\mathrm{BR}}`$ by fitting the nearest-neighbor level-spacing distribution, using the cumulative distribution as well as the so-called U-function . The obtained values of $`f_{\mathrm{hier}}`$ from both fitting procedures agree within the error bars and their average as a function of $`\mathrm{}`$ is shown in Fig. 2b.

The two methods for determining $`f_{\mathrm{hier}}`$, eigenfunction and level statistics, rely on several assumptions, and one may not expect identical absolute values. Both methods, however, show a clear power-law decay according to relation (2) for the fraction of hierarchical eigenstates as can be seen in Fig. 2b for three values of the kicking strength $`K`$, where $`\gamma `$ is extracted from the decay of the corresponding $`P(t)`$ shown in Fig. 2c.

We now want to derive Eq. (2). Let us remind that the classical power-law trapping $`P(t)t^\gamma `$ originates from partial transport barriers, e.g., Cantori, arranged in a hierarchical way in the chaotic part of phase space around islands of regular motion . The simplest model  yielding $`P(t)t^\gamma `$ describes this hierarchy of the chaotic sea as a chain of volumes $`\mathrm{\Omega }_n\omega ^n`$ ($`\omega <1,n=0,1,\mathrm{}`$) connected by a flux $`\mathrm{\Phi }_{n,n+1}\phi ^n`$ ($`\phi <1`$) as shown in Fig. 3a. Then $`\gamma =1/(1\mathrm{ln}\omega /\mathrm{ln}\phi )`$ holds . Quantum mechanically two neighboring volumes $`n`$ and $`n+1`$ are strongly coupled, as long as $`\mathrm{\Phi }_{n,n+1}>\mathrm{}`$, while they are weakly coupled in the opposite case . This determines a critical flux $`\mathrm{\Phi }_{n^{},n^{}+1}=\mathrm{}`$ with $`n^{}\mathrm{ln}\mathrm{}/\mathrm{ln}\phi `$. The volumes $`\mathrm{\Omega }_n`$ with $`n<n^{}`$ correspond to the main part of the chaotic sea and support the chaotic eigenstates. Regions with $`n>n^{}`$ correspond to the hierarchical part of the chaotic sea supporting the hierarchical states. Summation of their volumes, finally, yields the fraction

$$f_{\mathrm{hier}}\underset{n>n^{}}{}\mathrm{\Omega }_n\omega ^n^{}\mathrm{}^{11/\gamma }.$$
(8)
of hierarchical states. From $`\gamma >1`$ it follows, that in the limit $`\mathrm{}0`$ this fraction tends to zero, while the total number of hierarchical states goes to infinity. In case of a system with $`d2`$ pairs of conjugate variables, one needs to replace $`\mathrm{}`$ by $`\mathrm{}^{d1}`$ in the above derivation.

In order to check relation (8), we generated a random matrix model for the classical chain following the approach of Ref. . For each ensemble of these random matrices we fixed $`\omega ,\phi `$ and therefore also $`\gamma `$ and then varied $`\mathrm{}`$. We repeated the above eigenfunction and level statistics and extracted $`f_{\mathrm{hier}}`$ ($`f_{\mathrm{reg}}=\rho _{\mathrm{reg}}=0`$ in this model). The results for various $`\gamma `$ are shown in Fig. 3b. The numerical data agree very well with the expected power laws over three orders of magnitude in $`\mathrm{}`$ and therefore confirm relations (8) and (2).

In conclusion, we present numerical evidence that besides chaotic and regular eigenfunctions there exists a third class, namely hierarchical eigenstates. Their fraction is determined from eigenfunction and level statistics and shows a power law as a function of $`\mathrm{}`$. We explain the origin of the power law and relate its exponent to the well known power-law trapping of mixed systems.

We like to thank S. Tomsovic and L. Kaplan for helpful discussions.


text:
# Magnetic stress as a driving force of structural distortions: the case of CrN

## Abstract

We show that the observed transition from rocksalt to orthorhombic Pnma symmetry in CrN can be understood in terms of stress anisotropy. Using local spin density functional theory, we find that the imbalance between stress stored in spin-paired and spin-unpaired Cr nearest neighbors causes the rocksalt structure to be unstable against distortions and justifies the observed antiferromagnetic ordering. This stress has a purely magnetic origin, and may be important in any system where the coupling between spin ordering and structure is strong.

The simultaneous occurence of structural distortion and change of magnetic ordering is observed in several magnetic compounds. In general, the interplay between structural and spin degrees of freedom depends on the detailed configuration of the electronic and phonon states, and a unique picture that elucidates the coupling mechanism is missing. A large class of systems where this interplay is strong are materials containing ions with degenerate (usually $`e_g`$) orbitals (Jahn-Teller ions), such as Mn3+, Cr2+ or Cu2+. The degeneracy causes instabilities that can be relieved by structural distortions or orbital (and eventually charge) ordering. In this paper we focus on CrN, one of the most widely employed materials in coating technology, due to the extreme hardness of CrN films. It is also the prototype material of strong magneto-structural interactions, since here, more than in any other case, the coupling is clearly revealed: at room temperature it is paramagnetic (PM) in the RS structure, like most of the transition metal nitrides, but at Néel temperature $`T_N`$=273-286 K, a both structural and magnetic transition to an antiferromagnetic (AFM) orthorhombic Pnma phase is observed.

But other compounds presents indications of strong spin and structural coupling. To name a few of them, in the perovskite LaMnO3 the distortions are necessary to stabilize the observed A-type antiferromagnetic ordering. The doped cobalt perovskite La1-xSrxCoO3 exhibits unusual properties as a function of doping and temperature, and in the diluted magnetic semiconductors Zn1-xCrxSe, the magnetic Cr2+ impurities are strongly affected by Jahn-Teller distortions. Also large is the class of materials manifesting magnetoelastic effects, such as changes of magnetic anisotropy coupled to distortions and strains. For example, the martensitic phase transition of Ni2MnGa from cubic Heusler to tetragonal structure is accompanied by a change in magnetic anisotropy and driven by magnetoelastic coupling.. In the layered Cu/Co/Ni/Cu/Si(001) the magnetic anisotropy can be induced by strain and controlled through the film thickness, and in CeNiSn the magnetic ordering can be changed by application of uniaxial strain.

For all of these materials, the stress is an important quantity, being naturally able to bridge electronic and structural properties. The stress relief argument is very popular in surface physics to explain processes like surface reconstructions, adsorbtions and growth: after a structural perturbation, a stress excess is produced in the system (e.g. the surface stress), that is relieved by further structural rearrangements involving symmetry breaking and eventually change of atomic density at surface. In this paper we show that is possible to apply similar arguments to magnetic systems. Specifically, we use the concept of stress relief to explain the stability of one magnetic ordering with respect to another, and to show how the change of ordering is related to a structural distortion.

In analogy to the surface stress, we define magnetic stress to be the stress associated with a change of spin ordering. If $`O^1`$ and $`O^2`$ represent two magnetic phases of a system with the same structure, the magnetic stress $`\widehat{T}`$ associated with the transition is:

$`\widehat{T}=\widehat{T}^{O^1}\widehat{T}^{O^2}={\displaystyle \frac{E^{O^1}}{\widehat{ϵ}}}{\displaystyle \frac{E^{O^2}}{\widehat{ϵ}}},`$ (1)
where $`\widehat{ϵ}`$ is the strain tensor. $`\widehat{T}`$ has a purely magnetic origin, since it originates from changes in spin orientations, not structural transitions, but may be relieved by structural transitions. We also define the magnetic stress per bond by mapping $`\widehat{T}`$ onto the Heisenberg model,

$`\widehat{T}={\displaystyle \underset{ij}{}}S_iS_j\widehat{t_{ij}},`$ (2)
where $`\widehat{t_{ij}}`$ is the strain derivative of the usual energy exchange interaction parameter $`J_{ij}`$, and represents the change of stress caused by turning the spin orientation on atoms $`i`$ and $`j`$ from parallel to antiparallel. A positive $`\widehat{t_{ij}}`$ means that the single spin-flip produces a tensile stress that can be relieved by shortening the $`i`$-$`j`$ distance. A negative $`\widehat{t_{ij}}`$, instead, represents a compressive stress that is relieved by stretching the $`i`$-$`j`$ distance. The meaning of the definitions of Eqs.1 and 2 will be transparent when applied to the practical case of CrN.

Below $`T_N`$, CrN is orthorhombic, with an unusual magnetic ordering (see Fig.1, left) built up by (1$`\overline{1}`$0) FM layers whose spin direction alternates up and down after each 2 layers, so that the system is antiferromagnetic (AFM) overall (in the following this phase is labeled AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$). Along with the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$, we consider two other AFM phases as possible competitors: the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$, similar to the former but with a one-by-one compensation of FM layers along \[1$`\overline{1}`$0\] (Fig.1, right), and the AFM$`{}_{[111]}{}^{}{}_{}{}^{1}`$, made by (111) FM layers alternated along , common for some strongly ionic transition-metal oxides (e.g. MnO).

Our calculations are based on local spin density functional theory (LSD), and employ plane-wave basis and Troullier-Martins pseudopotentials (this methodology already described successfully some properties of CrN ). The microscopic stress tensor is calculated in LSD (within the framework of the ABINIT project) and problems of finite-size basis have been efficiently circumvented.

To rationalize the RS-to-orthorhombic transition, we can picture it as a three step process: first, an hydrostatic expansion of the volume accompanies the non spin polarized-to-spin polarized transition. This is the dominant step in terms of energy stabilization. Second, a further stabilization is obtained from the competition of different magnetic orderings with the same (cubic) structure. Third, a shear strain is applied onto the (001) plane, in accordance with the experimental finding. The observed distortion consists of a reduction of the lattice parameter $`a`$ and an increase of $`b`$ (of $``$ 2%) with respect to their ideal values ($`a=a_0\sqrt{2}`$, $`b=b_0/\sqrt{2}`$, and $`c=a_0`$, where $`a_0`$ is the cubic lattice constant), in a way that the Cr-N distance is unchanged. Although a second-order effect with respect to the hydrostatic expansion, we will show that the distortion is fundamental in stabilizing the observed AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ phase.

The effects related to the first step are described in Fig.2, where the energies and stresses of PM and FM orderings in the RS structure are shown. In the PM phase, the high density of states (DOS) at the Fermi energy (EF) strongly drives CrN towards a Stoner instability. Although the high DOS peak at EF disappears in the FM phase, the spin splitting does not manage to open a gap, thus, FM CrN is still a (poor) metal. It is also a robust magnetic compound, with a magnetic moment (M=2.09 $`\mu _B`$ per Cr atom) almost entirely due to the Cr $`t_{2g}`$ splitting, and a chemical configuration close to Cr2+N2-.

We find a difference of 0.182 eV/formula unit between PM and FM energies, calculated at the equilibrium $`a_0`$ of the respective phases, that are 7.72 a.u. and 7.84 a.u. for PM and FM, respectively. Thus, after polarization, the lattice constant increases $``$ 1.5%. This strong magnetovolume effect is related to a large magnetic stress. At the PM lattice constant (where TPM is zero by definition), the spin-polarization produces a compressive stress TFM=1.89 eV/formula unit, whose relief corresponds to an energy gain $`\mathrm{\Delta }E^{FM}`$=TFM$`\mathrm{\Delta }`$ a/a$``$ 0.03 eV/formula unit.

The increase in compressive stress is a consequence of some charge localization (mostly Cr $`d`$ $`t_{2g}`$ charge) around the Cr nucleus, due to the spin polarization. The charge localization increases the kinetic stress, and, in turn, decreases the tensile electron-ion interaction.

Within this expanded volume, but still retaining the RS symmetry, we consider the competition among different magnetic orderings (the calculated lattice constant is quite similar for FM and AFM, i.e. a change of magnetic ordering does not produce that significant hydrostatic compression). We find (Table I) all the AFM phases are favored over the FM one. Moreover, the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$, and not the observed AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$, is the stable phase in the RS structure.

The density of states (DOS) of a single Cr ion in the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ phase is shown in Fig. 3. The Cr configuration is nearly Cr2+, with magnetization M=2.12 $`\mu _B`$. Notice that in the (001) plane of the orthorhombic cell, the $`t_{2g}`$ orbitals (whose DOS is drawn by the shaded area in the Figure), not the $`e_g`$s, are directed along $`\widehat{x}`$ and $`\widehat{y}`$. The $`e_g`$ orbitals are non degenerate, due to the AFM symmetry, and each of them is occupied by nearly half an electron. Their total contribution to the magnetization is only $``$ 0.2 $`\mu _B`$. The exchange splitting for $`t_{2g}`$ orbitals is $``$ 2 eV, with the $`t_{2g}^{}`$ almost completely filled. The global $`t_{2g}`$ occupation is 2.76, and the contribution to the magnetic moment is 1.9 $`\mu _B`$.

We can understand the competition among spin orderings by a two-parameter Heisenberg model, taking into account the nearest neighbor Cr-Cr ($`J_1`$), and the next nearest neighbor Cr-N-Cr interactions ($`J_2`$). The calculated parameters are $`J_1`$=-9.5 meV, and $`J_2`$=4 meV. For the Cr-N-Cr interaction there is competition between superexchange (AFM) and double-exchange (FM) mechanisms involving Cr $`e_g`$ and N $`p`$ orbitals ($`dp\sigma `$ hybridization). The resulting coupling is slightly FM (i.e. $`J_2`$ is positive). The direct Cr-Cr coupling, instead, is mediated by $`t_{2g}`$ orbital interactions ($`dd\sigma `$), and is AFM (i.e. $`J_1`$ is negative), as we may argue assuming, as usual, the Hund’s rule. The formal equivalence of Cr2+ and Mn3+ ions may suggest some similarities between CrN and LaMnO3. Indeed, in both of them the magnetic ion is surrounded by an octhaedral cage of cations. However, in LaMnO3 the Mn coupling via oxygens is the only relevant interaction, whereas, for CrN the Cr-N-Cr coupling is smaller than the direct Cr-Cr AFM coupling. As a consequence, the stability increases with the number of AFM nearest neighbors (8 for the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$, 6 for AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ and AFM). Also, due to the symmetry of the Cr fcc sublattice, no other arrangement with more than 8 AFM nearest neighbors is possible. In other words, the signs of $`J_1`$ and $`J_2`$ rule out the possibility that other orderings, not taken into account here, might be more stable.

On top of the RS structure, we finally apply a shear strain $`ϵ`$ onto the (001) plane of the competing phases AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$ and AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ (Table II). First consider the stress values for the undistorted (i.e. $`ϵ=0`$) structure: the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ phase is under a condition of appreciable stress with respect to the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$. This stress originates from the spin ordering asymmetry on the (001) plane: in the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ phase, each Cr has two spin-antiparallel nearest neighbors along $`\widehat{x}`$, and two spin-parallel along $`\widehat{y}`$. The bonds between AFM neighbors are under condition of tensile stress (i.e. $`T_{xx}>0`$), whereas the bonds between FM neighbors store compressive stress (i.e. $`T_{yy}<0`$). Clearly, a planar distortion will be able to relieve the stress and lower the energy of the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$. In contrast, the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$, which is isotropic on (001), cannot undergo the same stabilization: each Cr has two spin-parallel nearest neighbors in both $`\widehat{x}`$ and $`\widehat{y}`$ directions, and the eventual energy gain due to the stretch along $`\widehat{y}`$ would be lost due to the compression along $`\widehat{x}`$. Also, $`T_{zz}`$ is different for the two phases: in the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$, the larger number of FM Cr-Cr bonds having an orthogonal component causes an excess of compressive stress along $`\widehat{z}`$.

Applying Eq.2 we can express the stress (at $`ϵ`$=0) in terms of exchange-interaction parameters. For a planar strain on (001), we have three independent parameters corresponding to the in-plane Cr-Cr bonds ($`t_1`$), and the planar and orthogonal components of the out-of-plane Cr-Cr bonds ($`t_{2x}`$ and $`t_{2z}`$). Our results gives $`t_1`$=1.04 eV, $`t_{2x}`$= 0.24 eV, and $`t_{2z}`$=0.20 eV. Thus, the parallel-to-antiparallel spin flip causes a significantly large tensile stress for the Cr-Cr bonds parallel to (001).

On the basis of the magnetic stress results at $`ϵ`$=0 we are able not only to predict the occurence of the distortion for AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$, but also to estimate the energy gain due to the stress relief. Indeed, a 2% distortion gives an energy gain $`\mathrm{\Delta }E=\left(T_{xx}T_{yy}\right)\frac{\delta a}{a}\mathrm{\hspace{0.17em}0.04}eV`$ which is enough for the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ to be stabilized over the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$. This estimation can be compared with the direct calculation of energies and stresses at $`ϵ0`$. As expected, the distortion progressively relieves the stress and stabilizes the observed phase. Our results gives a theoretical distortion slightly larger than the experimental value: at $`ϵ=4\%`$ most of the planar stress is relieved. In contrast, the distortion on the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$ increases the stress and is energetically unfavorable.

The stabilization of AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ ordering is ultimately due to different stress stored in the bonds between spin-antiparallel and spin-parallel Cr nearest neighbors. Stress asymmetry and ordering of $`t_{2g}`$ orbitals are intimately related: since the Cr-Cr interaction occurs via direct $`t_{2g}`$ coupling and is AFM, $`J_1`$ should vanish in the limit of large nearest neighbor distance, whereas a shortening of the distance should increase the overlap between $`t_{2g}`$ orbitals, making $`J_1`$ even more negative.

This explains why a bond contraction stabilizes the spin-antiparallel interactions, while a bond stretching stabilizes the spin-parallel interactions. In Fig.4 we show the DOS of $`d_{x^2}`$ and $`d_{y^2}`$ orbitals (that contribute most to the Cr-Cr planar coupling) for the undistorted AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ ordering (obviously they are degenerate for AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{1}`$). The tendency of these orbitals to increase the overlap (and, thus, the hybridization) along $`\widehat{x}`$, and to stay more localized along $`\widehat{y}`$ is clearly visible.

In conclusion, in this paper we have given a definition of magnetic stress, and proposed it as a tool for understanding and predicting the subtle relations between magnetic ordering and structural properties. For CrN, we show that the coupling between the AFM$`{}_{[1\overline{1}0]}{}^{}{}_{}{}^{2}`$ ordering and the orthorhombic $`P_{nma}`$ structure is driven by the relief of tensile stress stored in the AFM Cr nearest neighbors, and associated with the direct interaction of $`t_{2g}`$ orbitals. Our results for CrN may stimulate the application of arguments based on stress relief to many other magnetic compounds. Indeed, the predictive power of the stress should be important in every case where magnetism and structural properties are strongly related.

We acknowledge financial support from the National Science Foundation POWRE program under grant number DMR 9973859. Many thanks to Xavier Gonze and the ABINIT group.


text:
# Addendum to “Multichannel Kondo screening in a one-dimensional correlated electron system”

## Abstract

This is an addendum to our previous work cond-mat/9705048 (published in Europhysics Letters 41, 213 (1998)), clarifying the construction of the two-particle scattering matrices used for studying the magnetic impurity behavior in a multichannel correlated host.
Published in Europhysics Letters 50, 125 (2000).

We point out that our lattice host Hamiltonian in Eq. (1), cond-mat/9705048, (Europhysics Letters 41, 213 (1998)) does not directly produce the electron-electron scattering matrix $`\widehat{X}(p_ip_j)`$ as given in the subsequent text. To obtain the factorizable form of this matrix one must first linearize the dispersion relation of the itinerant electrons about the Fermi levels (with the Fermi velocity normalized to unity). The linearized impurity Hamiltonian is constructed from the lattice form, Eq. (3), in a similar way (for details, see P. Schlottmann and A. A. Zvyagin, Phys. Rev. B, 55, 5027 (1997)). With a linearized spectrum all our results for the behavior of the magnetic impurity are ensured to be valid.


text:
# Charge-density-wave formation by Van Hove nesting in the 𝛼-phase of Sn/Ge(111)

\[

## Abstract

We study the role of electron correlations in the formation of the surface charge-density-wave state in the Sn/Ge(111) interface. The Fermi energy of the overlayer is treated as a dynamical variable, which undergoes a substantial renormalization by the interaction. We show that the Fermi level turns out to be pinned to a Van Hove singularity in the density of states, which explains the formation of the charge-density-wave, the observation of a very flat band in photoemission experiments and the reduction of the spectral weight in the low-temperature phase.

\]

The formation of states with rearrangement of the electronic charge is accompanied in general by interesting properties in condensed matter systems. A recent example has been provided by the study of interfaces like Pb/Ge(111), Sn/Ge(111) or K/Si(111). All of them have in common the existence of a phase transition with the formation of a surface charge-density-wave (CDW) in the low-temperature phase. In the first two cases the adsorbate phases have at room temperature the same $`(\sqrt{3}\times \sqrt{3})R30^{}`$ structure, and they both show a reduction of the symmetry to a $`(3\times 3)`$ ground state below the critical temperature. However, while the Pb/Ge(111) overlayer becomes insulating in the charge-ordered state, the Sn/Ge(111) remains metallic below the transition. In the first case, the analysis of the surface band has revealed nesting along significant portions of the two-dimensional Fermi surface, which would account for the formation of the CDW. In the $`\alpha `$-phase of Sn/Ge(111), however, no similar effect seems to be at work, and the reason why the ground state of the system is led to reorder its charge is not fully understood yet.

It is believed that correlation effects may play an important role in the development of the charge-ordered state. The relatively small surface bandwidth of these systems ($`0.5`$ eV) leads to consider the on-site Coulomb repulsion as one of the relevant interactions. It has been also suggested that a nearest-neighbor repulsion may play an important role in producing the $`3\times 3`$ surface periodicity.

In this letter we develop the idea that the electron correlations are essential to understand the low-temperature phase of Sn/Ge(111), while relying on the fact that no strong-coupling effects need to be invoked for that purpose. This is plausible as long as the overlayer remains metallic below the transition. We will show that the $`(\sqrt{3}\times \sqrt{3})`$ phase is unstable towards the formation of saddle points near the Fermi surface in the electronic dispersion relation. This is in agreement with the observation of a very flat band below the Fermi level in the angle-resolved photoemission experiments carried out in the CDW phase. The Van Hove singularity controls the dynamics of the system since, as we will see, this finds energetically favorable to readjust the bands in order to place the Fermi surface closer to the saddle points. As well as nesting of the saddle points near the Fermi surface gives rise to a spin-density-wave instability in the repulsive $`tt^{}`$ Hubbard model, it is able to explain the formation of the CDW under the assumption of a more general interaction in the present instance.

We will take an extended Hubbard model with on-site $`U`$ and nearest-neighbor $`V`$ repulsive interactions on the triangular lattice of the Sn/Ge(111) overlayer to illustrate these ideas. The adatom dangling bonds give rise to a band that reaches a minimum at the corner of the hexagonal $`(\sqrt{3}\times \sqrt{3})`$ Brillouin zone, and a maximum close to the center. We use a tight-binding approximation with nearest neighbor hopping $`t`$ and next-to-nearest neighbor hopping $`t^{}`$ to parametrize the Fermi line of the $`(\sqrt{3}\times \sqrt{3})`$ phase as given in Ref. . Although the precision of the fit to the dispersion relation is limited with this approximation, we are able to capture in this way the main topological features, which play an essential role in the following considerations.

The topology of the dispersion relation is actually unstable under the smallest periodic lattice distortion. This has to be taken into account, as it has been recently observed that one of the Sn atoms in the $`(3\times 3)`$ unit cell is slightly displaced upwards with respect to the other two. By taking a hopping parameter equal to $`0.9t`$ between inequivalent Sn atoms in the $`(3\times 3)`$ unit cell, we already observe the appearance of two subbands of very different character. The original valence band for the $`(\sqrt{3}\times \sqrt{3})`$ lattice splits along the $`\mathrm{\Gamma }`$K direction of the $`(3\times 3)`$ Brillouin zone. Thus, there are an upper and a lower subband that join at the apex of the conical dispersion shown in Fig. 1, characteristic of the underlying honeycomb lattice, and a middle subband that develops a saddle point at the boundary of the Brillouin zone, as shown in Fig. 2. The density of states redistributes itself, having a linear dependence in energy near K in the upper and lower honeycomb-type subbands while becoming divergent at the M point in the middle subband.

Before switching the interaction, the Fermi level has a nominal location below the mentioned Van Hove singularity. However, the relative position of the Fermi level in the different subbands of the model may suffer significant changes by the effect of the electron correlations near the singularity. The forward scattering vertex $`\mathrm{\Gamma }_{\mathrm{fwd}}`$ giving the four-point interaction at small momentum transfer $`𝐤`$ gets corrected by the iteration of the particle-hole polarizability. For the saddle-point dispersion relation $`\epsilon (𝐤)k_x^2k_y^2`$, this object diverges in the static limit as $`\mathrm{log}(\mathrm{\Lambda }/\mathrm{max}\{|\epsilon (𝐤)|,|\epsilon _F|\})`$, where $`\epsilon _F`$ is the distance of the Fermi level to the Van Hove singularity and $`\mathrm{\Lambda }`$ is a cutoff of the order of the bandwidth. In the case of a spin-independent interaction, this leads to a reduction of the vertex at zero momentum-transfer

$$\mathrm{\Gamma }_{\mathrm{fwd}}(\epsilon _F)=\frac{U_{\mathrm{fwd}}}{1+(c/(2\pi ^2\mathrm{\Lambda }))U_{\mathrm{fwd}}\mathrm{log}|\mathrm{\Lambda }/\epsilon _F|}$$
(1)
where the bare coupling is $`U_{\mathrm{fwd}}U(𝐩_1,𝐩_2;𝐩_1,𝐩_2)U(𝐩_1,𝐩_2;𝐩_2,𝐩_1)`$, in terms of the original four-fermion interaction, and $`c`$ is the prefactor of the particle-hole polarizability. This renormalization of $`\mathrm{\Gamma }_{\mathrm{fwd}}`$ has important consequences, since it encodes the sum of dominant logarithmic contributions that give rise to the constant shift of the chemical potential in the self-energy corrections. In this respect, it has been shown that the Van Hove singularity has a natural tendency to pin the Fermi level in a model with bare repulsive interactions. In the present case, this can be understood in terms of the uneven screening of the forward scattering in the different subbands of the model, as we show in what follows.

Let us model the problem of the interaction between the two subbands in Figs. 1 and 2 by taking the divergent density of states near the saddle-point structure

$$n^{(1)}(\epsilon )=\frac{1}{\mathrm{\Lambda }}\mathrm{log}(|\epsilon |/\mathrm{\Lambda }),\mathrm{\Lambda }<\epsilon <\mathrm{\Lambda }$$
(2)
and a vanishing density of states for the lower subband

$$n^{(2)}(\epsilon )=\frac{\alpha }{\mathrm{\Lambda }^2}|\epsilon |,\beta \mathrm{\Lambda }<\epsilon <\beta \mathrm{\Lambda }$$
(3)
The physical value of the Fermi energy comes from the balance between the nominal value $`\mu `$ of the chemical potential and the shift of the one-particle levels to higher energies due to the repulsive interaction. The filling level of each subband is determined by the self-energy corrections, bearing in mind that the upward displacement of the levels depends in turn on the charge present in all the subbands. The filling level $`\epsilon _{F1}`$ in the middle subband is given by the equation

$$\epsilon _{F1}=\mu \mathrm{\Gamma }_{\mathrm{fwd}}(\epsilon _{F1})^{\epsilon _{F1}}𝑑\epsilon n^{(1)}(\epsilon )g_{\mathrm{fwd}}^{\epsilon _{F2}}𝑑\epsilon n^{(2)}(\epsilon )$$
(4)
where we have introduced a coupling constant $`g_{\mathrm{fwd}}`$ that parametrizes the Coulomb repulsion exerted on the middle subband by the charge present in the lower subband. Analogously, the filling level $`\epsilon _{F2}`$ in the latter is computed in the form

$$\epsilon _{F2}=\mu U_{\mathrm{fwd}}^{\epsilon _{F2}}𝑑\epsilon n^{(2)}(\epsilon )g_{\mathrm{fwd}}^{\epsilon _{F1}}𝑑\epsilon n^{(1)}(\epsilon )$$
(5)
We remark that $`\epsilon _{F1}`$ and $`\epsilon _{F2}`$ are measured in the reference frames in which the dependences $`n^{(1)}(\epsilon )`$ and $`n^{(2)}(\epsilon )`$ are as given in Eqs. (2) and (3). Thus, the fact that $`\epsilon _{F1}`$ and $`\epsilon _{F2}`$ are nominally different after renormalization is just an artifact of that convention. The physical picture is however the opposite, namely that the one-particle levels are shifted to higher energy by a different amount in each of the subbands, up to a point in which the respective Fermi levels reach the common chemical potential.

The coupled set of equations (4) and (5) gives rise to nontrivial physical effects, as a consequence of the nonlinearities introduced by the divergent density of states $`n^{(1)}(\epsilon )`$ and the renormalization of $`\mathrm{\Gamma }_{\mathrm{fwd}}(\epsilon )`$ close to the Van Hove singularity. It is convenient to solve for the location of $`\epsilon _{F1}`$ and $`\epsilon _{F2}`$ in terms of the total charge $`N`$ in the system, given by

$$N=^{\epsilon _{F1}}𝑑\epsilon n^{(1)}(\epsilon )+^{\epsilon _{F2}}𝑑\epsilon n^{(2)}(\epsilon )$$
(6)
The most remarkable effect is that there is not a one-to-one correspondence between $`N`$ and the respective filling levels $`\epsilon _{F1}`$ and $`\epsilon _{F2}`$ for the two subbands. The different branches of the solution are represented in Fig. 3 for the particular values $`U_{\mathrm{fwd}}=4.0\mathrm{\Lambda }`$ and $`g_{\mathrm{fwd}}=3.0\mathrm{\Lambda }`$. The parameter $`\beta `$ has been chosen equal to $`3.0`$, and $`\alpha `$ has been set so that there is the same number of states in the lower and in the middle subband.

At low values of $`N`$, the filling of the middle subband with the Van Hove singularity proceeds in a regular way, with a monotonous increase of $`\epsilon _{F1}`$. There is a point, however, in which two other locations of $`\epsilon _{F1}`$ become possible, closer to the singularity in the density of states. In these instances, the corresponding filling level $`\epsilon _{F2}`$ in the second subband suffers an abrupt decrease with respect to the expected value. It is interesting to discern what of the possible solutions is most favorable energetically. We have plotted in Fig. 4 the value of the total energy $`E`$ versus the total charge $`N`$. We see that above a certain value, that is $`N1.4`$ in our case, any of the two filling levels close to the Van Hove singularity give the most stable configuration of the system. This is in agreement with previous analyses of the pinning of the Fermi level of electrons near a Van Hove singularity.

The result that turns out to be valid under very general conditions is the existence of an intermediate range of filling levels that are forbidden below the Van Hove singularity. For the corresponding values of the charge $`N`$, this finds more favorable to fill the Fermi sea up to the Van Hove singularity, at the expense of the charge in the other subband. In general, there is a critical value of $`N`$ in which the filling level $`\epsilon _{F1}`$ jumps discontinously from the regular evolution upon doping to a position much closer to the Van Hove singularity.

Coming back to the physics of the Sn/Ge(111) overlayer, a description consistent with the plots in Figs. 1 and 2 leads to a nominal value of $`\epsilon _{F1}`$ that is about $`25\%`$ off the Van Hove singularity, relative to the energy difference between the level of the latter and the bottom of the subband. This nominal value falls certainly in the forbidden region, what leads us to conclude according to the above arguments that the Fermi level is pinned close to the Van Hove singularity developed in the system.

This abrupt renormalization of the Fermi level towards the Van Hove singularity is able to account for all the prominent experimental features found in the Sn/Ge(111) interface, as we show in what follows.

1. Strong correlations in the low-temperature phase.

The different photoemission experiments have in common the observation in the low-temperature phase of a lower band that disperses along the $`\mathrm{\Gamma }`$M direction of the $`(3\times 3)`$ Brillouin zone, together with another band with higher energy and much less dispersion. The point is that this intermediate band is found at an energy sensibly lower than predicted by LDA calculations. Quite remarkably, all the points along the $`\mathrm{\Gamma }`$M direction are found below the Fermi level. After taking into account the experimental overestimation of the binding energies remarked in Ref. , the M point turns out to be close but below the Fermi level, contrary to the estimates carried out neglecting many-body effects. This is actually the experimental signature of the pinning to the Van Hove singularity we have discussed above.

It is also clear that the width of the middle subband is sensibly smaller than predicted by LDA calculations. This reduction in the bandwidth is also consistent with the renormalization of the kinetic energy of electrons close to a saddle point, which leads to a flatter shape of the dispersion relation.

2. Formation of the surface CDW.

The saddle points of the middle subband are at the M points of the $`(3\times 3)`$ Brillouin zone, so that they are separated by momenta $`𝐊_i`$ that correspond actually to the wavevectors of the CDW structure observed experimentally. This Van Hove nesting by momenta $`𝐊_i`$ connecting the saddle points near the Fermi level leads to the CDW instability.

In the context of an extended Hubbard model and with the inclusion of the electron-phonon interaction, it has been shown that it is natural to have a negative coupling $`U_{𝐊_i}`$ for the four-fermion interaction with momentum transfer $`𝐊_i`$ (its bare value is actually $`U6V`$, in terms of the parameters of the extended Hubbard model). We remark that, with a divergent density of states close to the Fermi level, the slightest attractive interaction ($`U_{𝐊_i}<0`$) triggers the CDW instability. The particle-hole susceptibility with momentum $`𝐊_i`$ is itself logarithmically divergent, and the iteration of particle-hole diagrams leads to an effective vertex at momentum transfer $`𝐊_i`$

$$\mathrm{\Gamma }_{𝐊_i}(\epsilon )=\frac{U_{𝐊_i}}{1+(c^{}/(2\pi ^2\mathrm{\Lambda }))U_{𝐊_i}\mathrm{log}|\mathrm{\Lambda }/\epsilon |}$$
(7)
Thus, the attractive interaction is overscreened and renormalized towards more intense attraction at low energies.

Under these conditions, the only response function that shows a divergence at some critical frequency corresponds to a charge modulation with wavevectors $`𝐊_i`$. The gap that opens up at the saddle points can be actually estimated from the frequency at which the denominator in Eq. (7) vanishes. The order of magnitude thus obtained (tens of meV) is consistent with the temperature of the transition to the CDW.

3. Metallic behavior and depletion of the density of states in the low-temperature phase.

The metallic behavior of the $`(3\times 3)`$ phase follows naturally from our description of the different subbands in the system. Nesting of the saddle points at the Van Hove filling and the formation of the surface CDW lead to the opening of a gap and to a partial destruction of the Fermi line around the saddle points. The divergent density of states induces by itself a significant reduction of the quasiparticle weight, as appreciated in the computation of the self-energy corrections for the saddle point dispersion relation. This is consistent with the decrease of the photoemission intensity about the Fermi level observed in the CDW phase. The lower band is less sensitive to these renormalization effects, and the metallic properties come mainly from the band crossing the Fermi level below the cusps shown in Fig. 1.

In conclusion, the pinning of the Fermi level to a Van Hove singularity in the density of states leads to a phenomenology for the Sn/Ge(111) overlayer consistent with the experimental observations of CDW formation and photoemission measurements, which are challenging more conventional theoretical approaches. Further predictions of our approach should also be tested experimentally in Sn/Ge(111), in particular the renormalization of the lower band and the geometry of the Fermi line around the corners of the $`(3\times 3)`$ Brillouin zone.